www.gusucode.com > AR模型功率谱估计源码程序 > AR模型功率谱估计源码程序/功率谱估计/Mainfun.m

    clear;
clc;
x = [101 82 66 35 31 7 20 92 154 125 85 68 38 23 10 24 83 132 131 118 90 67 60 47 41 21 16 6 4 7 14 34 45 43 48 42 28 10 8 2 0 1 5 12 14 35 46 41 30 24 16 7 4 2 8 17 36 50 62 67 71 48 28 8 13 57 122 138 103 86 63 37 24 11 15 40 62 98 124 96 66 64 54 39 21 7 4 23 55 94 96 77 59 44 47 30 16 7 37 74];
N=length(x);
%调用自相关函数求得x的自相关R
R=ZXG(x);
H=1024;%FFT点数
%han=zeros(1,1024);
%定义一个窗长为N的汉明窗
for i=1:N
    han(i)=100*(0.54-0.46*cos(2*pi*i/N));
end

%用周期图法求数据的功率谱估计P
X=abs(fft(x,H));
for i=1:H
    P(i)=X(i)*X(i)/N;
end

%改进之后的周期图法,u为汉明窗
for i=1:99
    u(i)=R(i)*han(i);
end
UU=abs(fft(u,H));

%三阶AR模型参数得其功率谱
a_11=-R(2)/R(1);
D_1=(1-abs(a_11)*abs(a_11))*R(1);
a_22=-(R(3)+a_11*R(2))/D_1;
a_21=a_11+a_22*a_11;
D_2=(1-abs(a_22)*abs(a_22))*D_1;
a_33=-(R(4)+a_21*R(3)+a_22*R(2))/D_2;
a_32=a_22+a_33*a_21;
a_31=a_21+a_33*a_22;
D_3=(1-abs(a_33)*abs(a_33))*D_2;

A=[a_31 a_32 a_33];
D=abs(fft(A,H));
for i=1:H
    PP(i)=D_3/((1+D(i))*(1+D(i)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  因为AR谱估计高阶算法与三阶算法类似，而这种  %%%
%%%  方法与线性预测谱估计等效，故可以直接直接调  %%%
%%%  用求线性预测的系统函数LPC(X,N)求取高阶AR   %%%
%%%  谱估计                                   %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%取阶数为70
L=LPC(x,70);
LL=abs(fft(L,H));
for i=1:H
    PPP(i)=D_3/((1+LL(i))*(1+LL(i)));
end

%axis([0 128 0 5]);
%Burg法求解AR模型的参数，得到的功率谱
J=30;
%调用Burg算法
[K,D]=burg_Kp(x,J);
%J表示所取阶数

A=zeros(1,J+1);
A(1)=1;
A(2:J+1)=K(J,:);
D2=D(J+1);
FA=abs(fft(A,H));
for i=1:H
    PA(i)=D2/((1+FA(i))*(1+FA(i)));
end

%画周期图法谱估计
figure(1);
subplot(2,1,1);
plot(P);
axis([0 100 0 100000]);
title('周期图法谱估计(FFT点数:1024)');
xlabel('n');
ylabel('P(w)');
%axis([0 100 -15 15]);
grid on;

%画改进后的周期图法谱估计(加窗)：
subplot(2,1,2);
plot(UU);
axis([0 100 0 1000000]);
title('改进后的周期图法谱估计(FFT点数:1024)');
xlabel('n');
ylabel('P(w)');
grid on;
%自回归AR模型谱估计(Levinson-Durbin)法
%画三阶AR模型谱估计
figure(2);
subplot(2,1,1);
plot(PP);
title('三阶AR模型谱估计Levinson-Durbin法');
xlabel('n');
ylabel('P(w)');
grid on;

%画高阶(70阶)AR模型谱估计
subplot(2,1,2);
plot(PPP);
axis([0 100 0 400]);
title('高阶(70阶)AR模型谱估计Levinson-Durbin法');
xlabel('n');
ylabel('P(w)');
grid on;

%画Burg法求解AR模型的参数，得到的功率谱曲线
%阶数选30
figure(3);
subplot(2,1,1);
plot(PA);
axis([0 128 0 800]);
title('30阶AR模型谱估计Burg法(FFT点数:1024)');
xlabel('n');
ylabel('P(w)');
grid on;