ar.m

AR模型谱估计算法

  clear;
% 基本参数
p=10;
pi=3.1415;
f1=0.05;
f2=0.4;
f3=0.42;
a1=-0.850848;
a=0; 
b=sqrt(0.101043/2); 
% 噪声
u1=randn(1,32); 
u1=u1/std(u1); 
u1=u1-mean(u1); 
u1=a+b*u1; 

u2=randn(1,32); 
u2=u2/std(u2); 
u2=u2-mean(u2); 
u2=a+b*u2; 
u=u1+u2*i;

% 信号产生
z=zeros(1,32);
x=zeros(1,32);
z(1)=u(1);
for n=2:1:32,z(n)=-a1*z(n-1)+u(n);end,z;
for n=1:1:32,x(n)=2*cos(2*pi*f1*(n-1))+2*cos(2*pi*f2*(n-1))+2*cos(2*pi*f3*(n-1))+z(n);end,x;

% 求Cxx(j,k)矩阵
Cxx1=zeros(p,p);
Cxx2=zeros(p,p);
for j=1:1:p
    for k=1:1:p
        for n=p+1:1:32
            Cxx1(j,k)=Cxx1(j,k)+conj(x(n-j))*x(n-k);
        end
    end
end
for j=1:1:p
    for k=1:1:p
        for n=1:1:(32-p)
            Cxx2(j,k)=Cxx2(j,k)+x(n+j)*conj(x(n+k));
        end
    end
end
Cxx=Cxx1+Cxx2;
Cxx=Cxx/(2*(32-p));

% 求Cxx(k,0)
Cxx3=zeros(p,1);
Cxx4=zeros(p,1);
for j=1:p
    for n=(p+1):32
        Cxx3(j,1)=Cxx3(j,1)+conj(x(n-j))*x(n);
    end
end
for j=1:p
    for n=1:(32-p)
        Cxx4(j,1)=Cxx4(j,1)+x(n+j)*conj(x(n));
    end
end
Cxx5=Cxx3+Cxx4;
Cxx5=Cxx5/(2*(32-p));

% 求a(k)
A=inv(Cxx)*(-Cxx5);

% 求Cxx(0,0)
Cxx6=0;
for  n=(p+1):32
    Cxx6=Cxx6+conj(x(n))*x(n);
end
Cxx7=0;
for n=1:(32-p)
     Cxx7=Cxx7+x(n)*conj(x(n));
 end
 Cxx8=(Cxx6+Cxx7)/(2*(32-p));
 
%  求Cxx(0,k)
Cxx9=zeros(1,p);
Cxx10=zeros(1,p);
for k=1:p
    for n=(p+1):32
        Cxx9(k)=Cxx9(k)+conj(x(n))*x(n-k);
    end
end
for k=1:p
    for n=1:(32-p)
         Cxx10(k)=Cxx10(k)+x(n)*conj(x(n+k));
    end
end
Cxx11=(Cxx9+Cxx10)/(2*(32-p));

% 求噪声方差
Q=0;
for k=1:p
    Q=Q+A(k)*Cxx11(k);
end
Q=Q+Cxx8;
Q=abs(Q);

% 求功率谱密度、画图
 y=0;
 f=-0.5:0.001:0.5;
 for k=1:p
       y=y+A(k)*exp(-2*pi*f*k*i);
 end
 y=y+1;
 y=abs(y);
 y=y.^2;
 Paf=Q./y;
 Paf=10*log10(Paf);
 plot(f,Paf);