burg.m

用burg算法实现功率谱的估计

clc;
clear all
fs=1000; %设采样频率为100，
N=1000; %数据长度
P=50;%设置阶数
n=0:1/fs:(N-1)/fs; 

f0=100;f1=150%设置信号频率
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;%产生信号，并在信号中加入噪声

%(一)初始化:求出ef0,eb0和Pmin0（Pmin0=R0)
r=0;%计算自相关R0
for ni=1:N
    r=r+abs(x(ni))*abs(x(ni));
end
R0=r/N;
ef0=x;
eb0=x;

%(二)计算出一阶时(M=1)的反射系数K1，Pmin1及ef1,eb1
M=1;
fm=0;
fz=0;
for ni=M+1:N
    fz=fz+ef0(ni)*eb0(ni-1);
    fm=fm+(ef0(ni)*ef0(ni)+eb0(ni-1)*eb0(ni-1));
end
K1=-2*fz/fm;%算出反射系数K1
a11=K1;
Pmin1=(1-abs(K1)*abs(K1))*R0;%算出预测误差功率Pmin1
ef1(1)=ef0(1);%算出ef1和eb1
eb1(1)=K1*ef0(1);
for ni=2:N
    ef1(ni)=ef0(ni)+K1*eb0(ni-1);
    eb1(ni)=eb0(ni-1)+K1*ef0(ni);
end
a1(1)=a11;
M=M+1;%此时M=2

%(三)计算出两阶以上（M>=2)的参数。判断M<规定的阶数P。(若“小”时，M=M+1，继续;若“大”时结束)
while(M<P+1)
    fm=0;
    fz=0;
    for ni=M+1:N
        fz=fz+ef1(ni)*eb1(ni-1);
        fm=fm+(ef1(ni)*ef1(ni)+eb1(ni-1)*eb1(ni-1));
    end
    K2=-2*fz/fm;%算出反射系数K2
    
    for ni=1:M-1%计算M阶时的参数a(M,M)和a(M,K)(k=1,2,...,M-1)
        a2(ni)=a1(ni)+K2*a1(M-ni);
    end
    a2(M)=K2;
    
   Pmin2=(1-abs(K2)*abs(K2))*Pmin1;%算出预测误差功率Pmin2
   
    for ni=M+1:N%算出ef2和eb2
        ef2(ni)=ef1(ni)+K2*eb1(ni-1);
        eb2(ni)=eb1(ni-1)+K2*ef1(ni);
    end
    
    ef1=ef2;
    eb1=eb2;
    Pmin1=Pmin2;
    a1=a2;
    M=M+1;
    
end%结束while循环

%(四)参数确定后计算出功率谱
w=0:0.01:pi;
    s=0;
    for nj=1:P
        s=s+a2(nj)*exp(-j*w*nj);
    end
PSD=Pmin2./abs(s+1).^2;
F=w*fs/(pi*2); %角频率坐标换算成频率坐标
%plot(F,abs(PSD));
plot(F,10*log10(abs(PSD)));
xlabel('频率(Hz)'),ylabel('功率谱密度(dB)');
title('基于Burg算法的功率谱估计');
grid;