clear
clc
% NT=100;
% ex=zeros(NT,1);
% ei=zeros(NT,1);
% er=zeros(NT,1);
% for m=1:NT
L = 10000;
a = 0.95;
K = 50;
sigma_a2 = 1-a^2;
a_ = [ 1, -a];
while(1)
wn = sqrt(sigma_a2)*( randn(L,1));
sn = filter(1, a_, wn);
vn = randn(L,1);
xn = sn + vn;
r_xx = xcorr(xn,'unbiased');
r_xx_t = a.^abs([-K:K]);
r_xx_t(K+1)=r_xx_t(K+1)+1;
p = xcorr(sn,xn,'unbiased');
r_xs = p(L : L+K);
rou_xx = sum((r_xx(L-K:L+K)-r_xx_t').^2)/sum(r_xx_t.^2);
rou_xs = (sum(r_xs-a.^[0:K]').^2)/sum(a.^[0:K].^2);
if rou_xx<0.03 & rou_xs<0.01
break;
end
end
% figure(1),clf
% stem(r_xx(L-K:L+K),'*r')
% hold on
% stem(r_xx_t,'.b')
% title('xn的理论和实际自相关函数值');
% xlabel('t');
% ylabel('自相关函数值');
% legend('xn自相关函数实际值','xn自相关函数理论值');
% figure(2),clf
% x=[L-99:L];
% plot(x,sn(x),'--.r',x,xn(x),'--*b')
% title('sn和xn最后100个数值');
% xlabel('t');
% ylabel('数值');
% legend('sn最后100个值','xn最后100个值')
N=10;
R_xx=zeros(N,N);
R_xs=zeros(N,1);
hI=zeros(N,1);
for i=1:N
R_xx(i,:)=r_xx(L-i+1:L-i+N);
R_xs(i)=r_xs(i);
hI(i)=0.238*0.724^(i-1);
end
hF=inv(R_xx)*R_xs;
figure(3),clf
stem(hI,'.r')
hold on
stem(hF,'*b')
title('L=10000,N=10,理想的h(n)及其估值');
xlabel('t');
ylabel('函数值');
legend('理想的h(n)','估计的h(n)');
sI=zeros(L,1);
sI(1)=0.238*xn(1);
for i=1:L-1
sI(i+1)=0.724*sI(i)+0.238*xn(i+1);
end
% figure(4),clf
% x=[L-99:L];
% plot(x,sn(x),'--.r',x,sI(x),'--*b')
% title('sn和理想滤波sI最后100个数值');
% xlabel('t');
% ylabel('函数值');
% legend('sn','sI');
sR=zeros(L,1);
sR=conv(hF,xn);
% figure(5),clf
% x=[L-99:L];
% plot(x,sn(x),'--.r',x,sR(x),'--*b')
% title('sn和FIR滤波后的sR最后100个数值');
% xlabel('t');
% ylabel('函数值');
% legend('sn','sR');
% ex2=0;
% ei2=0;
% er2=0;
% ex2=sum((xn(1:L)-sn(1:L)).^2)/L;
% ei2=sum((sI(1:L)-sn(1:L)).^2)/L;
% er2=sum((sR(1:L)-sn(1:L)).^2)/L;
% ex(m)=ex2;ei(m)=ei2;er(m)=er2;
% end
% eex2=sum(ex)/NT
% eei2=sum(ei)/NT
% eer2=sum(er)/NT