zhouqitu.asv

用Matlab编写的功率谱估计

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);
R=ZXG(x);
%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,1024));
for i=1:1024
    P(i)=X(i)*X(i)/N;
end
plot(P);
%改进之后的周期图法,u为汉明窗
for i=1:99
    u(i)=R(i)*han(i);
end
UU=abs(fft(u,1024));
figure(2);
plot(UU);
axis([0 100 0 1000000]);
grid on

%三阶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,256));
for i=1:256
    PP(i)=D_3/((1+D(i))*(1+D(i)));
end
figure(3);
plot(PP);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%因为AR谱估计高阶算法与三阶算法类似， %
%%%而这种方法与线性预测谱估计等效，故                  %
%%%可以直接直接调用求线性预测的系统函数LPC(X,N)
L=LPC(x,70);
LL=abs(fft(L,256));
for i=1:256
    PPP(i)=D_3/((1+LL(i))*(1+LL(i)));
end
figure(4);
plot(PPP);
%axis([0 128 0 5]);