【噪声、干扰抑制】自适应滤波算法合集【附MATLAB代码】

1 自适应滤波器的基本概念

在实际应用中，常需要跟踪时变的、未知的信号衰减。自适应滤波器通过自 适应的改变滤波器系数对衰减信号进行补偿，从而达到降低噪声和失真的目的。针对不同的应用情况，有不同的自适应滤波器结构和自适应算法。本章讨论的自 适应滤波器基于FIR滤波器，并且利用了LMS(Least mean square，最小均方误差) 算法及其改进算法。 自适应滤波器由两个独立部分组成：自适应滤波器对期望信号进行处理，自 适应算法调整滤波器系数。自适应滤波器的常见模型如图1所示，d(n)是期望 信号（或基本的输入信号），y(n)是滤波器的输出信号，由参考输入信号x(n)驱动。误差信号e(n)是d(n)与y(n)的差值。自适应算法用于调整自适应滤波器系数，使 得e(n)的均方值最小。因此在一个接一个采样的基础上，滤波器系数得到更新， e(n)逐渐减小[36]。

图1 自适应滤波器

图2 FIR滤波器

一般来说有两种结构的自适应滤波器：IIR滤波器和FIR滤波器。IIR存在稳 定性问题，而FIR总是稳定的。在实时处理中FIR滤波器应用更为广泛[ 35]。最常 用的FIR自适应滤波器框图如2-2所示，给定长度为L的滤波器系数数组wl(n), l=0,1,…,L-1和数据序列{x(n)x(n-1)…x(n-L+1)}，滤波器输出信号如下式所示：

滤波器系数wl(n)是时变的，且通过自适应算法更新。 定义时刻n的输入数组如式：

时刻n的滤波器系数数组如式：

输出信号y(n)与期望响应d(n)相比较得出误差信号

自适应滤波器通过更新其系数，使某个预定的性能指标达到最佳，这个预定 的性能指标也称为代价函数[6]。应用最为广泛的自适应滤波器是维纳滤波器，其性 能函数是e(n)的均方误差(MSE)，定义为

最佳滤波器w0(n)使上式达到最小，维纳-霍夫方程给出了w0(n)的解法

其中，Rxx(n)和rxy(n)分别是信号x(n)的自相关矩阵以及x(n)与y(n)的互相关 矩阵，定义如下

求解上式可得滤波器系数的最佳解。但是该式的运算量很大，尤其是当自 适应滤波器阶数很高的情况下。为了避免直接计算自相关矩阵Rxx(n)， 可以采用基于梯度的最陡下降法递归地求得w0(n)。著名的LMS(least mean square， 最小均方误差)算法就是基于最陡下降法，下一节中将介绍LMS算法的原理。

2 自适应滤波算法

2.1 LMS算法

LMS(least mean square，最小均方误差)算法是应用最为广泛的随机梯度算法， 它较为简单且无需求出Rxx(n)和rxy(n)。最陡下降法的自适应滤波器系数更新公式 [36]为：

μ被称为步长因子，为一常数，用来控制作用 于滤波器系数矢量上的增量修正大小。在自适应滤波器系数迭代的过程中，μ的值必须满足下式以保证滤波器权系数矢量收敛：

λmax为x(n)的自相关矩阵Rxx的最大特征值。LMS算法是基于最陡下降法的一种自适应算法，根据最陡下降法进行合理近 似得到。LMS算法的导出过程为

即用误差信号平方值[e2(n)]的估计值代替性能函数，即误差信号平方值的统 计期望E{e2(n)}的估计值。因此可以得出

所以LMS算法的递推公式为：

MS算法的计算过程如下：

1. 确定L,μ和w(0)，L是滤波器阶数，μ是步长，w(0)是采样点n=0时刻滤波器系数的初始值，一般为零向量。

2. 计算滤波器输出

3. 计算误差信号

4. 利用LMS算法，从w(n)到w(n+1)更新自适应滤波器系数

LMS算法有很多优点，当输入信号是平稳信号时，算法能够迅速收敛到最优 解；鲁棒性较强，有限字长效应不会影响算法的稳定性；计算复杂度比较低，L阶 自适应滤波器只需2L+1次乘法和2L次加法。但是LMS算法本身也存在如下缺点： 第一，LMS算法的步长因子μ影响着算法的收敛速度和稳态误差[37]，比较小的μ会获得比较好的稳态性能，但收敛速度并不能得到保证；比较大的μ会获得比较快的收敛速度，但稳态性能无法保证。实际上，最佳步长μ是很难确定的， 因此不恰当的可能导致收敛速度不必要的下降，或者带来不佳的收敛性能。 第二，LMS算法的收敛速度受到输入信号x(n)自相关矩阵Rxx的特征值分布 范围λmax/λmin影响。一般来说，有色信号的特征值发散度大，当输入信号是语音 信号等有色信号时，LMS算法收敛缓慢；当输入信号是诸如白噪声的平稳信号时， LMS算法能够迅速收敛。 基于以上分析，对LMS算法的改进可以从以下两点入手：

1.调整步长因子取值。

由上面分析可知，在稳态性能与收敛速度之间需要折中考虑。在实际应用中， 可以采取可变步长的方法对二者进行折中。在算法运行初期采用较大的步长，加 快收敛速度；然后采用比较小的步长以获得较好的稳态性能。NLMS算法、VMLMS 算法、VSNLMS算法都采取了这样的方法，本论文将在下面进行介绍。

2.降低输入信号的相关性。

由于算法的收敛性与自相关矩阵特征值分布范围有关，而语音等有色信号的 自相关矩阵特征值分布范围较广，故可以从输入信号的角度来提高算法的性能， 采用正交变换的方法等来降低信号的相关性等等。

2.2 NLMS算法

归一化LMS算法(normalized LMS)[ 58]是LMS算法的改进。NLMS算法的稳定 性和收敛速度比LMS算法更佳。在NLMS算法中，为了避免输入信号过大而导致梯度噪声过大，将LMS算法的滤波器权系数更新方程中(n)x(n)e(n)项对输入信 号功率的当前能量Pxx(n)进行归一化处理。 LMS算法的权矢量更新方程为：

△w(n)是权系数更新时的调整量。 后验滤波误差e(n)为：

e(n)是先验滤波误差，取值一般不为0。为了使算法收敛最快，必须让后验滤波误差e(n)为0，即：

由于0<μ（n）<1/λM可得

可得NLMS算法的权矢量更新方程

而实际中上式并不常用，取而代之的是式：

式中μ是正的常量，δ是为了保证数值稳定性在分母中设置的一个很小的大于零的常数。当且仅当步长满足μ∈(0,2)时，算法收敛。

2.3 VSLMS算法

传统的变步长ＬＭＳ算法引入的变步长函数一般是基于ｓｉｇｍｏｉｄ函数的，但是，针对目前的研究现状能同时兼顾快的收敛速度和小的稳态误差的文献不多。为此，本文引入一种基于范数的双曲正切函数的变步长ＬＭＳ算法［６－７］。具体的计算公式如下所示：

式中，ｘ（ｎ）表示ｎ时刻的输入信号；ｗ（ｎ）表示ｎ时刻的加权系数；ｄ（ｎ）表示期望信号；ｅ（ｎ）表示误差信号；ｐ（ｎ）为新引入的非线性函数，且ｐ（０）＝０；μ（ｎ）表示步长函数；β控制步长形状；α控制步长大小；ａ和ｂ 用于调整收敛速度，且０＜ａ＜１，０＜ｂ＜１。从此算法可以看出，引入的变步长函数，不仅与滤波器权向量有关，还与输入信号有关。另外，通过ｐ２ 来调节步长，还可以提高抗干扰能力。因此该算法较现有的其它变步长ＬＭＳ算法在收敛速度、跟踪速度和稳态误差方面均有提高。

2.4 VSNLMS算法

VSNLMS(Variable step-size NLMS，变步长NLMS)[39]也是为了解决LMS算法 和NLMS算法的稳态误差和收敛速度之间的矛盾而提出的。在NLMS算法中，步长因子μ与输入信号功率成反比。VSNLMS算法为步长因子μ(n)设置了上限，在 每个采样周期都要计算μ(n)的最大值。 在VSNLMS算法中，第n个采样点的步长因子表示为μ(n)，每一时刻的滤波 器系数向量w(n)都对应着不同的μ(n)。VSNLMS算法过程如下

gi(n)是随机梯度项，L是自适应滤波器的阶数，ρ是一任意常数。

VSNLMS算法的自适应滤波器的收敛性能可由μ(n)控制，从而能够确保自适应滤波器的收敛速度。步长因子向量μ(n)的范围为

MATLAB代码：

LMS算法

% simLMS.m --- 
%
% Description
%   function [y,e,W] = simLMS(x,d,mu,W0)
​%   function [y,e,W] = simLMS(x,d,mu,W0,v)
%       A simple LMS Adaptive Filter
% Update equation
%       e(i) = d(i) - y(i);
% or    e(i) = d(i) + v(i) - y(i);
%
%       W(next) = W(now) + (mu*e(i)).*u;
% Parameters
%   x:          Input signal
%   d:          Desired output
%   mu:         Stepsize
%   W0:         Initial value of weights
% Return
%   y:          Output signal
%   e:          Error
%   W:          Weights
​
function [y,e,W] = simLMS(x,d,mu,W0,varargin)
error(nargchk(4,5,nargin));
if (nargin == 5)
    v = varargin{1};
else
    v = zeros(1,length(x));
end
W = W0;
L = length(W);
y = zeros(1,length(x));
e = zeros(1,length(x));
u = zeros(1,L); % input_sav vector
for i = 1 : length(x)
    u(2:L) = u(1:L-1);
    u(1) = x(i);
    y(i) = W * u';
    e(i) = d(i) + v(i) - y(i);
    W = W + (mu*e(i)).*u;
end
​

NLMS算法

% simNLMS.m --- 
%
% Description
%   function [y,e,W] = simNLMS(x,d,mu,W0)
%   function [y,e,W] = simNLMS(x,d,mu,W0,v)
%       A simple NLMS Adaptive Filter
% Update equation
%       e(i) = d(i) - y(i);
% or    e(i) = d(i) + v(i) - y(i);
%
%       W(next) = W(now) + (mu*e(i)/(u*u')).*u;
% Parameters
%   x:          Input signal
%   d:          Desired output
%   mu:         Normalized stepsize, 0 < mu < 2
%   W0:         Initial value of weights
% Return
%   y:          Output signal
%   e:          Error
%   W:          Weights
​
function [y,e,W] = simNLMS(x,d,mu,W0,varargin)
% 0 < mu < 2
error(nargchk(4,5,nargin));
if (nargin == 5)
    v = varargin{1};
else
    v = zeros(1,length(x));
end
W = W0;
L = length(W);
y = zeros(1,length(x));
e = zeros(1,length(x));
u = zeros(1,L); % input_sav vector
for i = 1 : length(x)
    u(2:L) = u(1:L-1);
    u(1) = x(i);
    y(i) = W * u';
    e(i) = d(i) +v(i) - y(i);
    W = W + (mu*e(i)/(u*u')).*u;
end
​

VSLMS算法

% simVSLMS.m --- 
%
% Description
%   function [y,e,mu,W] = simVSLMS(handle,x,d,W0)
%   function [y,e,mu,W] = simVSLMS(handle,x,d,W0,v)
%   function [y,e,mu,W] = simVSLMS(handle,x,d,W0,'USERPAR',userpar)
%   function [y,e,mu,W] = simVSLMS(handle,x,d,W0,'USERPAR',userpar,v)
%       A simple VS-LMS Adaptive Filter
% Update equation
%       e(i) = d(i) - y(i);
% or    e(i) = d(i) + v(i) - y(i);
%
%       W(next) = W(now) + (mu*e(i)).*u;
% Parameters
%   handle:     The handle of the step-size update function
%   x:          Input signal
%   d:          Desired output
%   mu:         Stepsize
%   W0:         Initial value of weights
% Return
%   y:          Output signal
%   e:          Error
​%   mu:         stepsize array
%   W:          Weights
​
function [y,e,mu,W] = simVSLMS(handle,x,d,W0,varargin)
error(nargchk(4,7,nargin));
v = zeros(1,length(x));
if (nargin > 4)
    if (nargin == 5)
        v = varargin{1};
    end
    if (varargin{1} == 'USERPAR')
        userpar = varargin{2};
        if (nargin == 7)
            v = varargin{3};
        end
    end
end
W = W0;
L = length(W);
y = zeros(1,length(x));
e = zeros(1,length(x));
mu = zeros(1,length(x));
u = zeros(1,L); % input_sav vector
for i = 1 : length(x)
    u(2:L) = u(1:L-1);
    u(1) = x(i);
    y(i) = W * u';
    e(i) = d(i) + v(i) - y(i);
    if (nargin>5)
        mu(i) = handle(u,W,i,y,e,mu,userpar);
    else
        mu(i) = handle(u,W,i,y,e,mu);
    end
    W = W + (mu(i)*e(i)).*u;
end
​

VSNLMS算法

% simVSNLMS.m --- 
%
​% Description
%   function [y,e,W] = simVSNLMS(handle,x,d,W0,delta)
%   function [y,e,W] = simVSNLMS(handle,x,d,W0,delta,v)
%   function [y,e,mu,W] = simVSNLMS(handle,x,d,W0,delta,'USERPAR',userpar)
%   function [y,e,mu,W] = simVSNLMS(handle,x,d,W0,delta,'USERPAR',userpar,v)
%       A simple VS-NLMS Adaptive Filter
% Update equation
%       e(i) = d(i) - y(i);
% or    e(i) = d(i) + v(i) - y(i);
%
%       W = W + (mu(i)*e(i)/(delta + u*u')).*u;
% Parameters
%   handle:     The handle of the step-size update function
%   x:          Input signal
%   d:          Desired output
%   mu:         Stepsize
%   W0:         Initial value of weights
%   delta:      Adjustion
% Return
%   y:          Output signal
%   e:          Error
%   mu:         stepsize array
%   W:          Weights
​
function [y,e,mu,W] = simVSNLMS(handle,x,d,W0,delta,varargin)
error(nargchk(5,8,nargin));
v = zeros(1,length(x));
if (nargin > 5)
    if (nargin == 6)
        v = varargin{1};
    end
    if (varargin{1} == 'USERPAR')
        userpar = varargin{2};
        if (nargin == 8)
            v = varargin{3};
        end
    end
end
W = W0;
L = length(W);
y = zeros(1,length(x));
e = zeros(1,length(x));
mu = zeros(1,length(x));
u = zeros(1,L); % input_sav vector
for i = 1 : length(x)
    u(2:L) = u(1:L-1);
    u(1) = x(i);
    y(i) = W * u';
    e(i) = d(i) + v(i) - y(i);
    if (nargin>6)
        mu(i) = handle(u,W,i,y,e,mu,userpar);
    else
        mu(i) = handle(u,W,i,y,e,mu);
    end
    W = W + (mu(i)*e(i)/(delta + u*u')).*u;
end
​

function stepsize = stdUpdateFunc(u,W,n,y,e,mu,varargin)
% u:         Filter input memory
% W:         last weights
% n:         current time
% y:         output array
​% e:         error array
% mu:        stepsize array
% varargin:  more parameters user wants to input
userpar = varargin{1};
a = userpar(1);
p = userpar(2);
stepsize = a*abs(e(n))^p/(u*u');

%% 需自行加入噪声n
s = randn(1,5000);
d1 = filter([0.1,0.2,0.3,0.4,0.4,0.3,0.2,0.1],1,s(1:2500));
d2 = filter(1,[1,-0.2,0.1],s(2501:5000));
d(1:2500) = d1;
d(2501:5000) = d2;
x = s + n;
w0 = zeros(1,12);
h1 = @testFunc; %设置函数句柄
[y1,e1,w1] = simNLMS(x,d,0.2,w0);
[y2,e2,mu2,w2] = simVSNLMS(h1,x,d,w0,0,'USERPAR',[0.2,0.5]);
[y3,e3,mu3,w3] = simVSNLMS(h1,x,d,w0,0,'USERPAR',[0.5,0.5]);