L = input('输入信号样本个数L=');
N = input('输入FIR滤波器阶数N=');
a = 0.95;
K = 50;
sigma_a2 = 1-a^2;
a_ = [1,-a];
while(1)
w = sqrt(sigma_a2)*(randn(L,1)); %获得方差为sigma_a2的高斯白噪声w
v = randn(L,1); %获得方差为1的高斯白噪声v
s(1) = w(1); %通过公式：s(n) = a*s(n-1)+wn(n)获得源信号s与加了噪声v的x
for i=1:L-1
s(i+1)=a*s(i)+w(i+1);
x(i) = s(i)+v(i);
end
x(L) = s(L)+v(L);
%%%%%%%%%%%%%%%%%估计相关函数r_xx和r_xs%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:K+1
rxx(i)=sum(x(i:L).*x(1:L-i+1))/(L-i+1);
end
r_xx_g=[rxx(K+1:-1:2),rxx(1:K+1)];
for i=1:K+1
rxs(i)=sum(x(i:L).*s(1:L-i+1))/(L-i+1);
end
r_xs_g=[rxs(K+1:-1:2),rxs(1:K+1)];
%%%%%%%%%%%%%%检验x的r_xx和r_xs是否与理论值相符%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r_xx_t = a.^abs([-K:K]);
r_xx_t(K+1) = r_xx_t(K+1)+1;
r_xs_t = a.^abs([-K:K]);
rou_xx=(sum((r_xx_g-r_xx_t).^2))/sum(r_xx_t.^2);
rou_xs=(sum(r_xs_g-r_xs_t).^2)/sum(r_xs_t.^2);
if rou_xx<0.03 & rou_xs<0.01
break;
end
end
%%%%%%%%%%%%%%x(n)的自相关函数的理论值(红色)和实际值(蓝色)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
stem(r_xx_t,'r');
hold on
stem(r_xx_g,'*','b');
title('x(n)的自相关函数的理论值和实际值(*)');
%%%%%%%%%%%%%%最后100个s(n)(红色)和x(n)(蓝色)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);
stem(s(L-100:L),'r');
hold on
stem(x(L-100:L),'*','b');
title('最后100个s(n)和x(n)(*)');
%%%%%%%%%%%%%%N个理想h(n)(h_t红色)估计h(n)(h_g蓝色)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = 0:N-1;
h_t = 0.238*(0.724).^n; %求得理想的h(n)
for i=1:N
xx(i)=sum(x(i:L).*x(1:L-i+1))/(L-i+1);
end
for i=1:N
Rxx(i,1:N)=[xx(i:-1:1),xx(2:N+1-i)];
end
for i=1:N
xs(i)=sum(x(i:L).*s(1:L-i+1))/(L-i+1);
end
invRxx=inv(Rxx);
h_g=invRxx*xs'; %估计h
figure(3);
stem(h_t,'r');
hold on
stem(h_g,'*','b');
title('N个理想h(n)(h_t)估计h(n)(*h_g)');
%%%%%%%%%%%%%%最后100个s(n)(红色)理想维纳滤波S_l(n)(蓝色)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_l(1) = s(1); %得到理想维纳滤波的S(n)
for i = 1:L-1
y_l(i+1) = 0.724*s(i)+0.238*x(i+1);
end
figure(4);
stem(s(L-100:L),'r');
hold on
stem(y_l(L-100:L),'*','b');
title('最后100个s(n)理想维纳滤波S_R(n)(*)');
%%%%%%%%%%%%%%最后100个s(n)(红色)根据式FIR滤波S_R(n)(蓝色)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:9 %得到FIR维纳滤波的S(n)
y_F(i)=sum(h_g(1:i)'.*x(i:-1:1));
end
for i=10:L
y_F(i)=sum(h_g'.*x(i:-1:i-9));
end
figure(5);
stem(s(L-100:L),'r');
hold on
stem(y_F(L-100:L),'*','b');
title('最后100个s(n)根据式FIR滤波S_R(n)(*)');
%%%%%%%%%%%%%%%%%求均方误差%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
e_x=sum((x-s(1:L)).^2)/L %求各自的均方误差
e_i=sum((y_l-s(1:L)).^2)/L
e_f=sum((y_F-s(1:L)).^2)/L
如果觉得《matlab实习维纳滤波 自编Matlab程序 维纳滤波器的Matlab实现》对你有帮助，请点赞、收藏，并留下你的观点哦！
