burg1.m

信号的谱估计的相关问题例子

close all;
clear;
clc;
fs=10000;                                       %采样频率设置（Hz）
f1=fs*0.2;
f2=fs*0.25;
N=100;                                      %信号的数据长度
t=0:1/fs:(N-1)/fs;
randn('state',0);
wn=randn(1,N);                                %产生高斯白噪声信号
SNR=10;                                       %信噪比
Pwn=sum(abs(wn).^2)/N;                        %噪声信号的功率
A1=sqrt(2*Pwn*10^(SNR/10));                   %由信噪比计算正弦信号的振幅
A2=A1;
sn1=A1*sin(2*pi*f1*t);                        %正弦序列
sn2=A2*sin(2*pi*f2*t);
sn=sn1+sn2;                                   %有用信号
Psn1=sum(abs(sn1).^2)/N;                      %正弦信号的功率
Psn2=sum(abs(sn2).^2)/N;
SNRreal1=10*log10(Psn1/Pwn);                  %计算实际信躁比
SNRreal2=10*log10(Psn2/Pwn);
xn=sn+wn;                                     %输入信号＝有用信号＋高斯白噪声

ef(1,:)=xn;
eb(1,:)=xn;

rou(1)=sum(abs(xn).^2)/N;

%p=1
add1=0;
add2=0;
for n=2:N,
    add1=add1+ef(1,n)*eb(1,n-1);
    add2=add2+abs(ef(1,n))^2+abs(eb(1,n-1))^2;
end
a(1,1)=-2*add1/add2;
rou(2)=(1-abs(a(1,1))^2)*rou(1);
sigma(1)=rou(2);
for n=3:N,
    ef(2,n)=ef(1,n)+a(1,1)*eb(1,n-1);
    eb(2,n)=eb(1,n-1)+a(1,1)*ef(1,n);
end

% %p=2
% add1=0;
% add2=0;
% for n=3:N,
%     add1=add1+ef(2,n)*eb(2,n-1);
%     add2=add2+abs(ef(2,n))^2+abs(eb(2,n-1))^2;
% end
% a(2,2)=-2*add1/add2;
% rou(3)=(1-abs(a(2,2))^2)*rou(2);
% sigma(2)=rou(3);
% a(2,1)=a(1,1)+a(2,2)*a(1,1);
% for n=4:N,
%     ef(3,n)=ef(2,n)+a(2,2)*eb(2,n-1);
%     eb(3,n)=eb(2,n-1)+a(2,2)*ef(2,n);
% end

%迭代
IP=20;
for p=2:IP,
    add1=0;
    add2=0;
    for n=p+1:N,
        add1=add1+ef(p,n)*eb(p,n-1);
        add2=add2+abs(ef(p,n))^2+abs(eb(p,n-1))^2;
    end
    a(p,p)=-2*add1/add2;
    for n=p+2:N,
        ef(p+1,n)=ef(p,n)+a(p,p)*eb(p,n-1);
        eb(p+1,n)=eb(p,n-1)+a(p,p)*ef(p,n);
    end    
    rou(p+1)=(1-abs(a(p,p))^2)*rou(p);
    sigma(p)=rou(p+1);
    for k=1:p-1,
        a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);
    end
end

nfft=1024;
if mod(nfft,2)==0
    Plength=nfft/2+1;
elseif mod(nfft,2)==1
    Plength=(nfft+1)/2;
end

F=linspace(0,fs/2,Plength);                  %频率坐标
W=F*2*pi/fs;                                 %数字角频率
for omiga=1:Plength
    add3(omiga)=0;
    for l=1:p
        add3(omiga)=add3(omiga)+a(p,l)*exp(-j*W(omiga)*l);
    end
    Pxx(omiga)=sigma(p)/abs(1+add3(omiga))^2;
end

subplot(211);
plot(F,10*log10(Pxx));
grid on;
xlabel('频率(Hz)');
ylabel('功率谱幅度(dB)');
title('编程实现Burg法功率谱估计');

[PxxB,FB]=pburg(xn,p,nfft,fs);
subplot(212);
plot(FB,10*log10(PxxB));
grid on;
xlabel('频率(Hz)');
ylabel('功率谱幅度(dB)');
title('命令实现Burg法功率谱估计');