深入理解自适应滤波与回声消除

深入理解 LMS 算法：自适应滤波与回声消除

在信号处理领域，自适应滤波是一种重要的技术，广泛应用于噪声消除、回声消除和信号恢复等任务。各种自适应滤波算法的原理，包括公式推导，并通过 Python 代码示例展示其在回声消除中的应用。

收敛效果对比

注：本文参数没有详细调优，效果仅供参考

1. LMS 算法介绍

1.1 算法原理

LMS 算法的目标是通过最小化输出信号与目标信号之间的均方误差（Mean Squared Error, MSE）来调整滤波器的系数。我们可以定义以下信号：

• 参考信号 $far(n)$：这是我们希望消除的回声信号（例如，来自扬声器的原始信号）。
• 经过系统的信号 $mic(n)$：这是通过扬声器和麦克风系统接收到的信号，通常包含了回声和噪声。
• 估计信号 $y(n)$：这是自适应滤波器的输出信号，用于估计回声。
• 残余回声 $e(n)$：这是输出信号与目标信号之间的误差，表示未能消除的回声部分。

1.2 公式推导

LMS 算法的目标是最小化以下均方误差：
$$E=\mathbb{E}\{[mic(n)-y(n){]}^2\}$$
其中，( y(n) ) 是由自适应滤波器生成的输出信号，可以表示为：
$$y(n)=w^T(n)\cdot far(n)$$
这里 $w(n)$是滤波器的系数向量。
在每次迭代中，LMS 算法执行以下步骤：

1. 计算输出：
   $$y(n)=w^T(n)\cdot far(n)$$
2. 计算误差：
   $$e(n)=mic(n)-y(n)$$
3. 更新滤波器系数：

$$w(n+1)=w(n)+\mu \cdot e(n)\cdot far(n)$$

其中，$\mu$ 是学习率，控制每次更新的幅度。

1.3 优缺点

优点：

• 简单易实现，适合实时应用。
• 能够在线学习，适应信号的变化。
• 计算复杂度低，适合资源受限的环境。

缺点：

• 收敛速度可能较慢，尤其在高噪声环境下。
• 学习率选择不当可能导致不稳定。
• 可能收敛到局部最优解，而非全局最优解。

2. 更好的替代算法

除了 LMS 算法，还有许多其他自适应滤波算法，它们在某些情况下可能表现得更好。以下是一些常见的替代算法及其特点。
是的，除了 LMS（Least Mean Squares）算法，还有许多其他自适应滤波算法，它们在某些情况下可能表现得更好。以下是一些常见的替代算法及其特点：

2.1. NLMS（Normalized Least Mean Squares）算法

• 概述：NLMS 是 LMS 的一种改进版本，通过归一化输入信号的能量来调整学习率。这有助于提高算法的稳定性和收敛速度。
• 优点：
  • 更加稳定，尤其是在输入信号能量变化较大的情况下。
  • 收敛速度通常比 LMS 更快。
• 公式：
  $$w(n+1)=w(n)+\frac{\mu }{\Vert x(n){\Vert }^2}e(n)far(n)$$

2.2. RLS（Recursive Least Squares）算法

• 概述：RLS 是一种基于最小二乘原理的自适应滤波算法，通过递归更新滤波器系数来最小化误差平方和。
• 优点：
  • 收敛速度快，通常优于 LMS 和 NLMS。
  • 能够处理非平稳信号，适应性强。
• 缺点：
  • 计算复杂度较高，尤其是在滤波器阶数较大时。
  • 需要更多的内存和计算资源。
• 公式：
  • 更新公式：
    $$w(n+1)=w(n)+\frac{P(n)\cdot far(n)\cdot e(n)}{\lambda +far(n{)}^{T}P(n)\cdot far(n)}$$
  • 其中，$P(n)$ 是协方差矩阵，$\lambda$ 是遗忘因子（通常接近于 1），用于控制过去数据的影响。

2.3. Affined LMS（A-LMS）算法

• 概述：A-LMS 是对 LMS 的一种变体，结合了线性预测和自适应滤波的思想。
• 优点：
  • 可以更好地处理噪声和信号相位的变化。
  • 在某些应用中表现出更好的性能。
• 公式：
  • 更新公式：
    $$w(n+1)=w(n)+\mu \cdot e(n)\cdot (far(n)-\hat{y}(n))$$
  • 其中，$\hat{y}(n)$是基于当前权重的预测信号。

2.4. Sign LMS（S-LMS）算法

• 概述：S-LMS 是 LMS 的一种简化版本，它只使用符号信息（正负）来更新权重。
• 优点：
  • 计算复杂度低，适合实时应用。
  • 在某些情况下，能够提供与 LMS 相似的性能。
• 缺点：
  • 收敛速度较慢，且对噪声的鲁棒性较差。
• 公式：
  • 更新公式：
    $$w(n+1)=w(n)+\mu \cdot \text{sign}(e(n))\cdot far(n)$$
  • 其中，$\text{sign}(e(n))$是误差信号的符号函数，返回 +1 或 -1。

2.5. Adaptive Filter with Kalman Filter

• 概述：卡尔曼滤波器是一种基于状态空间模型的滤波方法，适用于动态系统的状态估计。
• 优点：
  • 能够处理动态变化的系统，适应性强。
  • 提供最优估计，尤其是在高噪声环境下表现良好。
• 缺点：
  • 数学推导复杂，计算资源消耗大。
• 公式：
  • 状态更新公式：
    $$\hat{x}(n|n)=\hat{x}(n|n-1)+K(n)\cdot (y(n)-H\cdot \hat{x}(n|n-1))$$
  • 其中，$K(n)$ 是卡尔曼增益，$H$是观测矩阵。
  • 卡尔曼增益的计算公式：
    $$K(n)=P(n|n-1)\cdot H^T\cdot (H\cdot P(n|n-1)\cdot H^T+R{)}^{-1}$$
  • 其中，$P(n|n-1)$是预测误差协方差矩阵，$R$是观测噪声协方差。

import numpy as np
import matplotlib.pyplot as plt


def plot_signals(mic, far, y, e, algorithm_name):
    # 创建子图
    fig, axs = plt.subplots(4, 1, figsize=(16, 9))

    # 绘制目标信号
    axs[0].plot(mic, label='Mic Signal', color='blue', linestyle='-', alpha=0.5)  # 目标信号
    axs[0].set_title(f'{algorithm_name}')
    axs[0].set_xlabel('Sample Index')
    axs[0].set_ylabel('Amplitude')
    axs[0].legend()
    axs[0].grid()

    # 绘制原始信号
    axs[1].plot(far, label='Far Signal', color='green', linestyle='--', alpha=0.5)  # 原始信号
    axs[1].set_xlabel('Sample Index')
    axs[1].set_ylabel('Amplitude')
    axs[1].legend()
    axs[1].grid()

    # 绘制输出信号
    axs[2].plot(y, label='Estimated Signal', color='red', linestyle='-')  # 输出信号
    axs[2].set_xlabel('Sample Index')
    axs[2].set_ylabel('Amplitude')
    axs[2].legend()
    axs[2].grid()

    # 绘制残余回声 e(n)
    axs[3].plot(e, label='Residue Echo Signal', color='orange', linestyle='-', alpha=0.7)  # 残余回声
    axs[3].set_xlabel('Sample Index')
    axs[3].set_ylabel('Amplitude')
    axs[3].legend()
    axs[3].grid()

    # 调整布局
    plt.tight_layout()
    plt.show()


def lms_algorithm(far, mic, mu=0.06, order=10):
    N = len(far)
    w = np.zeros(order)  # 初始化权重为零
    y = np.zeros(N)  # 输出信号
    e = np.zeros(N)  # 误差信号

    # LMS算法迭代
    for i in range(order, N):  # 从 order 开始迭代
        # 获取最近的 order 个输入样本
        input_samples = far[i - order:i]  # 当前输入样本
        # 计算输出
        y[i] = np.dot(w, input_samples)  # 使用权重和输入样本计算输出
        # 计算误差
        e[i] = mic[i] - y[i]  # 计算误差信号
        # 更新滤波器系数
        w += mu * e[i] * input_samples  # 更新公式
    return y, e


def nlms_algorithm(far, mic, mu=0.06, order=10):
    N = len(far)
    w = np.zeros(order)  # 初始化权重为零
    y = np.zeros(N)  # 输出信号
    e = np.zeros(N)  # 误差信号

    # NLMS算法迭代
    for i in range(order, N):  # 从 order 开始迭代
        # 获取最近的 order 个输入样本
        input_samples = far[i - order:i]  # 当前输入样本
        # 计算输出
        y[i] = np.dot(w, input_samples)  # 使用权重和输入样本计算输出
        # 计算误差
        e[i] = mic[i] - y[i]  # 计算误差信号
        # 计算输入样本的能量
        input_energy = np.dot(input_samples, input_samples)  # 输入信号的能量
        # 更新滤波器系数，避免除以零
        if input_energy > 1e-6:  # 防止除以零
            w += (mu / input_energy) * e[i] * input_samples  # 更新公式
    return y, e


def rls_algorithm(far, mic, lam=0.99, mu=0.06, order=10):
    N = len(far)
    w = np.zeros(order)  # 初始化权重为零
    y = np.zeros(N)  # 输出信号
    e = np.zeros(N)  # 误差信号
    P = np.eye(order) * 1e6  # 初始化协方差矩阵，较大的值以确保初始更新
    # RLS算法迭代
    for i in range(order, N):  # 从 order 开始迭代
        # 获取最近的 order 个输入样本
        input_samples = far[i - order:i]  # 当前输入样本
        # 计算输出
        y[i] = np.dot(w, input_samples)  # 使用权重和输入样本计算输出
        # 计算误差
        e[i] = mic[i] - y[i]  # 计算误差信号
        # 计算增益
        Pi = np.dot(P, input_samples)  # 协方差矩阵与输入样本的乘积
        k = Pi / (lam + np.dot(input_samples, Pi))  # 计算增益
        # 更新滤波器系数
        w += k * e[i]  # 更新公式
        # 更新协方差矩阵
        P = (P - np.outer(k, Pi)) / lam  # 更新协方差矩阵
    return y, e


def alms_algorithm(far, mic, mu=0.06, order=10):
    N = len(far)
    w = np.zeros(order)  # 初始化权重为零
    y = np.zeros(N)  # 输出信号
    e = np.zeros(N)  # 误差信号
    # A-LMS算法迭代
    for i in range(order, N):  # 从 order 开始迭代
        # 获取最近的 order 个输入样本
        input_samples = far[i - order:i]  # 当前输入样本
        # 计算预测信号
        y[i] = np.dot(w, input_samples)  # 使用权重和输入样本计算预测信号
        # 计算误差
        e[i] = mic[i] - y[i]  # 计算误差信号
        # 更新滤波器系数
        w += mu * e[i] * (input_samples - y[i])  # 更新公式
    return y, e


def slms_algorithm(far, mic, mu=0.06, order=10):
    N = len(far)
    w = np.zeros(order)  # 初始化权重为零
    y = np.zeros(N)  # 输出信号
    e = np.zeros(N)  # 误差信号

    # S-LMS算法迭代
    for i in range(order, N):  # 从 order 开始迭代
        # 获取最近的 order 个输入样本
        input_samples = far[i - order:i]  # 当前输入样本
        # 计算输出
        y[i] = np.dot(w, input_samples)  # 使用权重和输入样本计算输出
        # 计算误差
        e[i] = mic[i] - y[i]  # 计算误差信号
        # 更新滤波器系数
        w += mu * np.sign(e[i]) * input_samples  # 更新公式
    return y, e


# 参数设置
N = 800  # 迭代次数
frequencies = [0.1, 0.07, 0.18]  # 正弦波频率列表
np.random.seed(0)  # 设置随机种子以便重现

# 生成原始参考信号（复杂信号）
n = np.arange(N)
far = sum(0.2 * np.sin(2 * np.pi * f * n) for f in frequencies)  # 复杂信号
mic = far * 2  # 目标信号

# 调用LMS算法
y_lms, e_lms = lms_algorithm(far, mic, mu=0.05, order=16)
# 绘制LMS信号
plot_signals(mic, far, y_lms, e_lms, algorithm_name='LMS')
# 调用NLMS算法
y_nlms, e_nlms = nlms_algorithm(far, mic, mu=0.05, order=16)
# 绘制NLMS信号
plot_signals(mic, far, y_nlms, e_nlms, algorithm_name='NLMS')
# 调用RLS算法
y_rls, e_rls = rls_algorithm(far, mic, lam=0.99, order=16)
# 绘制RLS信号
plot_signals(mic, far, y_rls, e_rls, algorithm_name='RLS')
# 调用A-LMS算法
y_alms, e_alms = alms_algorithm(far, mic, mu=0.05, order=16)
# 绘制A-LMS信号
plot_signals(mic, far, y_alms, e_alms, algorithm_name='A-LMS')
# 调用S-LMS算法
y_slms, e_slms = slms_algorithm(far, mic, mu=0.05, order=16)
# 绘制S-LMS信号
plot_signals(mic, far, y_slms, e_slms, algorithm_name='S-LMS')