AR模型估计功率谱

假设信号为，添加的噪声为方差0.1的白噪声。请采用AR模型估计其功率谱。

• 设计过程

1. 分析时间序列和回归一样，目的都是预测。在回归里面，我们有一元回归于多元回归，在时间序列里面，我们有自回归。与一元、多元一样，我们分为一阶与多阶自回归。其实还是那样的理念，只不过之前是变量与应变量，现在是存在时滞的序列之间的关系。

一阶自回归AR，也就是，Yt=b*Yt-1+ut。

1. 首先来生成一个时间序列，其自回归方程如下：

yt = 0.8 * yt-1 + c

其中c是残差项，我们用白噪音，也就是正态分布来表示。

1. AR模型阶数的选取：采用AIC准则进行选取。

 

选取最小的K，就是AR模型的阶次。

 

   三、程序

N=400;

%%%%%%%  生成x(n)的前200个样本  %%%%%%%%%%%%%%%%%%%%%%%%

R = normrnd(0,0.1,N);

for i=1:N

  x(i)=R(i)+sin(i/100)*sin(i/20);

end

 

theta=[];

for i=1:N

    sum=0;

    for k=1:N-i

        sum=sum+x(k)*x(i+k-1);

    end

    theta(i)=1/(N-i+1)*sum;

end

 

 

for k=-N:N

    w=pi*k/N;

    p=0;

    for i=2:N

        p=p+theta(i)*(exp(-j*w*(i-1))+exp(j*w*(i-1)));

    end

    q(k+N+1)=p+theta(1);

end

 

for k=-N:N

    w=pi*k/N;

     X(k+N+1)=0;

    for i=1:N

        X(k+N+1)= X(k+N+1)+x(i)*exp(-j*w*(i-1));

    end

end

for i=1:2*N+1

    X(i)=(1/N)*(abs(X(i)))^2;

end    

 

 

%%%%%%%%%%%%%%%%ARMA模型%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=100;

for l=1:s

   W=zeros(l+1,l+1);

   NI=zeros(l+1,l+1);

   %%%%%%%        构造矩阵M     %%%%%%%%%%%%%

   for e=1:l+1

      for f=e:l+1

        W(e,f)=theta(f-e+1);

        W(f,e)=W(e,f);

      end

   end

   %%%%%%  利用AIC准则判断最优阶次  %%%%%%%%%%

  NI=W^(-1);

 det=1/NI(1,1);

 AIC(l)=N*log(det)+2*l;

 

end

det1=min(AIC);

%%%%%%%%  找到最优阶次  %%%%%%%%%%%%%%%%%%%%%%

for i=1:s

    if AIC(i)==det1;

        disp('最优阶次是：')

        jc=i

    end

end

%%%%%%%%%%%%%%  利用AR来计算功率谱密度  %%%%%%

 bb=NI*det1;

for i=1:s

    aa(i,1)=bb(i+1,1);

end

for k=-N:N

     w=pi*k/N;

     SUM=0;

     for i=1:s

         SUM=SUM+aa(i,1)*exp(-j*w*i);

     end

     u(k+N+1)=det*(abs(1/(1+SUM)))^2;

end

plot(abs(q),'r');

hold on

plot(abs(X),'k')

hold on

plot(u)

legend('自相关法','周期图法','AR模型')

四、结果图

 

最优阶次是：jc =80