基于LMS和RLS的自适应滤波器的应用仿真

湖南大学计算机与通信学院课程作业2题目：基于LMS和RLS的自适应滤波器的应用仿真基于LMS和RLS的自适应滤波器应用仿真1.自适应滤波原理自适应滤波器是指利用前一时刻的结果，自动调节当前时刻的滤波器参数，以适应信号和噪声未知或随机变化的特性，得到有效的输出，主要由参数可调的数字滤波器和自适应算法两部分组成，如图1所示图1自适应滤波器原理图x(n)称为输入信号，y(n)称为输出信号，d(n)称为期望信号或者训练信号，e(n)为误差僖号，其中，e(n)=d(n)-y(n).自适应滤波器的系数(权值)根据误差信号e(n)，通过一定的自适应算法不断的进行改变，以达到使输出信号y(n)最接近期望信号图中参数可调的数字滤波器和自适应算法组成自适应滤波器。自适应滤波算法是滤波器系数权值更新的控制算法，根据输入信号与期望信号以及它们之间的误差信号，自适应滤波算法依据算法准则对滤波器的系数权值进行更新，使其能够使滤波器的输出趋向于期望信号。原理记数字滤波器脉冲响应为：h(k)=[h0(k)h1(k)…hn-1(k)]T输入采样信号为：x(k)=[x(k)x(k-1)…x(k-n-1)]误差信号为：)()()(^kykyke()()()()Tekykhkxk优化过程就是最小化性能指标J(k),它是误差的平方和：21()[()()()]kTiJkyihkxi求使J(k)最小的系数向量h(k),即使J(k)对h(k)的导数为零，也就是0)()(kdhkdJ。把J(k)的表达式代入，得：12[()()()]()0kTiyihkxixi和11()()()()()kkTTTiixiyihkxixi由此得出滤波器系数的最优向量：11()()()()()kTTikTixiyihkxixi这个表达式由输入信号自相关矩阵()xxcx和输入信号与参考信号的相关矩阵()yxck组成，如下所示，维数都为（n,n）:1()()()kTxxickxixi1()()()kTyxickxiyi系数最优向量也可以写成如下形式：1()()()Toptyxxxhkckck自相关和互相关矩阵的递归表达式如下：()(1)()()Txxxxckckxkxk()(1)()()Tyxyxckckykxk把()yxck的递归表达式代入系数向量表达式，得：1()()()Tyxxxhkckck即1()[(1)()()]()TTyxxxhkckxkykck考虑到(1)(1)(1)Tyxxxckhkck可以记1()()[(1)(1)()()]xxxxhkcxckhkykxk用前面得到的表达式求出(1)xxck，并代入上式：1()(){[()()()](1)()()}Txxxxhkcxckxkxkhkykxk或1()(1)()[()()()()(1)]Txxhkhkcxykxkxkxkhk则滤波器系数的递归关系式可以记作1()(1)()[()()()()(1)]Txxhkhkcxykxkxkxkhk其中()()()(1)Tekykxkhke(k)表示先验误差。只因为它是由前一个采样时刻的系数算出的，在实际中，很多时候由于h(k)计算的复杂度而不能应用于实时控制。用δ,I代换()xxck，其中：δ为自适应梯度，I为辨识矩阵（n，n）这时()(1)()()xhkhkkek这时就是一个最小均方准则问题。2.LMS自适应滤波器举例自回归过程的自适应预估器自回归过程是用来描述伴随一些可能性规律出现的统计现象的瞬时估计的随机过程。一阶自回归模型的公式如下：1()(1)()ykaykbk1a是模型的唯一参数，b(k)是零均值白噪声。用一个自适应滤波器生成一个可以对参数1a进行一步预测的一阶自适应预估器。LMS算法可由如下方程表示：)1(^)()(1kykykea)()1()(^)1(^11kekykkaa取N个点估计参数1a,为获取平均值重复M次。而且分别对δ=0.01，δ=0.05，δ=0.1进行计算。参数1a固定在-0.6。程序清单如下：N=500;M=20;n=1;a1=-0.8;h=zeros(M,n+1,3);e=zeros(M,n,3);ford=1:3ifd==1delta=0.01;elsedelta=0.05*(d-1);end;fork=1:Mb=0.2*randn(1,N);y(1)=1;fori=2:Ny(i)=-a1*y(i-1)+b(i);endfori=n+1:Ne(k,i,d)=y(i)-h(k,i,d)*y(i-1);h(k,i+1,d)=h(k,i,d)+delta*y(i-1)*e(k,i,d);endendendford=1:3fori=1:Nem(i,d)=0;hm(i,d)=0;forj=1:Mem(i,d)=em(i,d)+e(j,i,d)^2;hm(i,d)=hm(i,d)+h(j,i,d);endendendfigure(1)semilogy(1:150,em(1:150,1)),holdonsemilogy(1:150,em(1:150,2),'r'),holdonsemilogy(1:150,em(1:150,3),'g'),holdoffaxis([01500.011]),gridtitle('Meansquareerror')xlabel('Samples')gtext('\leftarrowd=0.01');gtext('\leftarrowd=0.05');gtext('\leftarrowd=0.1');figure(2),plot(1:N,hm(1:N,1)),holdonplot(1:N,hm(1:N,2),'r'),holdonplot(1:N,hm(1:N,3),'g'),holdoff,gridtitle('Filtercoeffcientevalution')xlabel('Samples'),gtext('d=0.01'),gtext('d=0.05'),gtext('d=0.1')得到的如下结果图：图2.平均方差误差图3.滤波器系数曲线系数以时间常数的指数曲线收敛，δ越大，时间常数越小3.RLS自适应滤波器的应用仿真从噪声中提取信号输入信号为：000()cos()()xkAwtbk其中b(k)是附加的白噪声。000()cos()rxkAwt应用于RLS自适应滤波器的算法可描述如下：11(1)()()1()(1)()xxTxxckxkgkxkckxk自适应增益行向量，大小（1，n）;()()(1)()Tekykhkxk先验误差()(1)()()hkhkgkek自适应滤波器系数行向量，大小（1,n）111()(1)()()(1)Txxxxxxckckgkxkck输入信号x(k)的自相关转制矩阵，大小（1,n）()()()Tykhkxk自适应滤波输出所研究的滤波器阶数为200，采样周期等于1ms.程序清单如下：N=1000;n=200;k=12;Ts=1e-1b=0.8*randn(1,N);fori=1:Nxr(1,i)=sin(k*2*pi*i/N);x(1,i)=xr(1,i)+b(i);endCxx=10000*eye(n);g=zeros(N,n);h=zeros(N,n);e=zeros(1,N);y=zeros(1,N);tr=zeros(1,N);fori=n+1:Ng(i,:)=(Cxx*x(i-n+1:i)'./(1+x(i-n+1:i)*Cxx*x(i-n+1:i)'))';e(1,i)=xr(i)-h(i-1,:)*x(i-n+1:i)';h(i,:)=h(i-1,:)+e(1,i)*g(i,:);Cxx=Cxx-g(i,:)'*x(i-n+1:i)*Cxx;y(1,i)=h(i,:)*x(i-n+1:i)';tr(1,i)=trace(Cxx);endfigure(1)plot(0:N-n,x(1,n:N)),gridtitle('x(k)inputsingnalinV')xlabel('Samples')figure(2)plot(0:N-n,xr(1,n:N),'r'),gridaxis([0800-1.21.2])title('xr(k)referencesingnalinV')xlabel('Samples')figure(3)plot(0:N-n,e(1,n:N)),holdonplot(0:N-n,y(1,n:N),'r'),holdongridtitle('e(k)errorandy(k)outputinV')xlabel('Samples')gtext('e(k)'),gtext('y(k)')figure(4)plot(0:N-n,h(n:N,1)),holdonplot(0:N-n,h(n:N,2),'r'),holdoffgridtitle('a(n-1)anda(n-2)coeffcientsevolution')xlabel('Samples')figure(5)num1=fliplr(h(N,:));sys1=tf(num1,1,Ts);bode(sys1),holdofftitle('Synthesizedfilter')xlabel('Frequencyinrad/s')ylabel('Phaseindegree;ModuleindB')figure(6)semilogy(0:N-n,tr(n:N)),gridtitle('Cxxmatrixtrace')xlabel('Samples')实验结果图如下：图4，输入信号x(k)图5参考信号xr(k)图6误差e(k)和输出信号y(k)图7.滤波器系数a(n-1)和a(n-1)变化曲线系数的变化曲线在200步时有一个超调，这是由于h(k)向量为零，所以200步以后仅代表x值。获得的滤波器的传递函数也类似于LMS滤波器的传递函数，相应的预测也类似。它的中心频率调整为正弦信号频率，即75rad/s，如下图所示图8.合成滤波器传递函数的幅频特性和相频特性图9Cxx矩阵曲线Cxx曲线在采样步数n=200时突变。200步以后曲线值变小，使滤波器不能再根据输入信号的统计变换进行调整。