自适应滤波器更新算法-EP1

  自适应滤波器是回声消除系统中非常重要的一个功能模块，而对于自适应滤波器来说，如果更新滤波器系数则是关键所在。本文将介绍几种现有的滤波器更新算法，并附上Matlab测试代码。

1、LMS算法

1.1算法原理

  LMS算法即最小均方算法（Least Mean Square），它是由Widrow和Hoff在1960年提出的。将LMS算法应用在滤波器上的结构中，其包含两个主要的部分，第一个部分是LMS算法的部分，另一个部分是常用的滤波器结构即横向有限冲击响应结构，假设ｎ时刻的输入信号向量表达式为
$$X(n)=[x(n),x(n-1),...,x(n-M+1)]$$
其中，M表示滤波器的阶数，是一个固定值，其根据实际需求而定。滤波器的权向量的表达式为

其中，i是滤波器权向量的下标，即表示 $w_i(n)$ 为第几个滤波器权向量，其取值范围为0~M-1，如 $w_0(n)$ 表示的是 $x(n)$ 所对应的权向量。滤波器输出信号表达式为：

n时刻期望信号 $d(n)$ 与算法实际输出信号之间的差值表示为：

其中误差信号 $e(n)$ 进入LMS算法的函数模型中，通过该函数模型控制着LMS算法的步长参数，随着每一次迭代系数向量 $W(n)$ 发生改变，最终 $W(n)$ 趋近于其最有解，从而使实际输出信号 $y(n)$ 与期望信号 $d(n)$ 的差值无限接近与0。

  传统LMS算法的具体步骤如下：

1. 初始化参数。将自适应滤波器权系数的初始值 $W(0)$ 设置为0，将自适应滤波器的阶数设置为M，选择合适的变化的步长值 $\mu$ ，为了使算法收敛，步长 $\mu$ 的取值应为：
   $$0<\mu <1/{\lambda }_{max}$$
   其中 ${\lambda }_{max}$ 表示输入信号 $X(n)$ 的自相关矩阵的最大特征值。

2. 计算滤波器的实际输出信号为 ：$y(n)=W^T(n)X(n)$

3. 计算自适应滤波器的误差为 ： $e(n)=d(n)-y(n)=d(n)-W^T(n)X(n)$

4. 更新滤波器权系数的值为： $W(n+1)=W(n)+\mu X(n)e(n)$

重复上述步骤2~4，直到滤波器权系数逼近最佳滤波为止。

1.2 算法代码

function [e, y, w] = LMS(d, x, mu, M)
% Inputs:
% d  - 麦克风语音
% x  - 远端语音
% mu - 步长，0.05
% M  - 滤波器阶数，也称为抽头数
%
% Outputs:
% e - 输出误差，近端语音估计
% y  - 输出系数，远端回声估计
% w - 滤波器参数d_length = length(d);if (d_length <= M)  print('error: 信号长度小于滤波器阶数！');return; endif (d_length ~= length(x))  print('error: 输入信号和参考信号长度不同！');return; endxx = zeros(M,1);w1 = zeros(M,1);    % 滤波器权重y = zeros(d_length,1);  % 远端回声估计e = zeros(d_length,1);  % 近端语音估计for n = M:d_lengthxx = x(n:-1:n-M+1);    % 纵向拼接  (40~1)-->(41~2)-->(42~3)....y(n) = w1' * xx;        % 远端回声估计 (40,1)'*(40,1)=1; (73113,1)e(n) = d(n) - y(n);     % 近端语音估计w1 = w1 + mu * e(n) * xx;   % (40,1)w(:,n) = w1;        % (40, 73113)end
end

1.3 算法分析

  LMS算法是一种随机梯度下降算法，随着迭代次数的增加，更新出的权重值无法精确收敛，但是其算法简单，易于实现，算法复杂度也较低。但是算法中的 $\mu$ 为步长因子，其对于算法的收敛速度、时变系统的跟踪速度和稳态误差这些特征起到关键的影响作用，较小的步长因子 $\mu$ 会获得较小的稳态误差从而提高算法的精度，同时也会获得较小的收敛速度，使算法收敛的很慢，对于时变系统来说较小的步长因子也会导致跟踪速度的减小。反之，较大的步长因子 $\mu$ 会带来较大的收敛速度使得算法收敛的很快，同时对于时变系统来说，会提高其跟踪速度，但是会带来较大的稳态误差，降低算法的精度。因此 $\mu$ 的选取很关键。

2、NLMS算法

2.1算法原理

  归一化最小均方(Normalized Least Mean Squares，NLMS)算法是改进的变步长的LMS算法，对于较大的输入 $X(n)$ ，会导致梯度噪声的放大，因此需要用输入向量的平方范数进行归一化。根据原LMS算法中误差信号与远端输入信号的乘积，对远端输入信号的平方(功率)进行归一化处理，将固定步长因子的LMS算法变为根据输入信号时变的变步长LMS算法，算法具体针对 $\mu$ 改进为函数 $\mu (n)$ 如下：
$$\mu (n)=\frac{\mu }{X^T(n)X(n)+\delta }$$
其中，$\mu$ 为自适应常量，取值为0~2，$\delta$ 为一个较小的正常数，是为了防止输入信号 $X(n)$ 过小导致分母趋近于0进而带来的算法不稳定性。则滤波器权系数更新公式为：
$$W(n+1)=W(n)+\mu (n)X(n)e(n)$$
  具体的实现步骤，除了$\mu (n)$ 的计算外，其他同LMS算法步骤相同。

2.2 算法代码

function[e,y,w] = NLMS(d,x,M)% 输入:% d  - 麦克风语音% x  - 远端语音% M  - 滤波器阶数%% 输出:% e - 近端语音估计% y - 远端回声估计% w - 滤波器参数d_length = length(d);if (d_length <= M)print('error: 信号长度小于滤波器阶数！');return; endif (d_length ~= length(x))  print('error: 输入信号和参考信号长度不同！');return; endxx = zeros(M,1);w1 = zeros(M,1);  %滤波器权重y = zeros(d_length,1); %近端语音e = zeros(d_length,1); %误差for n = 1:d_length%在输入信号x前补上M-1个0，使输出y与输入具有相同长度xx = [xx(2:M);x(n)]; %(39,1) +(1,1) = (40,1)y(n) = w1' *xx;mu = 1/(xx'*xx); %步长,这里简化了mu的计算，没有防止x（n）过小的情况e(n) = d(n) - y(n);  %误差w1 = w1 + mu * e(n) * xx;w(:,n) = w1;end
end

1.3 算法分析

  NLMS算法相对与LMS算法有着更快的收敛速度，同时保留了算法简洁、计算量小、易于实现的优点；但是NLMS算法没有考虑到滤波器阶数增长的问题。随着采样率的提升，滤波器的阶数已经达到千阶，算法复杂度也成倍增长。因此NLMS 算法仅适用于对处理精度和能力要求不是很高的一般场合。

3、PBFDAF算法

3.1 算法原理

参考文献一中对该算法的介绍：

3.2 算法代码

function [en, yk, W] = myFDAF(d,x,mu,mu_unconst, M, select)
% 参数:
% d                输入信号(麦克风语音)
% x                远端语音
% mu                约束 FDAF的步长
% mu_unconst        不受约束的 FDAF的步长
% M                 滤波器阶数
% select;            选择有约束或无约束FDAF算法x_new = zeros(M,1);     % 将新块的M个样本初始化为0
x_old = zeros(M,1);     % 将旧块的M个样本初始化为0AdaptStart = 2*M;       % 在获得2M个样本块后开始自适应
W = zeros(2*M,1);       % 将2M个滤波器权重初始化为0
d_block = zeros(M,1);   % 将期望信号的M个样本初始化为0power_alpha = 0.5;        % 常数以更新每个frequency bin的功率
power = zeros(2*M,1);   % 将每个bin的平均功率初始化为0
d_length = length(d);             % 输入序列的长度
en = [];                       % 误差向量
window_save_first_M = [ones(1,M), zeros(1,M)]';  % 设置向量以提取前M个元素 (2M,1)for k = 1:d_lengthx_new = [x_new(2:end,:); x(k)];         % 开始的输入信号块 (2M,1)d_block = [d_block(2:end,:); d(k)];     % 开始的期望信号快 (M,1)if mod(k,M)==0                        % If iteration == block length, x_blocks = [x_old; x_new];        % 2M样本的输入信号样本块 (2M,1)x_old = x_new;if k >= AdaptStart                % 频域自适应滤波器% 将参考信号转换到频域Xk = fft(x_blocks);     % (2M,1)% FFT[old block; new block]% Old block 包含M个先前的输入样本 (u_old)% New 包含M个新的输入样本 (u_new)% 计算滤波器估计信号Yk = Xk.*W;                  % 输入和权重向量的乘积(2M,1)*(2M,1)=(2M,1)temp = real(ifft(Yk));            % IFFT 输出的实部 (2M,1)yk = temp(M+1:2*M);               % 抛弃前M个元素，保留后M个元素 (M,1)% 计算误差信号error = d_block-yk;              % 误差信号块 (M,1)Ek = fft([zeros(1,M),error']');   % 在FFT之前插入零块以形成2M块(2M,1)% 更新信号功率估算power = (power_alpha.*power) + (1 - power_alpha).*(abs(Xk).^2); % (2M,1)% 计算频域中的梯度和权重更新if select == 1gradient = real(ifft((1./power).*conj(Xk).* Ek));   % (2M,1)gradient = gradient.*window_save_first_M;   % 去除后一个数据块，并且补零 (2M,1)W = W + mu.*(fft(gradient));    % 权重是频域的 (2M,1)elsegradient = conj(Xk).* Ek;   %  (2M,1)W = W + mu_unconst.*gradient;    % (2M,1)enden = [en; error];             % 更新误差块endend
end% 参考：https://github.com/CharlesThaCat/acoustic-interference-cancellation    

3.3 算法分析

  该算法利用FFT可快速实现卷积运算的特点，在频域进行滤波器的更新，同时使用逐块更新来代替逐点更新，从而大大降低了算法的计算复杂度。但同时其步长因子 $\mu$ 是由自适应常量和远端信号决定的，是固定的无法根据不同复杂的场景进行变换，该算法应用在WebRTC的AEC算法中。在AEC算法中，其用了延时估计模块来对齐语音信号，并采用了语音检测等模块来去除噪声和无用语音的影响，以提高滤波器更新算法的效果。

续：本文先介绍几种比较经典的自适应滤波器算法，后续将介绍读者阅读论文总结并复现的滤波器算法！！！

参考文献：

[1] 葛世超，吕强。《实时语音处理实践指南》。

[2] Soo J S , Pang K K . Multidelay block frequency domain adaptive filter[J]. Acoustics Speech & Signal Processing IEEE Transactions on, 1990, 38(2):373-376.