自适应滤波LMS与RLS的matlab实现

1、MATLAB仿真实现LMS和RLS算法的二阶AR模型及仿真结果分析一、题目概述：二阶AR模型如图1a所示，可以如下差分方程表示： （1）图1a其中，v(n)是均值为0、方差为0.965的高斯白噪声序列。a1，a2为描述性参数，设x(-1)=x(-2)=0，权值w10=w20=0，=0.04推导最优滤波权值（理论分析一下）。按此参数设置，由计算机仿真模拟权值收敛曲线并画出，改变步长在此模拟权值变化规律。对仿真结果进行说明。应用RLS算法再次模拟最优滤波权值。解答思路：（1）高斯白噪声用normrnd函数产生均值为0、方差为0.965的正态分布随机1*N矩阵来实现。随后的产生的信号用题目中的二阶A

2、R模型根据公式（1）产生，激励源是之前产生的高斯白噪声。（2）信号长度N取为2000点，用以观察滤波器权值变化从而估计滤波器系数，得到其收敛值。（3）仿真时分别仿真了单次LMS算法和RLS算法下的收敛性能以及100次取平均后的LMS和RLS算法的收敛性能，以便更好的比较观察二者的特性。（4）在用不同的分别取3个不同的值仿真LMS算法时，值分别取为0.001,0.003,0.006；用3个不同的值仿真RLS算法时值分别取为1,0.98,0.94，从而分析不同步长因子、不同遗忘因子对相应算法收敛效果的影响。二、 算法简介1自适应算法的基本原理自适应算法的基本信号关系如下图所示：图 1b 自适应滤波

3、器框图输入信号x(n)通过参数可调的数字滤波器后产生输出信号y(n)，将其与参考信号d(n)进行比较，形成误差信号e(n)。e(n)通过某种自适应算法对滤波器参数进行调整，最终是e(n)的均方值最小。当误差信号e(n)的均方误差达到最小的时候，可以证明信号y(n)是信号d(n)的最佳估计。2. LMS算法简介LMS算法采用平方误差最小的原则代替最小均方误差最小的原则，信号基本关系如下： （2）写成矩阵型式为： （3）式（3）中，W(n) 为n 时刻自适应滤波器的权值，N为自适应滤波器的阶数，本设计中取为2000；X( n) 为n 时刻自适应滤波器的参考输入矢量，由最近N 个信号采样值构成，；d

4、 ( n) 是期望的输出值；e ( n) 为自适应滤波器的输出误差调节信号(简称失调信号) ；是控制自适应速度与稳定性的增益常数，又叫收敛因子或步长因子。3 RLS算法简介RLS算法是用二乘方的时间平均的最小化准则取代最小均方准则，并按时间进行迭代计算。其基本原理如下所示：遗忘因子，它是小于等于1的正数。参考信号，也可称为期望信号。第n次迭代的权值。均方误差。RLS算法的准则为： （4）上式越旧的数据对的影响越小。通过计算推导得到系数的迭代方程为： wn=wn-1+kne*(n) （5）式(5)中，增量k(n)和误差e*(n)计算公式如下： (6) e*n=d*n-xTn*w(n-1) （7）

5、 式（6）中T(n)= R-1n，也就是当前时刻自相关矩阵的逆。按如下方式更新： = (-* )/ （8）由上边分析可知，RLS算法递推的步骤如下：1. 在时刻n，和也已经存储在滤波器的相应器件中2. 利用公式（5）、（6）、（7）和（8）计算T(n)、w(n)、k(n)、e*n，并得到滤波器的输出相应和误差即： （9） （10）3. 进入第次迭代这样做的优点是收敛速度快，而且适用于非平稳信号的自适应处理缺点是每次迭代时都要知道输入信号和参考信号，计算量比较大三、 仿真过程仿真过程按照如下过程进行（1） 信号产生：首先产生高斯白噪声序列w(n)，然后将此通过一个参数为a1=-0.195，a2=

6、0.95简单的二阶自回归滤波器生成信号。（2）将step（1）生成的信号通过LMS和RLS自适应滤波器进行处理（3）通过改变值对w1,w2收敛速度的影响来分析LMS算法的性能以及通过改变值对w1,w2收敛速度的影响来分析RLS算法的性能。（4）绘制相应图形曲线四、 仿真以及结果分析 信号和高斯白噪声波形如图2所示：图2 信号和高斯白噪声波形 图2中，上边的图形为信号波形，下边的为加性高斯白噪声。图3（a）LMS算法下单次收敛曲线图3(b) LMS算法下百次平均收敛曲线分析1：图3中，a、b两幅图分别为单次实现的LMS算法下最优权值变化过程和100次仿真实现后取平均值做的图，两个权值初始值由已知

7、条件设置为0，之后收敛到两个定值。 a图展现了滤波器权系数迭代更新的过程，可以发现其并不是平滑的变化，而是随机起伏的，跟最陡下降法不一样，这充分说明了其权向量是一个随机过程向量,梯度噪声对其产生了一定的影响。b图给出了100次独立实验得到的平均权向量的估计，即() ，其中是第t次独立实验中第n次迭代得到的权向量，T是独立实验次数。可以发现，多次独立实验得到的平均权向量的估计平滑了随机梯度引入的梯度噪声，十分接近理论最优权向量，曲线足够平滑，噪声影响很小。图3中a、b两图皆有较好的收敛特性，即使1次实现也能很好的逼近最优权值。图4（a）RLS算法下单次收敛曲线图4（b）RLS算法下百次平均收敛曲

8、线 分析2：图4中，a、b两幅图分别为单次实现的RLS算法下最优权值变化过程和100次仿真实现后取平均值做的图，a图中可看出权值已有较好收敛特性，两个权值初始值由已知条件设置为0，之后收敛到两个定值：-a1=0.195和-a2=-0.95，但是曲线不够平滑、噪声较大，b图经过取平均后噪声影响已经很小。单次实现和多次实现的联系与区别也与LMS算法基本相同。但是可以清晰地看出RLS算法收敛极。将两种算法整合于一个图中可以看出鲜明的对比。图5 LMS算法和RLS算法收敛曲线对比分析3：从图5可以看出，在和一定的情况下，RLS比LMS具有更快的收敛速度。在算法的前期收敛段，RLS算法的收敛速度要明显高

9、于LMS算法。 总的来说，LMS算法的收敛性和步长有关，受协方差矩阵的特征根影响，然而RLS算法一定是收敛的，且收敛速度很快，不过其迭代过程由于要求增益因子和逆矩阵等，使得其算法复杂度高，计算量比较大。接下来修改步长值，观察曲线的收敛情况。图6给出了=0.01, =0.03和=0.06三种情况下的最优权值变化曲线，由此分析不同步长对曲线收敛性产生的影响。图6 对LMS收敛速度的影响分析4： 可以看出LMS算法中的步长参数决定抽头权向量在每步迭代中的更新量，是影响算法收敛速度的关键参数，其决定了LMS算法学习过程的快慢。 图6展示了当收敛步长值变小时，其均方误差的收敛速度也相应减慢，降低对实变系

10、统的跟踪速度，不能及时调整至最优权值。 由前面的理论推导我们可以知道，在这个范围内，越大，均方误差收敛速度越快，收敛速度和步长因子成正比。可是如果大于这个范围会造成不稳定，较大步长会造成较大的稳态误差，带来算法的发散。图7 对RLS收敛速度的影响 分析5：RLS算法中遗忘因子的作用是对离n时刻越近的误差加比较大的权重，而对离n时刻越远的误差加比较小的比重。遗忘因子的选择对RLS算法的性能起决定性的作用。是个小于等于1的数，如果越小，能量信号就越接近最新的误差平方，对前面的误差遗忘的越快，跟踪效果就越好。但是，递推RLS算法中的误差是由期望信号决定的，如果很小则误差信号对期望信号的依赖性就会很大

11、，所以，输出信号就很接近期望信号。从图7可以看出，遗忘因子越小，系统的跟踪能力越强，收敛的越快（即变平稳得越快），但是收敛值比较大；越大，系统跟踪能力减弱，收敛时均方误差也越小，学习曲线收敛值越小。所以，我们在用普通递推RLS算法时，一定要对有个准确的取值，才能保证系统性能的最佳状态。 总结：LMS算法其优点是结构简单，算法复杂度低，易于实现，稳定性高，便于硬件实现，但是这种算法收敛速度慢，对快速变化的信号不适合。RLS算法是基于最小二乘准则的精确方法，它的收敛速度快，稳定性强， 因此常被应用于实时系统识别和快速启动的信道均衡。附件（程序代码）：%(1)信号序列与高斯白噪声的产生%参数初始化a

12、1=-0.195; %生成信号所用AR（2）滤波器的参数a2=0.95;N=2000; %信号长度 %信号及白噪声信号序列的初始化x=zeros(1,N)' %信号的初始化sigmasqu=0.0965;%噪声方差sigma=sqrt(sigmasqu);%噪声标准差w=normrnd(0,sigma,N)' %高斯白噪声的初始化，均值为0，方差为0.0965 x(1)=w(1); %信号前两点的初始赋值x(2)=-a1*x(1)+w(2); %信号序列的产生for i=3:N x(i)=-a1*x(i-1)-a2*x(i-2)+w(i);%信号由AR（2）产生end %绘制信

13、号和高斯白噪声波形figure(1)subplot(2,1,1),plot(1:N,x,'b-')axis(0,2000,-5,5)title('基于AR（2）模型产生的信号x');xlabel('信号长度n');ylabel('x(n)');subplot(2,1,2),plot(1:N,w,'r-');axis(0,2000,-5,5)title('基于AR（2）模型产生的高斯白噪声w(n)');xlabel('信号长度n');ylabel('w(n)'); %-

14、LMS100次叠加% N=2000;M=100;%计算的次数w1=zeros(N,M);w2=zeros(N,M);I1=eye(2);E1=zeros(N,M);wa1=zeros(N,1);wa2=zeros(N,1);en=zeros(N,1); for k=1:M e=zeros(1,N)'%定义误差向量 %根据RLS算法进行递推 x=zeros(1,N)' %信号的初始化 sigmasqu=0.0965;%噪声方差 sigma=sqrt(sigmasqu);%噪声标准差 w=normrnd(0,sigma,N)' %高斯白噪声的初始化，均值为0，方差为0.09

15、65 x(1)=w(1); %信号前两点的初始赋值 x(2)=-a1*x(1)+w(2); %信号序列的产生 for i=3:N x(i)=-a1*x(i-1)-a2*x(i-2)+w(i);%信号由AR（2）产生 end L=2; %滤波器长度 u=0.04; %LMS算法下自适应增益常数初始化 wL=zeros(L,N);%LMS滤波器的权值初始化 for i=(L+1):N X=x(i-1:-1:(i-L); y(i)=X'*wL(:,i); %i时刻输出信号 e(i)=x(i)-y(i); %i时刻误差信号 wL(:,(i+1)=wL(:,i)+2*u*e(i)*X; %i时刻

16、滤波器的权值 end;% a1R=wR(1,1:n); % a1在RLS算法下值的变化% a2R=wR(2,1:n); % a2在RLS算法下值的变化 w1(:,k)=wL(1,1:N)' w2(:,k)=wL(2,1:N)'%将每次计算得到的权矢量值储存 E1(:,k)=e(:,1);%将每次计算得到的误差储存 %求权矢量和误差的M次的平均值 wa1(:,1)=wa1(:,1)+w1(:,k); wa2(:,1)=wa2(:,1)+w2(:,k); en(:,1)=en(:,1)+E1(:,k); end %-止LMS100次叠加% %绘制LMS算法下a1、a2单次及100次

17、平均收敛曲线figure(2)plot(1:N,w1(1:N,M)','r-',1:N,w2(1:N,M)','b-',1:N,-a2,'k-',1:N,-a1,'k-');legend('LMS-W1变化','LMS-W2变化','W1收敛值','W2收敛值',0); % 图例title('LMS单次实现的权值收敛曲线');xlabel('信号长度n');ylabel('权值W'); figure(3)

18、plot(1:N,wa1(1:N,1)./M,'r-',1:N,wa2(1:N,1)./M,'b-',1:N,-a2,'k-',1:N,-a1,'k-');%作出100次计算权矢量的平均变化曲线xlabel('N');ylabel('W(N)');title('LMS算法下100次仿真W1(n)和W2(n)的平均变化曲线') % %RLS滤波 %-RLS100次叠加% N=2000;M=100;%计算的次数w3=zeros(N,M);w4=zeros(N,M);I2=eye(2);E

19、2=zeros(N,M);wa3=zeros(N,1);wa4=zeros(N,1);en2=zeros(N,1); for k=1:M e=zeros(1,N)'%定义误差向量 %根据RLS算法进行递推 x=zeros(1,N)' %信号的初始化 sigmasqu=0.0965;%噪声方差 sigma=sqrt(sigmasqu);%噪声标准差 w=normrnd(0,sigma,N)' %高斯白噪声的初始化，均值为0，方差为0.0965 x(1)=w(1); %信号前两点的初始赋值 x(2)=-a1*x(1)+w(2); %信号序列的产生 for i=3:N x(i

20、)=-a1*x(i-1)-a2*x(i-2)+w(i);%信号由AR（2）产生 end L=2; %滤波器长度 lam=0.98; %RLS算法下lambda取值 wR=zeros(L,N);%权系数，初值为0 T=eye(L,L)*10;% %RLS算法下T参数的初始化,T初始值为10 for i=(L+1):N X=x(i-1:-1:(i-L); K=(T*X)/(lam+X'*T*X);%i时刻增益值 e1=x(i)-wR(:,i-1)'*X; wR(:,i)=wR(:,i-1)+K*e1; %i时刻权值 y(i)=wR(:,i)'*X;%输出信号 e(i)=x(

21、i)-y(i);%预测误差 T=(T-K*X'*T)/lam; %i时刻的维纳解 end; w3(:,k)=wR(1,:)' w4(:,k)=wR(2,:)'%将每次计算得到的权矢量值储存 E2(:,k)=e(:,1);%将每次计算得到的误差储存 %求权矢量和误差的M次的平均值 wa3(:,1)=wa3(:,1)+w3(:,k); wa4(:,1)=wa4(:,1)+w4(:,k); en2(:,1)=en2(:,1)+E2(:,k); end %-止RLS100次叠加% figure(4)plot(1:N,w3(1:N,M)','r-',1:N

22、,w4(1:N,M)','b-',1:N,-a2,'k-',1:N,-a1,'k-');legend('RLS-W1变化','RLS-W2变化','W1收敛值','W2收敛值',0); % 图例title('RLS单次实现权值收敛曲线');xlabel('信号长度n');ylabel('权值W'); figure(5)plot(1:N,wa3(1:N,1)./M,'r-',1:N,wa4(1:N,1)./M,&#

23、39;b-',1:N,-a2,'k-',1:N,-a1,'k-');%作出100次计算权矢量的平均变化曲线xlabel('N');ylabel('W(N)');title('RLS算法下100次仿真W1(n)和W2(n)的平均变化曲线线') %绘制LMS与RLS算法下a1、a2收敛曲线figure(6)plot(1:N,w1(1:N,1)','r-.',1:N,w2(1:N,1)','b:',1:N,w3(1:N,1)','g-',1:

24、N,w4(1:N,1)','m-',1:N,-a2,'k-',1:N,-a1,'k-');legend('LMS-W1变化','LMS-W2变化','RLS-W1变化','RLS-W2变化','W1收敛值','W2收敛值',0); % 图例title('LMS与RLS算法单次实现的对比');xlabel('信号长度n');ylabel('权值W'); %(3)LMS算法下不同u值的参数收敛曲线wL

25、=zeros(L,N,3);eL=zeros(N,3); % LMS算法下误差初始化yL=zeros(N,3); % LMS算法下滤波器输出初始化u=0.01,0.03,0.06; %不同的u值 for j=1:3 for i=(L+1):N yL(i,j)=x(i-1:-1:i-2)'*wL(1:L,i-1,j); eL(i,j)=x(i)-yL(i,j); wL(1:L,i,j)=wL(1:L,i-1,j)+2*u(j)*eL(i,j)*x(i-1:-1:i-L); endenda1L1=wL(1,1:N,1);a1L2=wL(1,1:N,2);a1L3=wL(1,1:N,3);a

26、2L1=wL(2,1:N,1);a2L2=wL(2,1:N,2);a2L3=wL(2,1:N,3); figure(7)plot(1:N,a1L1,'b-',1:N,a1L2,'r:',1:N,a1L3,'g-.',1:N,-a1,'k',1:N,a2L1,'b-',1:N,a2L2,'r:',1:N,a2L3,'g-.',1:N,-a2,'k') % 画图显示不同u值下LMS算法性能差别legend('u=0.01','u=0.03','u=0.06','W收敛值',0) % 图例title('LMS算法下不同的u值对W