求维纳滤波Matlab代码

%实例一程序_维纳滤波的计算机实现
%初步处理，并接受输入数据，包括信号样本个数L和滤波器阶数N
clear all
close all
L=input('L=');
N=input('N=');
a=0.95;
%定义w，v，u
w=sqrt(3*(1-a^2))*(2*rand(1,L)-1);
v=sqrt(3)*(2*rand(1,L)-1);
u=ones(1,L);
%获得原信号s(n)和带噪声信号x(n)
s(1)=1;
for i=2:L
s(i)=a*s(i-1)+w(i);
end
for i=1:L
x(i)=s(i)+v(i);
end
%绘图比较原信号s(n)和带噪声信号x(n)
figure
k=(L-99):L;
plot(k,s(k),'r',k,x(k),'b');
legend('s(n)','x(n)',0);
title('comparation between s(n) and x(n)');
xlabel('n');ylabel('Input');
%计算信号x(n)的N阶自相关矩阵Rxx，x(n)与s(n)的互相关函数向量rxs，并获得估计FIR滤波器系数h1
phixx=xcorr(x,x);
for i=1:N
for j=1:N
Rxx(i,j)=phixx(i-j+L);
end
end
phixs=xcorr(x,s);
for i=1:N
rxs(i)=phixs(i+L);
end
h1=(inv(Rxx))*rxs';
%获得理想FIR滤波器系数h1
for i=1:N
h(i)=0.238*0.724^i*u(i);
end
%实例一程序
%绘图比较估计滤波器与理想滤波器
figure
k=1:N;
plot(k,h(k),'r',k,h1(k),'b');
title('Ideal h(n) & Calculated h(n)');
legend('Ideal h(n)',' Calculated h(n)');
xlabel('n');ylabel('h(n)');
%计算并绘图比较理想输出与实际输出
S=conv(h,v);
SI(1)=S(1);
for i=2:L
SI(i)=0.724*SI(i-1)+0.238*x(i);
end
figure
k=(L-99):L;
plot(k,s(k),'r',k,SI(k),'b');
title('s(n) VS. SI(n)');
legend('s(n)','SI(n)',0);
xlabel('n');ylabel('Ideal Output');

SR=conv(h1,x);
figure
k=(L-99):L;
plot(k,s(k),'r',k,SR(k),'b');
title('s(n)VS. SR(n)');
legend('s(n)','SR(n)',0);
xlabel('n');ylabel('Actual Output');

EX2=0;
EI2=0;
ER2=0;
%计算并输出所获得信号与原信号的均方误差，理想维纳滤波和估计维纳滤波均方误差
for i=1:L;
EX2=1/L*((x(i)-s(i))^2)+EX2;
end;
for i=1:L;
EI2=1/L*((SI(i)-s(i))^2)+EI2;
end;
for i=1:L;
ER2=1/L*((SR(i)-s(i))^2)+ER2;
end;

EX2
EI2
ER2