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

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

在信号处理领域，自适应滤波是一种重要的技术，广泛应用于噪声消除、回声消除和信号恢复等任务。LMS（Least Mean Squares）算法是实现自适应滤波的经典方法之一。本文将详细介绍 LMS 算法的原理，包括公式推导，并通过 Python 代码示例展示其在回声消除中的应用。我们还将介绍一些替代算法，比较它们的收敛效果。

收敛效果对比

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

1. LMS 算法介绍

1.1 算法原理

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

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

1.2 公式推导

LMS 算法的目标是最小化以下均方误差：

$$E=\mathbb{E}\{[d(n)-y(n){]}^2\}$$

其中，( y(n) ) 是由自适应滤波器生成的输出信号，可以表示为：

$$y(n)=w^T(n)\cdot x(n)$$

这里 ( w(n) ) 是滤波器的系数向量。

在每次迭代中，LMS 算法执行以下步骤：

1. 计算输出：

$$y(n)=w^T(n-1)\cdot x(n-1)$$

2. 计算误差：

$$e(n)=d(n)-y(n)$$

3. 更新滤波器系数：

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

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

1.3 优缺点

优点：

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

缺点：

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

1.4. LMS代码介绍

下面是使用 LMS 算法的 Python 示例，展示了如何通过输入信号和目标信号进行自适应滤波。

import numpy as np
import matplotlib.pyplot as plt# 参数设置
mu = 0.06  # 学习率
N = 300  # 迭代次数
f1 = 0.1  # 正弦波频率1
f2 = 0.07  # 正弦波频率2
np.random.seed(0)  # 设置随机种子以便重现# 生成原始参考信号（复杂信号）
n = np.arange(N)
x = 0.5 * np.sin(2 * np.pi * f1 * n) + 0.5 * np.sin(2 * np.pi * f2 * n)  # 复杂信号# 定义传递函数的参数
a = 1.0  # 基础传递函数增益
variation = 0.2  # 传递函数的微小变动范围# 生成经过传递函数的信号，添加随机变化
transfer_function = a + np.random.uniform(-variation, variation, N)  # 随机变化
d = x * transfer_function  # 目标信号# 初始化滤波器系数
order = 10  # 滤波器阶数
w = np.zeros(order)  # 初始化权重为零
y = np.zeros(N)  # 输出信号
e = np.zeros(N)  # 误差信号# LMS算法迭代
for i in range(order, N):  # 从 order 开始迭代# 获取最近的 order 个输入样本input_samples = x[i-order:i]  # 当前输入样本# 计算输出y[i] = np.dot(w, input_samples)  # 使用权重和输入样本计算输出# 计算误差e[i] = d[i] - y[i]  # 计算误差信号# 更新滤波器系数w += mu * e[i] * input_samples  # 更新公式# 创建子图
fig, axs = plt.subplots(2, 1, figsize=(12, 10))# 绘制目标信号、原始信号和输出信号
axs[0].plot(d, label='Received Signal (d)', color='blue', linestyle='-', alpha=0.5)  # 目标信号
axs[0].plot(x, label='Original Signal (x)', color='green', linestyle='--', alpha=0.5)  # 原始信号
axs[0].plot(y, label='Estimated Signal (y)', color='red', linestyle='-')  # 输出信号
axs[0].set_title('LMS Algorithm with Transfer Function Example (Order 8)')
axs[0].set_xlabel('Sample Index')
axs[0].set_ylabel('Amplitude')
axs[0].legend()
axs[0].grid()# 绘制误差信号
axs[1].plot(e * e, label='Error Signal rms(e)', color='orange', linestyle='-', alpha=0.7)  # 误差信号
axs[1].set_title('Error Signal After Echo Cancellation')
axs[1].set_xlabel('Sample Index')
axs[1].set_ylabel('Amplitude')
axs[1].legend()
axs[1].grid()# 调整布局
plt.tight_layout()
plt.show()

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)x(n)$$

2.2. RLS（Recursive Least Squares）算法

• 概述：RLS 是一种基于最小二乘原理的自适应滤波算法，通过递归更新滤波器系数来最小化误差平方和。
• 优点：
  • 收敛速度快，通常优于 LMS 和 NLMS。
  • 能够处理非平稳信号，适应性强。
• 缺点：
  • 计算复杂度较高，尤其是在滤波器阶数较大时。
  • 需要更多的内存和计算资源。
• 公式：
  • 更新公式较为复杂，涉及协方差矩阵的计算。

2.3. Affined LMS（A-LMS）算法

• 概述：A-LMS 是对 LMS 的一种变体，结合了线性预测和自适应滤波的思想。
• 优点：
  • 可以更好地处理噪声和信号相位的变化。
  • 在某些应用中表现出更好的性能。

2.4. Sign LMS（S-LMS）算法

• 概述：S-LMS 是 LMS 的一种简化版本，它只使用符号信息（正负）来更新权重。
• 优点：
  • 计算复杂度低，适合实时应用。
  • 在某些情况下，能够提供与 LMS 相似的性能。
• 缺点：
  • 收敛速度较慢，且对噪声的鲁棒性较差。

2.5. Adaptive Filter with Kalman Filter

• 概述：卡尔曼滤波器是一种基于状态空间模型的滤波方法，适用于动态系统的状态估计。
• 优点：
  • 能够处理动态变化的系统，适应性强。
  • 提供最优估计，尤其是在高噪声环境下表现良好。
• 缺点：
  • 数学推导复杂，计算资源消耗大。

2.6. 代码实现

1. NLMS（Normalized Least Mean Squares）

import numpy as np
import matplotlib.pyplot as plt# 参数设置
initial_mu = 0.06  # 初始学习率
N = 300  # 迭代次数
f1 = 0.1  # 正弦波频率1
f2 = 0.07  # 正弦波频率2
np.random.seed(0)  # 设置随机种子以便重现# 生成原始参考信号（复杂信号）
n = np.arange(N)
x = 0.5 * np.sin(2 * np.pi * f1 * n) + 0.5 * np.sin(2 * np.pi * f2 * n)  # 复杂信号# 定义传递函数的参数
a = 1.0  # 基础传递函数增益
variation = 0.2  # 传递函数的微小变动范围# 生成经过传递函数的信号，添加随机变化
transfer_function = a + np.random.uniform(-variation, variation, N)  # 随机变化
d = x * transfer_function  # 目标信号# 初始化滤波器系数
order = 10  # 滤波器阶数
w = np.zeros(order)  # 初始化权重为零
y = np.zeros(N)  # 输出信号
e = np.zeros(N)  # 误差信号# NLMS算法迭代
for i in range(order, N):  # 从 order 开始迭代# 获取最近的 order 个输入样本input_samples = x[i-order:i]  # 当前输入样本# 计算输出y[i] = np.dot(w, input_samples)  # 使用权重和输入样本计算输出# 计算误差e[i] = d[i] - y[i]  # 计算误差信号# 动态学习率mu = initial_mu / (1 + 0.01 * np.dot(input_samples, input_samples))  # 归一化学习率# 更新滤波器系数w += mu * e[i] * input_samples  # 更新公式# 创建子图
fig, axs = plt.subplots(2, 1, figsize=(12, 10))# 绘制目标信号、原始信号和输出信号
axs[0].plot(d, label='Received Signal (d)', color='blue', linestyle='-', alpha=0.5)  # 目标信号
axs[0].plot(x, label='Original Signal (x)', color='green', linestyle='--', alpha=0.5)  # 原始信号
axs[0].plot(y, label='Estimated Signal (y)', color='red', linestyle='-')  # 输出信号
axs[0].set_title('NLMS Algorithm with Transfer Function Example (Order 10)')
axs[0].set_xlabel('Sample Index')
axs[0].set_ylabel('Amplitude')
axs[0].legend()
axs[0].grid()# 绘制误差信号
axs[1].plot(e**2, label='Error Signal (rms[e])', color='orange', linestyle='-', alpha=0.7)  # 误差信号
axs[1].set_title('Error Signal After Echo Cancellation (NLMS)')
axs[1].set_xlabel('Sample Index')
axs[1].set_ylabel('Amplitude')
axs[1].legend()
axs[1].grid()# 调整布局
plt.tight_layout()
plt.show()

2. RLS（Recursive Least Squares）

# RLS算法实现
import numpy as np
import matplotlib.pyplot as plt# 参数设置
N = 300  # 迭代次数
f1 = 0.1  # 正弦波频率1
f2 = 0.07  # 正弦波频率2
np.random.seed(0)  # 设置随机种子以便重现# 生成原始参考信号（复杂信号）
n = np.arange(N)
x = 0.5 * np.sin(2 * np.pi * f1 * n) + 0.5 * np.sin(2 * np.pi * f2 * n)  # 复杂信号# 定义传递函数的参数
a = 1.0  # 基础传递函数增益
variation = 0.2  # 传递函数的微小变动范围# 生成经过传递函数的信号，添加随机变化
transfer_function = a + np.random.uniform(-variation, variation, N)  # 随机变化
d = x * transfer_function  # 目标信号# 初始化滤波器系数
M = 10 # 滤波器阶数
w = np.zeros(M)  # # 初始化权重为零
y = np.zeros(N)  # 输出信号
e = np.zeros(N)  # 误差信号
P = np.eye(M) * 1000  # 初始协方差矩阵# RLS算法迭代
for i in range(M, N):x_vec = x[i-M:i][::-1]  # 获取当前输入信号的向量y[i] = np.dot(w, x_vec)  # 计算输出信号e[i] = d[i] - y[i]  # 计算误差信号k = P @ x_vec / (1 + np.dot(x_vec.T, P @ x_vec))  # 计算增益w += k * e[i]  # 更新权重P = (P - np.outer(k, k.T) * (1 + np.dot(x_vec.T, k)))  # 更新协方差矩阵# 创建子图
fig, axs = plt.subplots(2, 1, figsize=(12, 10))# 绘制目标信号、原始信号和输出信号
axs[0].plot(d, label='Received Signal (d)', color='blue', linestyle='-', alpha=0.5)  # 目标信号
axs[0].plot(x, label='Original Signal (x)', color='green', linestyle='--', alpha=0.5)  # 原始信号
axs[0].plot(y, label='Estimated Signal (y)', color='red', linestyle='-')  # 输出信号
axs[0].set_title('RLS Algorithm with Transfer Function Example')
axs[0].set_xlabel('Sample Index')
axs[0].set_ylabel('Amplitude')
axs[0].legend()
axs[0].grid()# 绘制误差信号
axs[1].plot(e**2, label='Error Signal (rms[e])', color='orange', linestyle='-', alpha=0.7)  # 误差信号
axs[1].set_title('Error Signal After Echo Cancellation (RLS)')
axs[1].set_xlabel('Sample Index')
axs[1].set_ylabel('Amplitude')
axs[1].legend()
axs[1].grid()# 调整布局
plt.tight_layout()
plt.show()

3. Affined LMS

import numpy as np
import matplotlib.pyplot as plt# 参数设置
initial_mu = 0.06  # 初始学习率
N = 300  # 迭代次数
f1 = 0.1  # 正弦波频率1
f2 = 0.07  # 正弦波频率2
np.random.seed(0)  # 设置随机种子以便重现# 生成原始参考信号（复杂信号）
n = np.arange(N)
x = 0.5 * np.sin(2 * np.pi * f1 * n) + 0.5 * np.sin(2 * np.pi * f2 * n)  # 复杂信号# 定义传递函数的参数
a = 1.0  # 基础传递函数增益
variation = 0.2  # 传递函数的微小变动范围# 生成经过传递函数的信号，添加随机变化
transfer_function = a + np.random.uniform(-variation, variation, N)  # 随机变化
d = x * transfer_function  # 目标信号# 初始化滤波器系数
order = 10  # 滤波器阶数
w = np.zeros(order)  # 初始化权重为零
y = np.zeros(N)  # 输出信号
e = np.zeros(N)  # 误差信号# Affined LMS算法迭代
for i in range(order, N):  # 从 order 开始迭代# 获取最近的 order 个输入样本input_samples = x[i-order:i]  # 当前输入样本# 计算输出，增加偏置项y[i] = np.dot(w, input_samples) + 0.5 * input_samples[-1]  # 使用权重和输入样本计算输出，增加偏置# 计算误差e[i] = d[i] - y[i]  # 计算误差信号# 更新滤波器系数w += initial_mu * e[i] * input_samples  # 更新公式# 创建子图
fig, axs = plt.subplots(2, 1, figsize=(12, 10))# 绘制目标信号、原始信号和输出信号
axs[0].plot(d, label='Received Signal (d)', color='blue', linestyle='-', alpha=0.5)  # 目标信号
axs[0].plot(x, label='Original Signal (x)', color='green', linestyle='--', alpha=0.5)  # 原始信号
axs[0].plot(y, label='Estimated Signal (y)', color='red', linestyle='-')  # 输出信号
axs[0].set_title('Affined LMS Algorithm with Transfer Function Example (Order 10)')
axs[0].set_xlabel('Sample Index')
axs[0].set_ylabel('Amplitude')
axs[0].legend()
axs[0].grid()# 绘制误差信号
axs[1].plot(e**2, label='Error Signal (rms[e])', color='orange', linestyle='-', alpha=0.7)  # 误差信号
axs[1].set_title('Error Signal After Echo Cancellation (Affined LMS)')
axs[1].set_xlabel('Sample Index')
axs[1].set_ylabel('Amplitude')
axs[1].legend()
axs[1].grid()# 调整布局
plt.tight_layout()
plt.show()

4. Sign LMS（S-LMS）

import numpy as np
import matplotlib.pyplot as plt# 参数设置
initial_mu = 0.06  # 初始学习率
N = 300  # 迭代次数
f1 = 0.1  # 正弦波频率1
f2 = 0.07  # 正弦波频率2
np.random.seed(0)  # 设置随机种子以便重现# 生成原始参考信号（复杂信号）
n = np.arange(N)
x = 0.5 * np.sin(2 * np.pi * f1 * n) + 0.5 * np.sin(2 * np.pi * f2 * n)  # 复杂信号# 定义传递函数的参数
a = 1.0  # 基础传递函数增益
variation = 0.2  # 传递函数的微小变动范围# 生成经过传递函数的信号，添加随机变化
transfer_function = a + np.random.uniform(-variation, variation, N)  # 随机变化
d = x * transfer_function  # 目标信号# 初始化滤波器系数
order = 10  # 滤波器阶数
w = np.zeros(order)  # 初始化权重为零
y = np.zeros(N)  # 输出信号
e = np.zeros(N)  # 误差信号# Sign LMS算法迭代
for i in range(order, N):  # 从 order 开始迭代# 获取最近的 order 个输入样本input_samples = x[i-order:i]  # 当前输入样本# 计算输出y[i] = np.dot(w, input_samples)  # 使用权重和输入样本计算输出# 计算误差e[i] = d[i] - y[i]  # 计算误差信号# 更新滤波器系数（只使用符号信息）w += initial_mu * np.sign(e[i]) * input_samples  # 更新公式# 创建子图
fig, axs = plt.subplots(2, 1, figsize=(12, 10))# 绘制目标信号、原始信号和输出信号
axs[0].plot(d, label='Received Signal (d)', color='blue', linestyle='-', alpha=0.5)  # 目标信号
axs[0].plot(x, label='Original Signal (x)', color='green', linestyle='--', alpha=0.5)  # 原始信号
axs[0].plot(y, label='Estimated Signal (y)', color='red', linestyle='-')  # 输出信号
axs[0].set_title('Sign LMS Algorithm with Transfer Function Example (Order 10)')
axs[0].set_xlabel('Sample Index')
axs[0].set_ylabel('Amplitude')
axs[0].legend()
axs[0].grid()# 绘制误差信号
axs[1].plot(e**2, label='Error Signal (rms[e])', color='orange', linestyle='-', alpha=0.7)  # 误差信号
axs[1].set_title('Error Signal After Echo Cancellation (Sign LMS)')
axs[1].set_xlabel('Sample Index')
axs[1].set_ylabel('Amplitude')
axs[1].legend()
axs[1].grid()# 调整布局
plt.tight_layout()
plt.show()

5. 卡尔曼滤波器

import numpy as np
import matplotlib.pyplot as plt# 参数设置
N = 300  # 迭代次数
f1 = 0.1  # 正弦波频率1
f2 = 0.07  # 正弦波频率2
np.random.seed(0)  # 设置随机种子以便重现# 生成原始参考信号（复杂信号）
n = np.arange(N)
x = 0.5 * np.sin(2 * np.pi * f1 * n) + 0.5 * np.sin(2 * np.pi * f2 * n)  # 复杂信号# 定义传递函数的参数
a = 1.0  # 基础传递函数增益
variation = 0.2  # 传递函数的微小变动范围# 生成经过传递函数的信号，添加随机变化
transfer_function = a + np.random.uniform(-variation, variation, N)  # 随机变化
d = x * transfer_function  # 目标信号# 初始化卡尔曼滤波器参数
Q = 1e-6  # 过程噪声协方差
R = 0.2   # 测量噪声协方差
P = 0.2  # 估计误差协方差
M = 10  # 滤波器阶数
w = np.zeros(M)  # 初始化权重为零
y = np.zeros(N)  # 输出信号
e = np.zeros(N)  # 误差信号# 卡尔曼滤波器迭代
for i in range(M, N):  # 从 M 开始迭代# 获取最近的 M 个输入样本x_vec = x[i-M:i]  # 当前输入样本# 预测更新y[i] = np.dot(w, x_vec)  # 预测输出# 更新卡尔曼增益K = P / (P + R)  # 卡尔曼增益# 计算误差e[i] = d[i] - y[i]  # 计算误差信号# 更新估计w += K * e[i] * x_vec  # 更新权重# 更新估计误差协方差P = (1 - K) * P + Q  # 更新协方差# 创建子图
fig, axs = plt.subplots(2, 1, figsize=(12, 10))# 绘制目标信号、原始信号和输出信号
axs[0].plot(d, label='Received Signal (d)', color='blue', linestyle='-', alpha=0.5)  # 目标信号
axs[0].plot(x, label='Original Signal (x)', color='green', linestyle='--', alpha=0.5)  # 原始信号
axs[0].plot(y, label='Estimated Signal (y)', color='red', linestyle='-')  # 输出信号
axs[0].set_title('Kalman Filter Algorithm with Transfer Function Example (Order 10)')
axs[0].set_xlabel('Sample Index')
axs[0].set_ylabel('Amplitude')
axs[0].legend()
axs[0].grid()# 绘制误差信号
axs[1].plot(e**2, label='Error Signal (rms[e])', color='orange', linestyle='-', alpha=0.7)  # 误差信号
axs[1].set_title('Error Signal After Echo Cancellation (Kalman Filter)')
axs[1].set_xlabel('Sample Index')
axs[1].set_ylabel('Amplitude')
axs[1].legend()
axs[1].grid()# 调整布局
plt.tight_layout()
plt.show()