(完整word版)功率谱密度估计方法的MATLAB实现

功率谱密度估计方法的MATLAB实现

在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用概念，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）或者谱功率分布（spectral power distribution, SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者使用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。信号功率谱的概念和应用是电子工程的基础，尤其是在电子通信系统中，例如无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计十分重要，借助于Matlab工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修正周期图法、最大熵法、Levinson递推法和Burg法的功率谱密度估计方法进行程序设计及仿真并给出仿真结果。

以下程序运行平台：Matlab R2015a（8.5.0.197613）

一、周期图法谱估计程序

1、源程序

Fs=100000; %采样频率100kHz

N=1024; %数据长度N=1024

n=0:N-1;

t=n/Fs;

xn=sin(2000*2*pi*t); %正弦波,f=2000Hz

Y=awgn(xn,10); %加入信噪比为10db的高斯白噪声

subplot(2,1,1);

plot(n,Y)

title('信号')

xlabel('时间');ylabel('幅度');

window=boxcar(length(xn)); %矩形窗

nfft=N/4; %采样点数

[Pxx f]=periodogram(Y,window,nfft,Fs); %直接法

subplot(2,1,2);

plot(f,10*log10(Pxx));

grid on;

title(['周期图法谱估计,',int2str(N),'点']);

xlabel('频率（Hz）');ylabel('功率谱密度');

2、仿真结果

二、修正周期图法（加窗）谱估计程序

1、源程序

Fs=100000; %采样频率100kHz N=512; %数据长度

M=32; %汉明窗宽度

n=0:N-1;

xn=sin(2000*2*pi*t); %正弦波,f=2000Hz

Y=awgn(xn,10); %加入信噪比为10db的高斯白噪声subplot(2,1,1);subplot(2,1,1);

plot(n,Y)

title('信号')

xlabel('时间');ylabel('幅度');

grid on;

window=hamming(M); %汉明窗

[Pxx f]=pwelch(Y,window,10,256,Fs);

subplot(2,1,2);

plot(f,10*log10(Pxx));

grid on;

title(['修正周期图法谱估计N=',int2str(N),' M=',int2str(M)]); xlabel('频率（Hz）');ylabel('功率谱密度');

2、仿真结果

三、最大熵法谱估计程序

1、源程序

fs=1; %设采样频率

N=128; %数据长度改变数据长度会导致分辨率的变化；

f1=0.2*fs; %第一个sin信号的频率,f1/fs=0.2

f2=0.3*fs; %第二个sin信号的频率，f2/fs=0.2或者0.3

P=10; %滤波器阶数

n=1:N;

s=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs); %s为原始信号

x=awgn(s,10); %x为观测信号，即对原始信号加入白噪声，信噪比10dB figure(1); %画出原始信号和观测信号

subplot(2,1,1);

plot(s,'b'),xlabel('时间'),ylabel('幅度'),title('原始信号s');

grid;

subplot(2,1,2);

plot(x,'r'),xlabel('时间'),ylabel('幅度'),title('观测信号x');

[Pxx1,f]=pmem(x,P,N,fs); %最大熵谱估计

figure(2);

plot(f,10*log10(Pxx1));

xlabel('频率(Hz) ');ylabel('功率谱(dB) ');

title(['最大熵法谱估计模型阶数P=',int2str(P),' 数据长度N=',int2str(N)]);

2、仿真结果

四、L evinson递推法谱估计程序

1、源程序

fs=1; %设采样频率为1

N=1000; %数据长度改变数据长度会导致分辨率的变化；

f1=0.2*fs; %第一个sin信号的频率,f1/fs=0.2

f2=0.3*fs; %第二个sin信号的频率，f1/fs=0.2或者0.3

M=16; %滤波器阶数的最大取值,超过则认为代价太大而放弃

L=2*N; %有限长序列进行离散傅里叶变换前，序列补零的长度

n=1:N;

s=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%s为原始信号

x=awgn(s,10);%x为观测信号，即对原始信号加入白噪声，信噪比10dB figure(1); %画出原始信号和观测信号

subplot(2,1,1);

plot(s,'b'),axis([0 100 -3 3]),xlabel('时间'),ylabel('幅度'),title('原始信号s'); grid;

subplot(2,1,2);

plot(x,'r'),axis([0 100 -3 3]),xlabel('时间'),ylabel('幅度'),title('观测信号x'); grid;

%计算自相关函数

rxx = xcorr(x,x,M,'biased');%计算有偏估计自相关函数，长度为-M到M，

%共2M+1

r0 = rxx(M+1); %r0为零点上的自相关函数，相对于-M，第M+1个点为零点R = rxx(M+2:2*M+1);% R为从1到第M个点的自相关函数矩阵

%确定矩阵大小

a = zeros(M,M);

FPE = zeros(1,M);%FPE：最终预测误差，用来估计模型的阶次

var = zeros(1,M);

%求初值

a(1,1) = -R(1)/r0;%一阶模型参数

var(1) = (1-(abs(a(1,1)))^2)*r0;%一阶方差