wiener_pz_var_1.m

维纳滤波的一个应用 是关于AR(4)模型的

clc;
clear all;
maxlag=100;
%%%%%%%%%%%%%%%%%%%%%%%%
N=100;   
x=zeros(N,1);
y=zeros(N,1);
var=1;
%%%%%%%%%%%%%%%%%%%%%%%%%列出状态方程
x(1)=randn(1,1);   %令x(-1)=x(-2)=x(-3)=x(-4)=0
x(2)=randn(1,1)+1.352*x(1);
x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);
x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);
for n=5:N   
    x(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1); %x为真实值
end;
v=randn(N,1);
y=x+v;  %z_x为观测样本值=真值+噪声
%%%%%%%%%%%%%%%%%%%%%%%%滤波
x = x';
y = y';
xk_s(1)=y(1);    %赋初值
xk_s(2)=y(2);
xk_s(3)=y(3);
xk_s(4)=y(4);
xk=[y(1);y(2);y(3);y(4)];  
%%%%%%%%%%%%%%%%%%%%%%%%
%维纳滤波器的生成
[rx,lags]=xcorr(y,maxlag,'biased');%观测信号的自相关函数
rx1=toeplitz(rx(101:end));%对称化自相关函数矩阵使之成为方阵，滤波器的阶数为101阶
rx2=xcorr(x,y,maxlag,'biased');%观测信号与期望信号的互相关函数
rx2=rx2(101:end);
h=inv(rx1)*rx2';%维纳-霍夫方程
xk_s=filter(h,1,y);%加噪信号通过滤波器后的输出
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
e_x=0;
eq_x=0;
e_x1=N:1;
%计算滤波的均值,计算滤波误差的均值
for i=1:N
    e_x(i)=x(i)-xk_s(i); %误差=真实值-滤波估计值
end
%%%%%%%%%%%%%%%%%%%%%%%作图
t=1:N;
figure(1);
plot(t,x,'r-',t,y,'g:',t,xk_s,'b-.');
legend('真实轨迹','观测样本','估计轨迹'); 
figure(2);
plot(e_x);          
legend('平均误差');