p_morden.m

用Matlab实现现代功率谱估计

clear;close all;clc;
L=1024;%采样点数为1024
SNR1=-10;SNR2=-10;%SNR1和SNR2分别为两路正弦信号的信噪比，单位为dB
A1=sqrt(2*10^(SNR1/10));%根据dB单位的信噪比，换算出两路信号的振幅
A2=sqrt(2*10^(SNR2/10));
f1=100;f2=110;fs=1000;%两路信号的频率
Num_fft=1024;%进行FFT变换时的点数
N=40;%谱估计的阶数
randn('state',0);
w=randn(1,L);%产生噪声序列
t=0:1/fs:(L-1)/fs;%根据采样率划分采样时间
x=A1*sin(2*pi*f1*t)+A2*sin(2*pi*f2*t)+w;%产生观测信号
%xlim([0,500]);

f_begin=0;
f_end=200;
figure(1)
[Pxx1,f]=pyulear(x,N,Num_fft,fs);%用Levinson法进行功率谱估计
Pxx1=Pxx1/min(Pxx1);%找出谱密度中的最小值，作归一化处理
Pxx1=10*log10(Pxx1);%产生谱密度的dB数据
subplot(2,2,1);
plot(f,Pxx1);grid on;xlim([f_begin,f_end]);
title('Levinson法估的功率谱密度曲线');
xlabel('频率f');
ylabel('功率谱密度 (dB)');

[Pxx2,f]=pburg(x,N,Num_fft,fs);  %用Burg法进行功率谱估计
Pxx2=Pxx2/min(Pxx2);%找出谱密度中的最小值，作归一化处理
Pxx2=10*log10(Pxx2);%产生谱密度的dB数据
subplot(2,2,2);
plot(f,Pxx2);grid on;xlim([f_begin,f_end]);
title('Bug法估的功率谱密度曲线');
xlabel('频率f');
ylabel('功率谱密度 (dB)');

%figure(2)
[Pxx3,f]=pcov(x,N,Num_fft,fs);%用协方差法进行功率谱估计
Pxx3=Pxx3/min(Pxx3);%找出谱密度中的最小值，作归一化处理
Pxx3=10*log10(Pxx3);%产生谱密度的dB数据
subplot(2,2,3);
plot(f,Pxx3);grid on;xlim([f_begin,f_end]);
title('协方差法估的功率谱密度曲线');
xlabel('频率f');
ylabel('功率谱密度 (dB)');

[Pxx4,f]=pmcov(x,N,Num_fft,fs);  %用修正协方差法进行功率谱估计
Pxx4=Pxx4/min(Pxx4);%找出谱密度中的最小值，作归一化处理
Pxx4=10*log10(Pxx4);%产生谱密度的dB数据
subplot(2,2,4);
plot(f,Pxx4);grid on;xlim([f_begin,f_end]);
title('修正协方差法估的功率谱密度曲线');
xlabel('频率f');
ylabel('功率谱密度 (dB)');