[Matlab]LMS滤波器设计

[Matlab]自适应（LMS）滤波器设计

自适应滤波是近年以来发展起来的一种最佳滤波方法。它是在维纳滤波,Kalman滤波等线性滤波基础上发展起来的一种最佳滤波方法。由于它具有更强的适应性和更优的滤波性能。从而在工程实际中,尤其在信息处理技术中得到了广泛的应用。自适应滤波存在于信号处理、控制、图像处理等许多不同领域，它是一种智能更有针对性的滤波方法，通常用于去噪。

背景介绍

自适应滤波的研究对象是具有不确定的系统或信息过程。这里的“不确定性”是指所研究的处理信息过程及其环境的数学模型不是完全确定的。其中包含一些未知因素和随机因素。

任何一个实际的信息过程都具有不同程度的不确定性,这些不确定性有时表现在过程内部,有时表现在过程外部。从过程内部来讲,描述研究对象即信息动态过程的数学模型的结构和参数是设计者事先并不一定能确切知道的。作为外部环境对信息过程的影响,可以等效地用扰动来表示。这些扰动通常是不可测的,它们可能是确定性的,也可能是随机的。此外,还有一些测量噪音 [1] 也以不同的途径影响信息过程。这些扰动和噪声的统计特性常常是未知的。面对这些客观存在的各式各样的不确定性,如何综合处理该信息过程,并使得某一些指定的性能指标达到最优或近似最优,这就是自适应滤波所要解决的问题。

LMS与维纳滤波器（Wiener Filter）的区别

这里介绍的LMS/NLMS，通常逐点处理，对应思路是：随机梯度下降；
对于Wiener Filter，给定准则函数J，随机/批量梯度都可以得出最优解；
LMS虽然基于梯度下降，但准则仅仅是统计意义且通常引入误差，可以定义为J0J_0J0，简而言之J通常不等于J0J_0J0，得出的最优解wow_{o}wo,自然也通常不等于维纳最优解；
分析LMS通常会分析稳定性，稳定性是基于Wiener解，之前已给出分析。但LMS是Wiener解的近似，所以：迭代步长的稳定性，严格适用于Wiener解，对于LMS只是一种近似参考，并没有充分的理论依据。

LMS与维纳滤波器原理

自适应滤波器的原理如图1所示。

图中x（k）表示k时刻的输入信号值，y（k）表示 k时刻的输出信号值，d（k）表示k 的参考信号值或所期望响应信号值，误差信号e（k）为d（k）与y（k）之差。自适应数字滤波器的滤波参数受误差信号e（k）的控制，根据e（k）的值而自动调整，使之适合下一时刻的输入x（k+1），以便使输出y（k+1）接近于所期望的参考信号d(k+1)。

自适应滤波器可以分为线性自适应滤波器和非线性自适应滤波器。非线性自适应滤波器包括Voetlrra滤波器和基于神经网络的自适应滤波器。非线性自适应滤波器具有更强的信号处理 [3] 能力。但是，由于非线性自适应滤波器的计算较复杂，实际用得最多的仍然是线性自适应滤波器。

典型算法

对自适应滤波算法 [4] 的研究是当今自适应信号处理中最为活跃的研究课题之一。自适应滤波算法广泛应用于系统辨识、回波消除、自适应谱线增强、自适应信道均衡、语音线性预测、自适应天线阵等诸多领域中。总之,寻求收敛速度快,计算复杂性低,数值稳定性好的自适应滤波算法是研究人员不断努力追求的目标。虽然线性自适应滤波器和相应的算法具有结构简单、计算复杂性低的优点而广泛应用于实际,但由于对信号的处理能力有限而在应用中受到限制。由于非线性自适应滤波器,如Voletrra滤波器和基于神经网络的自适应滤波器,具有更强的信号处理能力,已成为自适应信号处理中的一个研究热点。其中较典型的几种算法包括:

LMS自适应滤波算法
RLS自适应滤波算法
变换域自适应滤波算法
仿射投影算法
共扼梯度算法
基于子带分解的自适应滤波算法
基于QR分解的自适应滤波算法

算法性能评价

变步长的自适应滤波算法 [6] 虽然解决了收敛速度、时变系统跟踪速度与收敛精度方面对算法调整步长因子u的矛盾,但变步长中的其它参数的选取还需实验来确定,应用起来不太方便。对RLS算法的各种改进,其目的均是保留RLS算法收敛速度快的特点而降低其计算复杂性。变换域类算法亦是想通过作某些正交变换使输入信号自相关矩阵的特征值发散程度变小,提高收敛速度。而仿射投影算法的性能介于LMS算法和RLS算法之间。共扼梯度自适应滤波算法的提出是为了降低RLS类算法的复杂性和克服某些快速RLS算法存在的数值稳定性问题。信号的子带分解能降低输入信号的自相关矩阵的特征值发散程度,从而加快自适应滤波算法的收敛速度,同时便于并行处理,带来了一定的灵活性。矩阵的QR分解具有良好的数值稳定性。

自适应（LMS）滤波器设计(1):

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输入参数:
%     xn   输入的信号序列      (列向量)
%     dn   所期望的响应序列    (列向量)
%     M    滤波器的阶数        (标量)
%     mu   收敛因子(步长)      (标量) 要求大于0,小于xn的相关矩阵最大特征值的倒数
% 输出参数:
%     W    滤波器的权值矩阵     (矩阵)
%          大小为M x itr,
%     en   误差序列(itr x 1)    (列向量)
%     yn   实际输出序列         (列向量)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [yn,W,en]=LMS1(xn,dn,M,mu)
itr = length(xn);
en = zeros(itr,1);     % 误差序列,en(k)表示第k次迭代时预期输出与实际输入的误差
W  = zeros(M,itr);     % 每一行代表一个加权参量,每一列代表-次迭代,初始为0
% 迭代计算
for k = M:itr              % 第k次迭代x = xn(k:-1:k-M+1);    % 滤波器M个抽头的输入y = W(:,k-1).' * x;    % 滤波器的输出en(k) = dn(k) - y ;    % 第k次迭代的误差W(:,k) = W(:,k-1) + 2*mu*en(k)*x; % 滤波器权值计算的迭代式
end
% 求最优时滤波器的输出序列  r如果没有yn返回参数可以不要下面的
yn = inf * ones(size(xn)); % inf 是无穷大的意思
for k = M:length(xn)x = xn(k:-1:k-M+1);yn(k) = W(:,end).'* x;%用最后得到的最佳估计得到输出
end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%案例分析
clc;clear all; close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 周期信号的产生
t=0:0.5:100;
xs=30*sin(t);
figure;
subplot(2,1,1);
plot(t,xs);grid;
ylabel('幅值/v');
xlabel('时间/t');
title('{输入原始性信号}');
% 噪声信号的产生
randn('state',sum(200*clock));
xn=randn(1,201);
subplot(2,1,2);
plot(t,xn);grid;
ylabel('幅值/v');
xlabel('时间/t');
title('{随机噪声信号}');% 信号滤波
xn = xs+xn;
xn = xn.' ;   % 输入信号序列
dn = xs.' ;   % 预期结果序列
M  = 20 ;   % 滤波器的阶数rho_max = max(eig(xn*xn.'));   % 输入信号相关矩阵的最大特征值
mu = (1/rho_max) ;    % 收敛因子 0 < mu < 1/rho
[yn,W,en] = LMS1(xn,dn,M,mu);% 绘制滤波器输入信号
figure;
subplot(2,1,1);
plot(t,xn);grid;
ylabel('幅值');
xlabel('时间');
title('{滤波器输入信号}');% 绘制自适应滤波器输出信号
subplot(2,1,2);
plot(t,yn);grid;
ylabel('幅值');
xlabel('时间');
title('{自适应滤波器输出信号}');% 绘制自适应滤波器输出信号,预期输出信号和两者的误差
figure
plot(t,yn,'b.-',t,dn,'g--',t,dn-yn,'r');grid;
legend('自适应滤波器输出','预期输出','误差');
ylabel('幅值');
xlabel('时间');
title('{自适应滤波器}');

自适应（LMS）滤波器设计(2)–DSP-2 最陡下降法:

%LMS算法
% 输入参数:
%     Y   输入的信号序列      (列向量)
%     X   所期望的响应序列    (列向量)
%     COUNT 输入信号长度数     (标量)
%     mu   收敛因子(步长)      (标量) 要求大于0,小于xn的相关矩阵最大特征值的倒数
% 输出参数:
%     W    滤波器的权值矩阵     (矩阵)
%          大小为M x itr,
%     en   误差序列(itr x 1)    (列向量)
%     yn   实际输出序列         (列向量)
function [e,w1,w2]=LMS(COUNT,X,Y,u)
e=zeros(1,COUNT);
w1=zeros(1,COUNT);
w2=zeros(1,COUNT);
w1(1)=3;
w2(1)=-4;
for i=1:COUNTif(i-1==0)yy=w1(i)*X(i)+w2(i)*X(16); elseyy=w1(i)*X(i)+w2(i)*X(i-1);end;e(i)=Y(i)-yy;if(i<COUNT)w1(i+1)=w1(i)+u*e(i)*X(i);end;if(0==i-1)w2(i+1)=w2(i)+u*e(i)*X(16);elseif(i<COUNT)w2(i+1)=w2(i)+u*e(i)*X(i-1);end;end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Expectation函数实现
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function EX=Expectation(x,y,N)
EX=0;
for n=1:1:NEX=EX+x(n)*y(n)/N;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%最陡下降法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [w1,w2]=Steepest_Algorithm(R,P,u,COUNT)
VGn=zeros(2,1);
H=zeros(2,1);
w1=zeros(1,COUNT);
w2=zeros(1,COUNT);
w1(1)=3;
w2(1)=-4;
for i=1:COUNTH=[w1(i),w2(i)]';VGn=2*R*H-2*P';w1(i+1)=w1(i)-0.5*u*VGn(1,1);w2(i+1)=w2(i)-0.5*u*VGn(2,1);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear all; close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
COUNT=300;
n=1:COUNT;
N=sin(2*n*pi/16+pi/10);
X=2^0.5*sin(2*n*pi/16);
X1=2^0.5*sin(2*(n-1)*pi/16);
S=sqrt(0.05)*randn(1,COUNT);
%随机信号
Y=S+X;
figure
subplot(311)
plot(n,X);
title('信号图');
subplot(312)
plot(n,S);
title('噪声图');
subplot(313)
plot(n,Y);
title('信号加噪声图');
%误差性能曲面及等值线
ESS=0.05;
ENN=Expectation(N,N,16);
EXX=Expectation(X,X,16);
R(1,1)=EXX;
R(2,2)=EXX;
EXX1=Expectation(X,X1,16);
R(1,2)=EXX1;
R(2,1)=EXX1;
EYY=ENN+ESS;
EYX=Expectation(N,X,16);
EYX1=Expectation(N,X1,16);
P=zeros(1,2);
P(1)=EYX;
P(2)=EYX1;
x = -4:0.05:6;
y = -5:0.05:5;
[h0,h1] = meshgrid(x,y);
z=EYY+R(1,1)*h0.*h0+2*R(1,2)*h0.*h1+R(2,2)*h1.*h1-2*P(1)*h0-2*P(2)*h1;
figure
subplot(1,2,1);
mesh(h0,h1,z);
xlabel('h0');
ylabel('h1');
title('误差性能曲面图');
subplot(1,2,2);
V=0.2:0.2:3;
contour(h0,h1,z,V);
xlabel('h0');
ylabel('h1');
title('等值线');
hold on;
%LMS算法
u=0.4;
[e,w1,w2]=LMS(COUNT,X,Y,u);
length(w1);
length(w2);
plot(w1,w2);
hold on;
% DSP-2 最陡下降法
[w1,w2]=Steepest_Algorithm(R,P,u,COUNT);
plot(w1,w2);
%LMS算法中一次和多次实验中梯度估计和平均值随时间n的变化情况
Jn=zeros(1,COUNT);
Jn=e.^2;
figure
subplot(1,2,1);
plot(n,Jn);
title('LMS算法中1次实验梯度估计');
%200次实验
e_avr=zeros(1,COUNT);
w1_avr=zeros(1,COUNT);
w2_avr=zeros(1,COUNT);
for i=1:200S=sqrt(0.05)*randn(1,COUNT);Y=S+N;[e,w1,w2]=LMS(COUNT,X,Y,u);e_avr=e_avr+e./100;w1_avr=w1_avr+w1./100;w2_avr=w2_avr+w2./100;
end;
subplot(1,2,2);
Jn=e_avr.^2;
plot(n,Jn);
title('LMS算法中200次实验梯度估计曲线图');
figure
plot(w1_avr,w2_avr);
xlabel('h0');
ylabel('h1');
title('LMS算法200次实验中H(n)的平均轨迹曲线图');