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

功率谱密度估计方法的MATLAB实现 在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用概念，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）或者谱功率分布（spectral power distribution, SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者运用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。假如信号不是平稳过程，那么自相关函数肯定是两个变量的函数，这样就不存在功率谱密度，但是可以运用类似的技术估计时变谱密度。信号功率谱的概念和应用是电子工程的基础，尤其是在电子通信系统中，例如无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计非常重要，借助于Matlab工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修正周期图法、最大熵法、Levinson递推法和Burg法的功率谱密度估计方法进行程序设计及仿真并给出仿真结果。 以下程序运行平台Matlab R2015a（8.5.0.197613） 一、 周期图法谱估计程序 1、 源程序 Fs100000; 采样频率100kHz N1024; 数据长度N1024 n0N-1; tn/Fs; xnsin2000*2*pi*t; 正弦波,f2000Hz Yawgnxn,10; 加入信噪比为10db的高斯白噪声 subplot2,1,1; plotn,Y title信号 xlabel时间;ylabel幅度; grid on; windowboxcarlengthxn; 矩形窗 nfftN/4; 采样点数 [Pxx f]periodogramY,window,nfft,Fs; 干脆法 subplot2,1,2; plotf,10*log10Pxx; grid on; title[周期图法谱估计,,int2strN,点]; xlabel频率（Hz）;ylabel功率谱密度; 2、 仿真结果 二、 修正周期图法（加窗）谱估计程序 1、源程序 Fs100000; 采样频率100kHz N512; 数据长度 M32; 汉明窗宽度 n0N-1; tn/Fs; xnsin2000*2*pi*t; 正弦波,f2000Hz Yawgnxn,10; 加入信噪比为10db的高斯白噪声 subplot2,1,1;subplot2,1,1; plotn,Y title信号 xlabel时间;ylabel幅度; grid on; windowhammingM; 汉明窗 [Pxx f]pwelchY,window,10,256,Fs; subplot2,1,2; plotf,10*log10Pxx; grid on; title[修正周期图法谱估计 N,int2strN, M,int2strM]; xlabel频率（Hz）;ylabel功率谱密度; 2、 仿真结果 三、 最大熵法谱估计程序 1、源程序 fs1; 设采样频率 N128; 数据长度 变更数据长度会导致辨别率的变更； f10.2*fs; 第一个sin信号的频率,f1/fs0.2 f20.3*fs; 其次个sin信号的频率，f2/fs0.2或者0.3 P10; 滤波器阶数 n1N; ssin2*pi*f1*n/fssin2*pi*f2*n/fs; s为原始信号 xawgns,10; x为观测信号，即对原始信号加入白噪声，信噪比10dB figure1; 画出原始信号和观测信号 subplot2,1,1; plots,b,xlabel时间,ylabel幅度,title原始信号s; grid; subplot2,1,2; plotx,r,xlabel时间,ylabel幅度,title观测信号x; [Pxx1,f]pmemx,P,N,fs; 最大熵谱估计 figure2; plotf,10*log10Pxx1; xlabel频率Hz ;ylabel功率谱dB ; title[最大熵法谱估计 模型阶数P,int2strP, 数据长度N,int2strN]; 2、 仿真结果 四、 Levinson递推法谱估计程序 1、 源程序 fs1; 设采样频率为1 N1000; 数据长度 变更数据长度会导致辨别率的变更； f10.2*fs; 第一个sin信号的频率,f1/fs0.2 f20.3*fs; 其次个sin信号的频率，f1/fs0.2或者0.3 M16; 滤波器阶数的最大取值,超过则认为代价太大而放弃 L2*N; 有限长序列进行离散傅里叶变换前，序列补零的长度 n1N; ssin2*pi*f1*n/fssin2*pi*f2*n/fs;s为原始信号 xawgns,10;x为观测信号，即对原始信号加入白噪声，信噪比10dB figure1; 画出原始信号和观测信号 subplot2,1,1; plots,b,axis[0 100 -3 3],xlabel时间,ylabel幅度,title原始信号s; grid; subplot2,1,2; plotx,r,axis[0 100 -3 3],xlabel时间,ylabel幅度,title观测信号x; grid; 计算自相关函数 rxx xcorrx,x,M,biased;计算有偏估计自相关函数，长度为-M到M， 共2M1 r0 rxxM1; r0为零点