自适应噪声抵消LMS 算法Matlab 仿真传统的宽带信号中抑制正弦干扰的方法是采用陷波器(notch filter),为此我们需要精确知道干扰正弦的频率.然而当干扰正弦频率是缓慢变化时,且选频率特性要求十分尖锐时,则最好采用自适应噪声抵消的方法.下图是用一个二阶FIR 的LMS 自适应滤波器消除正弦干扰的一个方案。

1） 借助MATLAB 画出误差性能曲面和误差性能曲面的等值曲线； 2） 写出最陡下降法， LMS 算法的计算公式（δ=0.4）；3） 用MATLAB 产生方差为0.05,均值为0白噪音S(n)，并画出其中一次实现的波形据2）中的公式，并利用3）中产生的S(n)，在1）中的误差性能曲面的等值曲n 的值曲线上叠加画出LMS 法时100情况确定，一般选取足够大以使算法达到基)(n y 宽带信号+正弦干扰0()()()y n S n N n =+图；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()sin()16ss S n r N n n N n n πππ==+是均匀分布的白噪音不相关和)(),()(10n N n N n S)(n x x 1()())(n e n N n =1、用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);echo on;v=0:0.1:2;%axis([-4 4 -4 4 0 100]);figure(1);%误差曲面surf(h0,h1,J);xlabel('h0');ylabel('h1');title('误差性能曲面');figure(2);contour(h0,h1,J,v); %等值曲线xlabel('h0');ylabel('h1');title('误差性能曲面等值曲线');运行结果如下图示：2、①最陡下降法计算公式：）（n 21)()1(H G V n H n δ−=+ 其中δ取0.4，H(0)=[3 -4]，T ⎟⎠⎞⎜⎝⎛+=⎟⎠⎞⎜⎝⎛−⎥⎦⎤⎢⎣⎡⎟⎠⎞⎜⎝⎛+==⎟⎠⎞⎜⎝⎛−⎟⎠⎞⎜⎝⎛=⎥⎦⎤⎢⎣⎡−⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=−=∑∑==1016k 2cos 2116)(2sin 210162sin 2161)(r 16k2cos 16)(2sin 2162sin 2161)(r )1()0(2)()()0()1()1()0(22)(2)(V 15015010G ππππππππi yx i xx yx yx xx xx xx xx yxxx k i i k k i i k r r n h n h r r r r r n H R n 而故⎥⎦⎤⎢⎣⎡−−=⎥⎦⎤⎢⎣⎡−⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=5320.37362.2)0(5377.06725.02)()(19239.09239.012)(10G G V n h n h n V②LMS 算法计算公式：,...2,1,0),1()1()()1()1()()1()1(e =+++=++−+=+n n X n e n H n H n X n H n y n T δ其中δ取0.4。

3、 用MATLAB 产生方差为0.05,均值为0白噪音s(n) 程序如下： n=1:1024;s(n)=sqrt(0.05)*randn(1,length(n)); plot(s);axis([0 1024 -1 1]); title('白噪声s(n)');其一次实现图形如下：4、 最陡下降法和LMS 算法程序： clear all;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=[[1 0.9239]',[0.9239 1]'];Ryx=[0.6725 0.5377]';echo off;for i=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;for i=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/4 0);v=0:0.1:2;%等值曲线取值figure(2);%surf(h0,h1,J); %误差曲面contour(h0,h1,J,v); %等值曲线xlabel('h0');ylabel('h1');hold on;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)特性程序如下：clear all;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次统计试验for j=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]';for i=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)的统计平均值for i=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)曲线程序如下：clear all;g=100;%统计试验次数N=1000;%信号点数k=2;%滤波器阶数q=0.4;hh=zeros(2,N,g);%存放g次迭代下LMS算法的H(n)for j=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=[[1 0.9239]',[0.9239 1]'];Ryx=[0.6725 0.5377]';echo off;for i=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]';for i=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;for i=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/4 0);v=0:0.1:2;figure(1);contour(h0,h1,J,v); %等值曲线xlabel('h0');ylabel('h1');hold on;plot(hn(1,:),hn(2,:));%最陡下降法轨迹曲线plot(hhh(1,:),hhh(2,:),'r');%LMS法时H(n)100次统计平均轨迹曲线title('最陡下降法及 LMS法时H(n)100次统计平均轨迹曲线');legend('等值曲线','最陡下降法轨迹','LMS算法100次平均轨迹','NorthEastOutside');图形如下：。