www.pudn.com > 2wiener.rar > 2wiener.m

clc; 
clear all; 
maxlag=100; 
 
N=100;   %观测点数取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('平均误差');