levinson.m

AR模型参数确定后计算出功率谱, 用Levision_Durbin 算法

clc;
clear all
fs=100; %设采样频率为100，
N=150 %数据长度,改变数据长度会导致分辨率的变化；
n=0:1/fs:(N-1)/fs; 

f0=20;f1=25;%设置信号频率
SNR=10;%设置信噪比
A=sqrt(10.^(SNR/10));%信号幅度值
vx=randn(1,length(n));%噪声
x=A*cos(2*pi*f0*n)+A*sin(2*pi*f1*n)+vx;%产生信号，并在信号中添加噪声

%计算自相关函数
for m=1:N
    R(m)=0;
end
for m=1:N
    s=0;
    for n=1:N-(m-1)
        s=s+x(n)*x(n+m-1);
    end
    R(m)=s/N;
end 

%迭代
for M=1:N-5%阶数M要小于数据长度N，这里取M=N-5
    a(M,M)=0; 
    FPE(M)=0;
   P(M)=0;
end
%（一)初始化：由Levinson公式计算一阶模型参数a(1，1)和误差功率P（1）及FPE（1)；
a(1,1)=-R(2)/R(1); %计算一阶模型参数a(1,1)；
P(1)=R(1)+a(1,1)*R(2); %计算一阶模型的误差功率P(1)；
FPE(1)=P(1)*(N+1+1)/(N-1-1); %计算一阶模型的最终预测误差准则函数FPE(1)；

%（二）递推：设已经计算出M-1阶时a(M-1,k)和P(M-1),计算M阶模型参数a(M,M),a(M,k)和误差功率P（M）及FPE（M)；
 sum=0;
for M=2:N-5 
    for k=1:M-1
        sum=sum+a(M-1,k)*R(M-k+1);
    end
    a(M,M)=-(R(M+1)+sum)/P(M-1);%计算反射系数
    for k=1:M-1
        a(M,k)=a(M-1,k)+a(M,M)*a(M-1,M-k);
    end
   P(M)=(1-(abs(a(k,k)))^2)*P(M-1);
    FPE(M)=P(M)*(N+M+1)/(N-M-1);
    sum=0;
end %递推计算完毕；

%(三)Model阶次确定：求出使FPE最小时的阶数M1；
min=FPE(1); 
for M=2:N-5
    if FPE(M)<min
        min=FPE(M);
        M1=M;
    end
end 
disp('阶数M为：')
disp(M1);
disp('模型参数a为：');
for i=1:M1
    disp(a(M1,i));
end
disp('误差功率为：');
disp(P(M1));

%(四)AR模型参数确定后计算出功率谱
Z=0;
W=0:0.01:pi; %功率谱以2*pi为周期，又信号为实信号，只需输出0到PI即可；
for k=1:M1
    Z=Z+(a(M1,k).*exp(-j*k*W));
end
PSD=P(M1)./((abs(1+Z)).^2); %算出功率谱

F=W*fs/(pi*2); %角频率坐标换算成频率坐标
plot(F,abs(PSD));
xlabel('Frequency(Hz)'),ylabel('Power Spectral Density');
title('基于Levison Durbin方法的功率谱估计');
grid