gusucode.com > 自相关法求AR模型谱估计源码程序 > 自相关法求AR模型谱估计源码程序/code/ARzxg.m

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   自相关法求AR模型谱估计 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
close all;
clear all;


N=32;%信号长度
f1=0.05;
f2=0.4;
f3=0.42;
A1=-0.850848;
p=10; %AR模型阶次



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     生成高斯白噪声序列 V(n)    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

V1=randn(1,N);
V2=randn(1,N);
U=0; %噪声均值
Q=0.101043;%噪声方差
b=sqrt(Q/2);
V1=U+b*V1;%生成1*N阶均值为U,方差为Q/2的高斯白噪声序列
V2=U+b*V2;%生成1*N阶均值为U,方差为Q/2的高斯白噪声序列

V=V1+j*V2;%生成1*N阶均值为U,方差为Q的复高斯白噪声序列


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        生成一阶自回归过程 z(n)序列       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

z(1)=V(1,1);
for n=2:1:N
    z(n)=-A1*z(n-1)+V(1,n);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        生成x(n)序列       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x(1)=6;
for n=2:1:N
    x(n)=2*cos(2*pi*f1*(n-1))+z(n-1);
end


%%%%%%%%%%%%%%%%%%%%%%%%
%    自相关法求Rxx(k)   %
%%%%%%%%%%%%%%%%%%%%%%%%

for k=0:1:p
    t5=0;
    for n=0:1:N-k-1
        t5=t5+conj(x(n+1))*x(n+1+k);
    end
    Rxx(k+1)=t5/N;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    由自相关函数Rxx(k) 求a(k,k)   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a(1,1)=-Rxx(2)/Rxx(1);
p1(1)=(1-abs(a(1,1))^2)*Rxx(1);

for k=2:1:p
    t=0;
    for l=1:1:k-1
        t=a(k-1,l).*Rxx(k-l+1)+t;
    end
    a(k,k)=-(Rxx(k+1)+t)./p1(k-1);
    for i=1:1:k-1
        a(k,i)=a(k-1,i)+a(k,k)*conj(a(k-1,k-i));
    end
    p1(k)=(1-(abs(a(k,k)))^2).*p1(k-1);
end


%%%%%%%%%%%%%%%%%%%%%%%%
%    由a(p,k)  求a(k)  %
%%%%%%%%%%%%%%%%%%%%%%%%

for k=1:1:p
    a(k)=a(p,k);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    生成x(n)的功率谱 X(f)   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f=-0.5:0.0001:0.5;
f0=length(f);
for t=1:f0
    s=0;
    for k=1:p
        s=s+a(k)*exp(-j*2*pi*f(t)*k);
    end
    X(t)=Q/(abs(1+s))^2;
end

plot(f,10*log10(X))
xlabel('频率');
ylabel('PSD(dB)');
title('自相关法求AR模型谱估计')
legend('p=10')