AR.rar_walker_yuler-walker

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_) %求平均值