wienerfilter.m

设计一维纳滤波器,估计AR模型参数以及无码率

function WienerFilter
%第二章上机作业二：设计Wiener滤波器
%信研0501班 李鹏
clc;
clear all;
a=0.6;
input1=inputdlg('请选择信号长度L=(32,256,1024)');%选择输入信号长度和模型阶数
L=str2num(char(input1));
SN1=6;                                         %不同的信噪比
SN2=10;
SN3=20;
W=random('norm',0,1,L,1);                    %产生高斯白噪声序列(L*1)
S(1)=0;
for n=2:L
      S(n)=a*S(n-1)+W(n);                    %产生原始信号
  end
  Am=sum(abs(S).^2)/L;
  P1=Am/(10^(SN1/20));                     %根据信噪比计算噪声功率，也就是高斯噪声的方差
  P2=Am/(10^(SN2/20));
  P3=Am/(10^(SN3/20));
  V1=random('norm',0,P1,L,1);
  V2=random('norm',0,P2,L,1);
  V3=random('norm',0,P3,L,1);
  for n=1:L                                %产生三种观测信号
      X1(n)=S(n)+V1(n);
      X2(n)=S(n)+V2(n);
      X3(n)=S(n)+V3(n);
  end
  fprintf('AR模型阶数N依次取1,2,3：\n')
i=1:1:L;                                    %绘制信号波形图
plot(i,X1(i),'g',i,X2(i),'b',i,X3(i),'r');
title('三种不同信噪比下的信号波形对比');
xlabel('信号序列点');
ylabel('幅度');
legend('S/N=6dB','S/N=10dB','S/N=20dB');
grid on;
fprintf('\n对X1信号来说：\n')
 AR(X1,1);                              %产生N阶AR模型参数a,和误差e;
 AR(X1,2); 
 AR(X1,3);   
 fprintf('\n对X2信号来说：\n')
 AR(X2,1);                           
 AR(X2,2); 
 AR(X2,3);
fprintf('\n对X3信号来说：\n')
 AR(X3,1);                              
 AR(X3,2); 
 AR(X3,3);

 function AR(X,N)                         %子函数实现AR模型；
for m=1:N                              %产生Yule-Walker方程矩阵形式
    for n=m:N                    
        rx(m,n)=Rxx(X,n-m);
        rx(n,m)=rx(m,n);
    end
    rx(m,m)=Rxx(X,0);
    zx(m)=Rxx(X,m);
end
ap=inv(rx)*(-zx');
a=[1,ap'];
e=Rxx(X,0)+zx*ap;
disp(['AR系数ap： ',num2str(a)]);
disp(['均方误差e：',num2str(e)]);



function Rxx=Rxx(X,m)                       %生成Rxx自相关矩阵
N=length(X);
for n=1:N-abs(m)
    s(n)=X(n)*X(n+m);
end
Rxx=sum(s)/N;