wiener filter.txt

通过计算机MATLAB仿真分析和研究了维纳滤波器的阶数、信号的噪声方差、随机信号的采样点数、经过维纳滤波的均方误差之间的关系

clc;clear;
%初始化变量
L=300;N1=150; %滤波器阶数
m=3;    %噪声方差
x1=randn(1, L);x1=x1/std(x1);x1=x1-mean(x1);
wn=0+sqrt(1)*x1;wn=wn';
sn=zeros(L,1);
for i=1:1:51,
   xn(:,i)=zeros(L,1);
end
x=zeros(L,1);
%根据给定公式生成随机信号sn
sn(1,1)=wn(1,1); sn(2,1)=wn(2,1)+0.5833*sn(1,1);sn(3,1)=wn(3,1)+0.5833*sn(2,1)+0.3750*sn(1,1);
for i=4:1:L,
sn(i,1)=0.5833*sn(i-1,1)+0.3750*sn(i-2,1)-0.04167*sn(i-3,1)+wn(i,1); 
end
%原始信号图像
subplot(2,2,1);plot(sn),axis([200 320 -10 10]),xlabel('时间'),ylabel('幅度'),title('原始随机信号时间谱');
x2=randn(1,L); x2=x2/std(x2); x2=x2-mean(x2); a=0; 
for p=1:1:51,
	b(1,p)=sqrt(p); %噪声方差
    A=a+b(1,p)*x2;
	vn(:,p)=A';%生成高斯噪声信号	
	xn(:,p)=sn+vn(:,p);%受噪后信号
   end  
%高斯白噪声信号图像    
subplot(2,2,2);plot(vn(:,m));axis([200 320 -10 10]);xlabel('时间'),ylabel('幅度'),title('高斯白噪声时间谱');
%受噪后信号和原始信号图像对比
subplot(2,2,3);plot(xn(:,m));axis([200 320 -10 10]);hold on;
plot(sn,'r');axis([200 320 -10 10]);xlabel('时间'),ylabel('幅度'),title('红:原始随机信号 蓝:受噪后信号');
%生成滤波器系数h(n)
for N=1:1:N1,
  hn(:,N)=zeros(N1,1);	P1=zeros(N1,N1);
end
for N=1:1:N1,
	[W1,R,P1(N,[1:N+1])]=firwiener(N,xn(:,m),sn);	hn([1:N+1],N)=W1';
end
%实现滤波
for N=1:1:N1
    y(:,N)=conv(xn(:,m),hn(:,N));
end
%画出理想信号与输出信号对比图
subplot(2,2,4);
plot(sn,'r'),axis([0 200 -10 10]),xlabel('时间'),ylabel('幅度'),title('输入与输出信号对比');hold on;
plot(y(:,N1),'g'),axis([200 320 -10 10]);hold on;
plot(xn(:,m),'b'),axis([200 320 -10 10]);
title('红:原始随机信号 蓝:受噪后信号(即滤波器输入信号) 绿:维纳滤波输出');
j=zeros(N1,1);s=std(sn);s2=s*s; %原始随机信号的方差
for N=1:1:N1,
	j(N,1)=s2-P1(N,[1:1+N])*hn([1:1+N],N);  %不同阶数对应的均方误差
end
%滤波器阶数从1到160变化时对应均方误差与阶数的关系图
n=(1:1:N1);
figure;plot(n,j');xlabel('滤波器阶数'),ylabel('均方误差'),title('均方误差与滤波器阶数的关系图');
e=zeros(L,N1);
for N=1:1:N1,
    for i=1:1:L,
    e(i,N)=y(i,N)-sn(i,1);  %误差
    end 
end
figure;plot(e(:,N1));  title('阶数为N1时的滤波误差');
%滤波器阶数为50时高斯白噪声方差从1到51时均方误差与之的关系
P2=zeros(51,51);
for i=1:1:51,
    [W2,R2,P2(i,:)]=firwiener(50,xn(:,i),sn);		hn1(:,i)=W2';
end
for i=1:1:51
    y1(:,i)=conv(xn(:,i),hn1(:,i));
end
figure;plot(y1(:,30),'g'),axis([200 320 -20 20]);hold on;
plot(sn,'r'),axis([0 200 -20 20]),xlabel('时间'),ylabel('幅度');hold on;
plot(xn(:,30),'b'),axis([200 320 -20 20]);title('滤波器阶数为50噪声方差为30时的效果图');
j1=zeros(51,1);
for i=1:1:51,
	j1(i,1)=s2-P2(i,:)*hn1(:,i);
end
n1=(1:1:51);
figure;plot(n1,j1');title('噪声方差与均方误差的关系图');  xlabel('噪声方差'),ylabel('均方误差');
e1=zeros(L,51);
for i=1:1:51,
    for j=1:1:L,
        e1(j,i)=y1(j,i)-sn(j,1);
    end
end
figure;plot(e1(:,4));title('噪声方差为4时的滤波误差');