ly1.m

MATLAB实现功率谱估计的程序最大熵谱法,Brug法和Levinson法。

clear;                      
clc;
N=1024;                                                     %点数
Nfft=128;                                                   %定义变换点数
n=[0:N-1];
randn('state',0);
wn=randn(1,N);                                             %定义一均值为0,方差为r的白噪声wn
xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(20)*sin(2*pi*0.3*n)+wn;   %定义初始信号，取fs=1Hz,f1=0.2Hz.f2=0.3Hz
[Pxx1,f]=pmem(xn,15,Nfft,1);                               %用最大熵谱估计法估计,阶数为15,点数为1024
Pxx1=10*log10(Pxx1);                                       %取对数坐标
hold on                      
subplot(2,2,1);plot(f,Pxx1,'b');                           %画出功率谱估计的图形
xlabel('Frequency');
ylabel('Power Spectrum (dB)');
title('最大熵法 order=15,N=1024');
grid on;

[Pxx2,f]=pmem(xn,20,Nfft,1);                              %用最大熵谱估计法估计,阶数为20,点数为1024
Pxx2=10*log10(Pxx2);
hold on
subplot(2,2,2);plot(f,Pxx2,'b');
xlabel('Frequency');
ylabel('Power Spectrum (dB)');
title('最大熵法 order=20,N=1024');
grid on;

N=512;
Nfft=128;
n=[0:N-1];
randn('state',0);
wn=randn(1,N);
xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(20)*sin(2*pi*0.3*n)+wn;
[Pxx3,f]=pmem(xn,10,Nfft,1);                            %用最大熵谱估计法估计,阶数为10,点数为512
Pxx2=10*log10(Pxx2);
Pxx3=10*log10(Pxx3);
hold on
subplot(2,2,3);plot(f,Pxx3,'b');
xlabel('Frequency');
ylabel('Power Spectrum (dB)');
title('最大熵法  order=10,N=512');
grid on;

[Pxx4,f]=pmem(xn,10,Nfft,1);                            %用最大熵谱估计法估计,阶数为10,点数为1024
Pxx4=10*log10(Pxx4);
hold on
subplot(2,2,4);plot(f,Pxx4,'b');
xlabel('Frequency');
ylabel('Power Spectrum (dB)');
title('最大熵法 order=10,N=1024');
grid on;