LS.rar_LS加窗算法_加窗LS_滤波协方差_滤波算法比较_相关滤波

%%%%%%%%%%%%%%%%%%%%%%%%%LS算法%%%%%%%%%%%%%%%%%%%% %%%刘建卫 %%%2009.6.9 N=2000; M=1000; interval=4*pi; n=0:interval/N:interval-interval/N; vn=randn(1,N); yn=(1+cos(2*pi*n*0.5)).*cos(2*n*pi*2); xn=yn+vn; x=[xn,zeros(1,M)]; A=[;;]; for i=1:N+M-1 for j=1:M if i>=j A(i,j)=x(i-j+1); else A(i,j)=0; end end end figure plot(n,yn),title('原始信号波形') grid on figure plot(n,x(1:N)),title('叠加噪声以后的波形') grid on %%%%%%%自相关算法%%%%%%%%%%% h=zeros(1,M); Acovariance=[;;]; for i=M:N; for j=1:M Acovariance(i-M+1,j)=A(i,j); end end AHcovariance=Acovariance'; d=[]; for i=M:N d(i+1-M)=x(i)-vn(i); end h=inv(AHcovariance*Acovariance)*AHcovariance*d'; xLSc=conv(x,h); figure plot(n,xLSc(1:N)),title('用自相关算法得到的滤波波形') grid on %%%%%%%%协方差算法%%%%%%%%%%%%%%%%%%% hs=zeros(1,M); Aselfcorrelation=A; ds=zeros(1,N+M-1); for i=1:N ds(i)=x(i)-vn(i); end AHselfcorrelation=Aselfcorrelation'; hs=inv(AHselfcorrelation*Aselfcorrelation)*AHselfcorrelation*ds'; xLSs=conv(x,hs); figure plot(n,xLSs(1:N),'r',n,yn,'.-'),title('N=1000,M=500时的滤波波形') grid on legend('LS滤波以后的波形','原始信号波形') figure plot(n,xLSs(1:N)-yn) grid on title('误差信号波形') %%%%%%%%%%%%%%%%%%后向加窗算法%%%%%%%%%%%%%%%%%%%%%%% hba=zeros(1,M); dba=zeros(1,N); Abaw=[;;]; for i=M:N+M-1 for j=1:M Abaw(i-M+1,j)=A(i,j); end end AHbaw=Abaw'; for i=M:N dba(i-M+1)=yn(i); end hba=inv(AHbaw*Abaw)*AHbaw*dba'; xLSba=conv(hba,x); figure plot(n,xLSba(1:N)),title('用后向加窗算法滤波以后的波形') grid on %%%%%%%%%%%前加窗算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% hfa=zeros(1,M); dfa=[]; Afaw=[;;]; for i=1:N for j=1:M Afaw(i,j)=A(i,j); end end AHfaw=Afaw'; for i=1:N dfa(i)=yn(i); end hfa=inv(AHfaw*Afaw)*AHfaw*dfa'; xLSfa=conv(hfa,x); figure plot(n,xLSfa(1:N)),title('用前向加窗算法滤波以后的波形') grid on %%%%%%几种滤波方法的比较%%%%%%%%%%%%%%%% figure plot(n,yn,'r',n,xLSc(1:N),n,xLSs(1:N),':m',n,xLSba(1:N),'--k',n,xLSfa(1:N),'-.g'),title('几种加窗方法的比较') grid on legend('原始信号波形','自相关算法','协方差算法','后向加窗算法','前向加窗算法') %%%%%%%%%%前向加窗的QR解法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [Q,R]=qr(Afaw); Rfa=[;;]; Qfa=[;;]; for i=1:M for j=1:M Rfa(i,j)=R(i,j); end end for i=1:N; for j=1:M Qfa(i,j)=Q(i,j); end end QHfa=Qfa'; zfa=QHfa*yn'; hQRfa=inv(Rfa)*zfa; xLSQRfa=conv(hQRfa,x); figure plot(n,xLSQRfa(1:N)),title('用前向加窗算法QR分解法滤波以后的波形') grid on %%%%%%%%%%%%%%两种解法的比较%%%%%%%%%%%%%%%%%%%%% figure plot(n,xLSfa(1:N),'-*',n,xLSQRfa(1:N),'k'),title('两种解法的比较') grid on legend('常规方程算法','QR分解算法') %%%%%%%%%%%%%%%%%the end%%%%%%%%%%%%%%%%%%