wiener.m

matlab编程

function wiener(a,db,L,N)%db是信噪比，N是AR模型阶数，L是信号长度
bs=10^(db/10);%信号与噪声的比值
w=randn(1,L);%产生均值为0的随机信号
s=zeros(1,L);
x=zeros(1,L);
r_ss=zeros(1,L);% r_ss是s的相关函数，r_ss（1）其实是r_ss(0),是方差
r_vv=zeros(1,L);
r_xx=zeros(1,L);
r_xs=zeros(1,L);
s(1)=w(1);%s(0)=0
for i=2:L
    s(i)=a*s(i-1)+w(i);
end
for i=1:L
    r_ss(1)=r_ss(1)+s(i)*s(i);
end
r_ss(1)=r_ss(1)/L;
r_vv(1)=r_ss(1)/bs;% r_vv(1)是白噪声信号的方差
v=sqrt(r_vv(1))*randn(1,L);%要产生均值为0且方差等于r_vv(0)的白噪声
x(1:L)=s(1:L)+v(1:L)% x(n)是观测数据也是输入数组
r_xx(1)=r_ss(1)+r_vv(1);% r_xx是x(n)的自相关函数
r_xs(1)=r_ss(1);% r_xs是x(n)与s(n)的互相关函数
for i=2:N
    %r_ss(i)=a^(1-i)*r_ss(1);
     r_ss(i)=a^(i-1)*r_ss(1);
    r_vv(i)=0;
    r_xx(i)=r_ss(i)+r_vv(i);
    r_xs(i)=r_ss(i);
end
R_xs=zeros(N,1);
R_xx=zeros(N,N);
for i=1:N
    k=1;
    for j=i:N
        R_xx(i,j)=r_xx(k);
        k=k+1;
    end
   if i>1
       k=2;
      for j=i-1:1
          R_xx(i,j)=r_xx(k);
          k=k+1;
      end
  end
  R_xs(i)=r_xs(i);
end
hopt=zeros(N,1);
hopt=inv(R_xx)*R_xs
%E=r_ss(1)-(conj(R_xs))'*hopt
Emin=abs(r_ss(1)-(conj(R_xs))'*hopt)

y=conv(hopt',x);
n=1:L;
plot(n,s,'-bo');
hold on;
plot(y,'-r');
hold off;
%xlabel('n'),ylabel('观测数据x(n)');
%grid on;
%title('观测数据x的波形');