wiener.m

研究生时的现代信号处理作业

clear;
clc;
N=200;
v=randn(1,N);
a1=-1.352;a2=1.338;a3=-0.662;a4=0.240;
x(1)=v(1);
x(2)=-a1*x(1)+v(2);
x(3)=-a1*x(2)-a2*x(1)+v(3);
x(4)=-a1*x(3)-a2*x(2)-a3*x(1)+v(4);
for n=1:N-4;
    x(n+4)=-a1*x(n+3)-a2*x(n+2)-a3*x(n+1)-a4*x(n)+v(n+4);   %产生真实数据
end
v2=randn(1,N);%观测噪声
z=x+v2;%产生观测数据
%维纳滤波
for m=1:N
   r(m)=mean(z(m)*z(m));
   p(m)=mean(x(m)*z(m));
   h(m)=p(m)/r(m);
   y(m)=h(m)*z(m);%就是16到20行这一段换上正确的维纳滤波的程序就可以了,当然可能别人还设了其他的变量,你看看怎么改,前面的部分就是产生真实数据和观测数据 
end
%R=xcorr(z);
%估计值与真实值之差
%error_px=X(1,:)-x;%x距离误差
%绘图
t=1:N;
figure(1);
plot(t,x,'k-',t,z,'g:',t,y,'r-.');
legend('真实轨迹','观测样本','估计轨迹');