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