自适应噪声抵消LMS算法Matlab仿真

自适应噪声抵消LMS算法Matlab仿真传统的宽带信号中抑制正弦干扰的方法是采用陷波器(notchfilter),为此我们需要精确知道干扰正弦的频率.然而当干扰正弦频率是缓慢变化时,且选频率特性要求十分尖锐时,则最好采用自适应噪声抵消的方法.下图是用一个二阶FIR的LMS自适应滤波器消除正弦干扰的一个方案。1）借助MATLAB画出误差性能曲面和误差性能曲面的等值曲线；2）写出最陡下降法，LMS算法的计算公式（δ=0.4）；3）用MATLAB产生方差为0.05,均值为0白噪音S(n)，并画出其中一次实现的波形据2）中的公式，并利用3）中产生的S(n)，在1）中的误差性能曲面的等值曲n的值曲线上叠加画出LMS法时100情况确定，一般选取足够大以使算法达到基)(ny宽带信号+正弦干扰0()()()ynSnNn=+图；4）根线上叠加画出采用最陡下降法，LMS法时H(n)的在叠代过程中的轨迹曲线。5）用MATLAB计算并画出LMS法时随时间变化曲线（对应S(n)的某一次的一次实现）和e(n)波形；某一次实现的结果并不能从统计的角度反映实验的结果的正确性，为得到具有统计特性的实验结果，可用足够多次的实验结果的平均值作为实验的结果。用MATLAB计算并画出LMS法时J(n)的100次实验结果的平均值随时间n的变化曲线。6）用MATLAB计算并在1）中的误差性能曲面的等次实验中的H(n)的平均值的轨迹曲线；（在实验中n=1,,…..N,N的取值根据实验本收敛）01(),(0)0.052()sin()16102()2sin()16ssSnrNnnNnnπππ==+=是均匀分布的白噪音不相关和)(),()(10nNnNnS⊕ADF)(nx)(~ny)(nyx1()())(ne正弦干扰相关成分nNn=FIR阶次为2(N=2)−WidebandSignal1、用Matlab画误差性能曲面和误差性能曲面的等值曲线的程序如下：[h0,h1]=meshgrid(-2:0.1:4,-4:0.1:2);J=0.55+h0.*h0+h1.*h1+2*cos(pi/8)*h1.*h0-sqrt(2)*h0*cos(pi/10)-sqrt(2)*h1*cos(9*pi/40);echoon;v=0:0.1:2;%axis([-44-440100]);figure(1);surf(h0,h1,J);%误差曲面xlabel('h0');ylabel('h1');title('误差性能曲面');figure(2);contour(h0,h1,J,v);%等值曲线xlabel('h0');ylabel('h1');title('误差性能曲面等值曲线');运行结果如下图示：2、①最陡下降法计算公式：）（n21)()1(HGVnHnδ−=+其中δ取0.4，H(0)=[3-4]，T⎟⎠⎞⎜⎝⎛+=⎟⎠⎞⎜⎝⎛−⎥⎦⎤⎢⎣⎡⎟⎠⎞⎜⎝⎛+==⎟⎠⎞⎜⎝⎛−⎟⎠⎞⎜⎝⎛=⎥⎦⎤⎢⎣⎡−⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=−=∑∑==1016k2cos2116)(2sin210162sin2161)(r16k2cos16)(2sin2162sin2161)(r)1()0(2)()()0()1()1()0(22)(2)(V15015010GππππππππiyxixxyxyxxxxxxxxxyxxxkiikkiikrrnhnhrrrrrnHRn而故⎥⎦⎤⎢⎣⎡−−=⎥⎦⎤⎢⎣⎡−⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=5320.37362.2)0(5377.06725.02)()(19239.09239.012)(10GGVnhnhnV②LMS算法计算公式：,...2,1,0),1()1()()1()1()()1()1(e=+++=++−+=+nnXnenHnHnXnHnynTδ其中δ取0.4。3、用MATLAB产生方差为0.05,均值为0白噪音s(n)程序如下：n=1:1024;s(n)=sqrt(0.05)*randn(1,length(n));plot(s);axis([01024-11]);title('白噪声s(n)');其一次实现图形如下：4、最陡下降法和LMS算法程序：clearall;N=1000;%信号点数k=2;%滤波器阶数q=0.4;%步长s=sqrt(0.05)*randn(1,N);%原始信号n=1:N;y=s(n)+sin(2*pi*n/16+pi/10);%叠加噪声后的期望输出信号x=sqrt(2)*sin(2*pi*n/16);%相关噪声输入figure(1);subplot(3,1,1);plot(n,s);title('原始宽带信号S(n)');subplot(3,1,2);plot(n,y);title('叠加噪声后输出信号y(n)');subplot(3,1,3);plot(n,x);title('相关噪声输入信号x(n)');%设置初值并进行两种算法的迭代hn=zeros(2,N);%存放最陡下降法下H(n)迭代数据hn(:,1)=[3-4]';vg=[-2.7362-3.5320]';Rxx=[[10.9239]',[0.92391]'];Ryx=[0.67250.5377]';echooff;fori=1:N-1%最陡下降法，第n=i-1次迭代vg=2*Rxx*hn(:,i)-2*Ryx;hn(:,i+1)=hn(:,i)-1/2*q*vg;end;hm=zeros(2,N);%存放LMS算法下H(n)迭代数据h=[3-4]';hm(:,1)=h;fori=2:N%XN=x(i+1:i)';u=x(i:-1:i-1);en=y(i)-u*h;h=h+q*en*u';hm(:,i)=h;end;%绘制图像[h0,h1]=meshgrid(-2:.1:4,-4:.1:2);J=0.55+h0.*h0+h1.*h1+2*cos(pi/8)*h1.*h0-sqrt(2)*h0*cos(pi/10)-sqrt(2)*h1*cos(9*pi/40);v=0:0.1:2;%等值曲线取值figure(2);%surf(h0,h1,J);%误差曲面contour(h0,h1,J,v);%等值曲线xlabel('h0');ylabel('h1');holdon;plot(hn(1,:),hn(2,:),'g');plot(hm(1,:),hm(2,:),'r');title('最陡下降法，LMS法时H(n)的在叠代过程中的轨迹曲线');legend('等值曲线','最陡下降法轨迹','LMS算法1次轨迹','NorthEastOutside');图形如下：5、单次实现及统计平均(100次试验)下J(n)及e(n)特性程序如下：clearall;g=100;%统计平均次数N=400;k=2;%滤波器阶数q=0.4;enn=zeros(1,N-k+1);%存放单次实现下的J(n)=e(n)^2en0=zeros(1,N-k+1);%存放单次实现下的e(n)ejn=zeros(1,N-k+1,g);%存放100次实现下的J(n)=e(n)^2%g次统计试验forj=1:gs=sqrt(0.05)*randn(1,N);%原始信号n=1:N;y=s(n)+sin(2*pi*n/16+pi/10);%叠加噪声后的期望输出信号x=sqrt(2)*sin(2*pi*n/16);%相关噪声输入%设置初值，并进行一次LMS算法迭代h=[3-4]';fori=k:Nu=x(i:-1:i-1);en=y(i)-u*h;h=h+q*en*u';en0(i-k+1)=en;%保存单次e(n)，由于进行g次试验，仅保存最后一次enn(i-k+1)=en^2;%保存单次J(n)，由于进行g次试验，仅保存最后一次ejn(:,i-k+1,j)=en^2;%保存g次试验的J(n)end;end;%g次试验结束%计算g次试验J(n)的统计平均值fori=1:N-k+1jn(i)=sum(ejn(1,i,:))/g;end;figure(1);subplot(2,1,1);plot(enn);%单次实现下Jntitle('LMS算法单次实现下J(n)');subplot(2,1,2);plot(en0,'r');title('LMS算法单次实现下e(n)');figure(2);plot(jn);%统计平均下J(n)title('LMS算法100次实现下平均J(n)');单次实现下J(n)和e(n)如下图所示：100次统计平均下的J(n)如下图所示：6、LMS算法下100次实验H(n)曲线程序如下：clearall;g=100;%统计试验次数N=1000;%信号点数k=2;%滤波器阶数q=0.4;hh=zeros(2,N,g);%存放g次迭代下LMS算法的H(n)forj=1:g%g次统计试验s=sqrt(0.05)*randn(1,N);%原始信号n=1:N;y=s(n)+sin(2*pi*n/16+pi/10);%叠加噪声后的期望输出信号x=sqrt(2)*sin(2*pi*n/16);%相关噪声输入%设置初值hn=zeros(2,N);%存放最陡下降法下H(n)迭代数据，仅保存最后一次试验hn(:,1)=[3-4]';vg=[-2.7362-3.5320]';Rxx=[[10.9239]',[0.92391]'];Ryx=[0.67250.5377]';echooff;fori=1:N-1%最陡下降法，第n=i-1次迭代vg=2*Rxx*hn(:,i)-2*Ryx;hn(:,i+1)=hn(:,i)-1/2*q*vg;end;%LMS算法一次试验hm=zeros(2,N);%存放LMS算法下H(n)迭代数据h=[3-4]';hm(:,1)=[3-4]';fori=k:N%迭代N-k+1次u=x(i:-1:i-1);en=y(i)-u*h;h=h+q*en*u';hm(:,i)=h;hh(:,i,j)=h;%保存g次试验下H(n)end;end;fori=1:N-k+1%统计平均hhh(1,i+1)=sum(hh(1,i+1,:))/g;hhh(2,i+1)=sum(hh(2,i+1,:))/g;end;hhh(:,1)=[3-4]';[h0,h1]=meshgrid(-2:.1:4,-4:.1:2);J=0.55+h0.*h0+h1.*h1+2*cos(pi/8)*h1.*h0-sqrt(2)*h0*cos(pi/10)-sqrt(2)*h1*cos(9*pi/40);v=0:0.1:2;figure(1);contour(h0,h1,J,v);%等值曲线xlabel('h0');ylabel('h1');holdon;plot(hn(1,:),hn(2,:));%最陡下降法轨迹曲线plot(hhh(1,:),hhh(2,:),'r');%LMS法时H(n)100次统计平均轨迹曲线title('最陡下降法及LMS法时H(n)100次统计平均轨迹曲线');legend('等值曲线','最陡下降法轨迹','LMS算法100次平均轨迹','NorthEastOutside');图形如下：