统计信号处理实验四东南大学

统计信号处理实验四《统计信号处理》实验四目的：掌握自适应滤波的原理；内容一：假设一个接收到的信号为：x(t)=s(t)+n(t),其中s(t)=A*cos(wt+a),已知信号的频率w=1KHz，而信号的幅度和相位未知，n(t)是一个服从N(0,1)分布的白噪声。为了利用计算机对信号进行处理，将信号按10KHz的频率进行采样。1)通过对x(t)进行自适应信号处理，从接收信号中滤出有用信号s(t);2）观察自适应信号处理的权系数；3）观察的滤波结果在不同的收敛因子u下的结果，并进行分析；4）观察不同的抽头数N对滤波结果的影响，并进行分析；内容二：在实验一的基础上，假设信号的频率也未知，重复实验一；内容三：假设s(t)是任意一个峰峰值不超过1的信号（取幅度为0.5的方波），n(t)是一个加在信号中的幅度和相位未知的，频率已知的50Hz单频干扰信号（假设幅度为1）。信号取样频率1KHz，试通过自适应信号处理从接收信号中滤出有用信号s(t)。要求：1）给出自适应滤波器结构图；2）设计仿真计算的Matlab程序，给出软件清单；3）完成实验报告，对实验过程进行描述，并给出试验结果，对实验数据进行分析。实验过程：1、假设一个接收到的信号为：d(t)=s(t)+n(t),其中s(t)=A*cos(wt+a),已知信号的频率w=1KHz，而信号的幅度和相位未知，n(t)是一个服从N(0,1)分布的白噪声。为了利用计算机对信号进行处理，将信号按10KHz的频率进行采样。1）参考信号d(k)=s(k)+n(k)，s(k)=A*cos(wk+a),产生一个与载波信号具有相同频率的正弦信号作为输入信号()xk,即x(k)=cos(wk)。经过自适应处理后，就可以在输出信号()yk端得到正确的载波信号（包含相位和幅度）。框图如下：Σ与载波同频正弦信号d(k)＝s(k)+n(k)+-自适应滤波器e(k)y(k)x(k)s^(k))()()()1()()()()()()()()()(kxkekwkwkxkwkdkykdkekxkwkyTT2）改变收敛因子，观察滤波结果。3）改变滤波器抽头数N，观察滤波结果。2、在实验一的基础上，假设信号的频率也未知，重复实验一。参考信号d(k)=s(k)+n(k)，s(k)=A*cos(wk+a),将参考信号延时一段时间后得到的信号作为输入信号()xk,即x(k)=d(k-m)。经过自适应处理后，就可以在误差输出端y(k)得到正确的载波信号（包含频率、相位和幅度）。3、假设s(t)是任意一个峰峰值不超过1的信号（取幅度为0.5的方波），n(t)是一个加在信号中的幅度和相位未知的，频率已知的50Hz单频干扰信号（可以假设幅度为1）。信号取样频率1KHz，试通过自适应信号处理从接收信号中滤出有用信号s(t)。我们可以使用陷波滤波器对噪声进行滤除，但普通滤波器一旦做成，其陷波频率难以调整。如果使用自适应陷波滤波器，不仅可以消除单频干扰，而且可以跟踪干扰的频率变化，持续消噪。自适应陷波滤波器的原理框图如下图所示：参考信号d(t)输入信号Ccos(ω0t+φ)90移相器Σd(k)x1(k)x2(k)w1(k)w2(k)Σ+-y1(k)y2(k)y(k)LMS算法+e(k)同步取样假如输入信号是一个纯余弦信号0cosCt，则可将其分为两路，将其中一路进行90度相移，然后同步采样，得到：1020cossinxkCkTxkCkT它们通过相关抵消回路以最小均方算法去控制1wk和2wk加权，然后两个加权输出1yk和2yk相加得到yk，称为自适应滤波器的输出。权值计算如下：1112221212wkwkekxkwkwkekxk写成矩阵形式为：(1)()2()()kkekkwwx也可以简单的采用与内容1一样的LMS滤波器，本实验选择了只有一路参考信号的自适应滤波器。实验结果及分析：1.已知信号频率。可以看出，滤波器输出信号（蓝色波形）与有用信号（红色波形）逐渐逼近，加权系数也逐渐稳定。改变收敛因子：可以看出，收敛因子太小会导致收敛速度过慢，收敛因子太大会导致最终波动幅度太大，无法精确地逼近有用信号。改变滤波器阶数：可以看出，阶数越高，收敛地越快，但更高得阶数需要适当改变收敛因子来达到最理想的滤波效果。2.信号频率未知。可以看出，输出信号可以逼近有用信号s(t)，但性能不如第一题的结果那么好。改变收敛因子：可以看出，收敛因子太小会导致收敛速度过慢，收敛因子太大会导致最终波动幅度太大，无法精确地逼近有用信号。改变滤波器阶数：可以看出，阶数越高，收敛地越快，但更高得阶数需要适当改变收敛因子来达到最理想的滤波效果。3.消除干扰采用30阶滤波器，收敛因子μ取0.0015，由输出波形可以看出，该滤波器消除了频率干扰，得到了矩形波信号。源程序：%lms.mfunction[e,wl,y]=lms(N,u,n,x,d)wl=zeros(N,n);e=zeros(1,n);y=zeros(1,n);fori=N:n-1X=x(i:-1:i-N+1);y(i)=X*wl(:,i);%i时刻输出信号e(i)=d(i)-y(i);%i时刻误差信号wl(:,i+1)=wl(:,i)+(u.*e(i).*X)';%i时刻滤波器的权值end;%1clearall;n=500;%信号点数f=1000;fs=10000;%采样点数N=1;%滤波器阶数t=(0:n-1)/fs;s=cos(2*pi*f*t+pi/8);%有用信号nt=randn(1,n);%噪声信号x=cos(2*pi*f*t);%参考信号d=s+nt;%观测信号u=0.016;%LMS算法下自适应增益常数[~,wl,y]=lms(N,u,n,x,d);figure(1);subplot(3,1,1);plot(t,d);grid;title('观测信号d');subplot(3,1,2);plot(t,y,'b',t,s,'r');grid;title('输出信号y');subplot(3,1,3);plot(t,wl);grid;title('加权因子w');u=[0.0050.0150.1];figure(2);[~,~,y]=lms(N,u(1),n,x,d);subplot(3,1,1);plot(t,y,'b',t,s,'r');grid;title('u=0.005');[~,~,y]=lms(N,u(2),n,x,d);subplot(3,1,2);plot(t,y,'b',t,s,'r');grid;title('u=0.015');[~,~,y]=lms(N,u(3),n,x,d);subplot(3,1,3);plot(t,y,'b',t,s,'r');grid;title('u=0.1');%改变阶数u=0.016;N=[151015];figure(3);[~,~,y]=lms(N(1),u,n,x,d);subplot(4,1,1);plot(t,y,'b',t,s,'r');grid;title('N=1');[~,~,y]=lms(N(1),u,n,x,d);subplot(4,1,2);plot(t,y,'b',t,s,'r');grid;title('N=5');[~,~,y]=lms(N(1),u,n,x,d);subplot(4,1,3);plot(t,y,'b',t,s,'r');grid;title('N=10');[~,~,y]=lms(N(1),u,n,x,d);subplot(4,1,4);plot(t,y,'b',t,s,'r');grid;title('N=15');%2clearall;n=600;%信号点数f=1000;fs=10000;t=(0:n-1)/fs;s=cos(2*pi*f*t+pi/8);nt=randn(1,n);d=s+nt;d1=[dzeros(1,10)];%将观测信号补零x=[zeros(1,10)d];%延时后的信号作为参考信号u=0.001;%LMS算法下自适应增益常数N=12;m=n+10;[~,wl,y]=lms(N,u,m,x,d1);figure(1);subplot(3,1,1);plot(t,d);grid;title('观测信号');t=(0:m-1)/fs;s=[szeros(1,10)];subplot(3,1,2);plot(t,y,'b',t,s,'r');grid;title('输出信号');ylim([-1.5,1.5]);xlim([0,0.06])subplot(3,1,3);plot(t,wl);grid;title('加权因子');xlim([0,0.06])%改变收敛因子u=[0.00020.0010.01];figure(2);[~,~,y]=lms(N,u(1),m,x,d1);subplot(3,1,1);plot(t,y,'b',t,s,'r');grid;ylim([-1.5,1.5]);xlim([0,0.06]);title('u=0.0002');[~,~,y]=lms(N,u(2),m,x,d1);subplot(3,1,2);plot(t,y,'b',t,s,'r');grid;ylim([-1.5,1.5]);xlim([0,0.06]);title('u=0.001');[~,~,y]=lms(N,u(3),m,x,d1);subplot(3,1,3);plot(t,y,'b',t,s,'r');grid;ylim([-1.5,1.5]);xlim([0,0.06]);title('u=0.01');%改变阶数u=0.001;N=[11030];figure(3);[~,~,y]=lms(N(1),u,m,x,d1);subplot(3,1,1);plot(t,y,'b',t,s,'r');grid;ylim([-1.5,1.5]);xlim([0,0.06]);title('N=1');[~,~,y]=lms(N(2),u,m,x,d1);subplot(3,1,2);plot(t,y,'b',t,s,'r');grid;ylim([-1.5,1.5]);xlim([0,0.06]);title('N=10');[~,~,y]=lms(N(3),u,m,x,d1);subplot(3,1,3);plot(t,y,'b',t,s,'r');grid;ylim([-1.5,1.5]);xlim([0,0.06]);title('N=30');%3clearall;n=1001;%信号点数N=30;%阶数f=50;%t=1:0.01:n;T=5;t=-5*T:0.05:5*T;s=0.5*square(t,50);%产生周期性方波信号nt=sin(f*t);%工频干扰x=sin(f*t);%参考信号d=s+nt;%观测信号u=0.0015;%LMS算法下自适应增益常数[e,wl,y]=lms(N,u,n,x,d);figure(1);subplot(3,1,1);plot(t,d);grid;title('观测信号d');subplot(3,1,2);plot(t,e,'b',t,s,'r');grid;title('输出信号y');ylim([-1.5,1.5]);subplot(3,1,3);plot(t,wl);grid;title('加权因子w');