ardj.m

AR模型定阶

%AIC定阶，书上193页信号
clear all
tic,
M=1000;
sg=zeros(1,10);
Ag=zeros(10,10);
rr=zeros(M,10);
for Q=1:M
    
    N=500;
    a=[1.16 -0.37 -0.11 0.18];
    wn=randn(1,N);
    x=zeros(1,N);
    x(2)=wn(2);
    x(3)=a(1)*x(2)+a(2)*x(1)+wn(3);
    x(4)=a(1)*x(3)+a(2)*x(2)+a(3)*x(1)+wn(4);
   for j=5:N
    
       x(j)=+a(1)*x(j-1)+a(2)*x(j-2)+a(3)*x(j-3)+a(4)*x(j-4)+wn(j);
  
   end
   xx=mean(x);
   gama=zeros(1,N);
   for k=0:N-1
       for j=1:N-k
        gama(k+1)=gama(k+1)+1/N*(x(j)-xx)*(x(j+k)-xx);
       end
  end
%对阶数计算协方差函数和自协方差矩阵，并利用yule walker方程计算相应的参数
 ag=zeros(10,10);
   for p=1:10
       G=zeros(p,p);    
       for ss=1:p
           for tt=1:p
               G(ss,tt)=gama(abs(ss-tt)+1);
           end
       end
       rr(Q,1:p)=gama(2:p+1);  
       ag(p,1:p)=(inv(G)*rr(Q,1:p)')';
       sigma(p)=gama(1)-dot(ag(p,1:p),rr(Q,1:p));
       AIC(p)=log(sigma(p))+2*p/N;
   end
   Ag=Ag+ag;
   sg=sg+sigma;    
   pp(Q)=find(AIC==min(AIC));
end
Agj=Ag/M;
sgj=sg/M;   
for h=1:10
   js(h)=length(find(pp==h));
end
js
toc,