clc;
clear all;
M=2;
p=2;
L=100;
mu=0.01;
sigma_2=0.01;
a1(1)=0;
a2(1)=0;
a1_i=-1.3;
a2_i=0.8;
for k=1:500
wn = sqrt(sigma_2)*( randn(L,1));
A_i=[1 a1_i a2_i];
yn_i=filter(1,A_i,wn);
yn_ii=[zeros(L,1)
yn_i];
for i=1:L
yn(i+L)=-a1(i)*yn_ii(i-1+L)-a2(i)*yn_ii(i-2+L);
e(i)=yn_ii(i+L)-yn(i+L);
a1(i+1)=a1(i)-2*mu*e(i)*yn_ii(i-1+L); %得到a1的AR估计
a2(i+1)=a2(i)-2*mu*e(i)*yn_ii(i-2+L); %得到a2的AR估计
end
for i=1:L
a1_ii(i)=a1_i;
a2_ii(i)=a2_i;
end
figure(1),clf
subplot(2,1,1)
plot(a1_ii,'*r')
hold on
plot(a1,'ob')
legend('a1的理想值','a1的估计值');
grid on; axis tight;
subplot(2,1,2)
plot(a2_ii,'*r')
hold on
plot(a2,'ob')
legend('a2的理想值','a2的估计值');
grid on; axis tight;
r_yy=xcorr(yn_i,'unbiased');
Ryy=zeros(M+1,M+1);
ep=zeros(M+1,1);
ep(1)=1;
wn=zeros(M+1,1);
for i=1:p+1
Ryy(i,:)=r_yy(L-i+1:L-i+M+1);
end
Ryy_=inv(Ryy);
sigma2=1/Ryy_(1,1);
awn=Ryy_*(sigma2*ep);
a1wn_(k)=awn(2); %得到a1的Y-W方程估计
a2wn_(k)=awn(3); %得到a2的Y-W方程估计
a1_(k)=a1(L);
a2_(k)=a2(L);
end
figure(2),clf
subplot(2,1,1)
plot(a1_,'*-r')
hold on
plot(a1wn_,'-ob')
legend('a1的AR估计值','a1的Yule-Walker估计值');
grid on; axis tight;
subplot(2,1,2)
plot(a2_,'*-r')
hold on
plot(a2wn_,'-ob')
legend('a2的AR估计值','a2的Yule-Walker估计值');
grid on; axis tight;
a1a=mean(a1_)
a2a=mean(a2_)
a1awn=mean(a1wn_)
a2awn=mean(a2wn_) %求平均值
