功率谱估计及其MATLAB仿真

1、功率谱估计及其MATLAB仿真王凤瑛,张丽丽(山东科技大学信息与电气工程学院,山东青岛 266510摘要:从介绍功率谱的估计原理入手分析了经典谱估计和现代谱估计两类估计方法的原理、各自特点及在Matlab中的实现方法。关键词:功率谱估计;周期图法;AR参数法;Matlab中图分类号:TP391.9 文献标识码:APower Spectrum Density Estimation and the simulation in MatlabWANG Feng-ying ,Zhang Lili(College of Info and Elect. Eng., SUST, Qingdao, Shando

2、ng 266510, ChinaAbstract:Mainly introduces the principles of classical PSD estimation and modern PSD estimation, discusses the characteristics of the methods of realization in Matlab.Moreover,It gives an example of each part in realization using Matlab functions.Keywords:PSD Pstimation, Periodogram

3、method, AR Parameter method, Matlab现代信号分析中,对于常见的具有各态历经的平稳随机信号,不可能用清楚的数学关系式来描述,但可以利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD。它是数字信号处理的重要研究内容之一。功率谱估计可以分为经典功率谱估计(非参数估计和现代功率谱估计(参数估计。功率谱估计在实际工程中有重要应用价值,如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域,发挥了重要作用。Matlab是Math

4、Works公司于1982年推出的一套高性能的数值计算和可视化软件,人称矩阵实验室,它集数值分析、矩阵运算、信号处理和图形显示于一体,构成了一个方便的、界面友好的用户环境,成为目前极为流行的工程数学分析软件。也为数字信号处理进行理论学习、工程设计分析提供了相当便捷的途径。本文的仿真实验中,全部在Matlab6.5环境下调试通过;随机序列由频率不同的正弦信号加高斯白噪声组成。1 经典功率谱估计经典功率谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗。经典功率谱估计方法分为:相关函数法(BT法、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法,其中周期图法应用较多,具有代表

5、性。1.1 相关函数法(BT法该方法先由序列x(n估计出自相关函数R(n,然后对R(n进行傅立叶变换,便得到x(n的功率谱估计。当延迟与数据长度相比很小时,可以有良好的估计精度。Matlab代码示例1:Fs=500;%采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*40*n+3*cos(2*pi*90*n+randn(size(n;nfft=512;cxn=xcorr(xn,'unbiased' %计算序列的自相关函数CXk=fft(cxn,nfft;Pxx=abs(CXk;index=0:round(nfft/2-1;k=index*Fs/nfft;

6、plot_Pxx=10*log10(Pxx(index+1;figure(1%plot(k,plot_Pxx;1.2 周期图法(periodogram周期图法是把随机序列x(n的N个观测数据视为一能量有限的序列,直接计算x(n的离散傅立叶变换,得X(k,然后再取其幅值的平方,并除以N,作为序列x(n真实功率谱的估计。Matlab代码示例2:Fs=600; %采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*40*n+3*cos(2*pi*90*n+0.1*randn(size(n;window=boxcar(length(xn;%矩形窗nfft=512;Pxx,f=pe

7、riodogram(xn,window,nfft,Fs;%直接法plot(f,10*log10(Pxx;window=boxcar(length(xn;%矩形窗nfft=1024;Pxx,f=periodogram(xn,window,nfft,Fs; %直接法figure(1plot(f,10*log10(Pxx;1.3 平均周期图法和平滑平均周期图法对于周期图的功率谱估计,当数据长度N太大时,谱曲线起伏加剧,若N太小,谱的分辨率又不好,因此需要改进。两种改进的估计法是平均周期图法和平滑平均周期图法。Bartlett法:Bartlett平均周期图的方法是将N点的有限长序列x(n分段求周期图再

8、平均。Matlab代码示例3:fs=600;n=0:1/fs:1;xn=cos(2*pi*20*n+3*cos(2*pi*90*n+randn(size(n;nfft=512;window=hamming(nfft; %矩形窗noverlap=0;%数据无重叠p=0.9;%置信概率Pxx,Pxxc=psd(xn,nfft,fs,window,noverlap,p;index=0:round(nfft/2-1;k=index*fs/nfft;plot_Pxx=10*log10(Pxx(index+1;plot_Pxxc=10*log10(Pxxc(index+1;figure(1plot(k,p

9、lot_Pxx;figure(2plot(k,plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc;Welch法:Welch法对Bartlett法进行了两方面的修正,一是选择适当的窗函数w(n,并在周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时,可使各段之间有重叠,这样会使方差减小。Matlab 代码示例4: Fs=600; n=0:1/Fs:1;xn=cos(2*pi*40*n+3*cos(2*pi*90*n+randn(size(n; nfft=512;window=boxcar(100;%矩形窗 window1

10、=hamming(100;%海明窗window2=blackman(100;%blackman 窗 noverlap=20; %数据无重叠range='half' %频率间隔为0 Fs/2,计算一半的频率 Pxx,f=pwelch(xn,window,noverlap,nfft,Fs,range; Pxx1,f=pwelch(xn,window1,noverlap,nfft,Fs,range; Pxx2,f=pwelch(xn,window2,noverlap,nfft,Fs,range; plot_Pxx=10*log10(Pxx; plot_Pxx1=10*log10(Px

11、x1; plot_Pxx2=10*log10(Pxx2; figure(1plot(f,plot_Pxx; figure(2plot(f,plot_Pxx1; figure(3plot(f,plot_Pxx2;2 现代功率谱估计现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。主要方法有最大嫡谱分析法(AR 模型法、Pisarenko 谐波分解法、Prony 提取极点法、Prony 谱线分解法以及Capon 最大似然法等。其中AR 模型应用较多,具有代表性。常用的模型有ARMA 模型、

12、AR 模型、MA 模型。2.1模型的公式表达ARMA 模型功率谱数学式子表达为:222(1/1ppj kj kkj kxk lk leb ea eP=+其中2是激励白噪声的方差,(j x e P 为功率谱密度,a k 和b k 为模型参数。如果ARMA 模型参数b 1,b 2b q 全为0,就演化为AR 模型:22(/1pj k j kxk l ea eP =+如果ARMA 模型参数a 1,a 2a q 全为0,就演化为MA 模型:22(/1pj kj kxk le beP =+在实际中,AR 模型的参数估计比较简单,对其有充分的研究,而对于ARMA 模型,其参数比较复杂,对其算法的研究和改进

13、还在完善中。2.2 AR 模型功率谱估计AR模型功率谱估计又称为自回归模型,它是一个全极点的模型,要利用AR模型进行功率谱估计须通过levinson_dubin递推算法由Yule-Walker方程求得AR的参数:212p。在Matlab 仿真中可调用Pburg函数直接画出基于burg算法的功率谱估计的曲线图。Matlab代码示例5:用周期图法求出的功率谱曲线和burg算法求出的AR功率谱曲线(p=50 fs=200;n=0:1/fs:1;xn=cos(2*pi*40*n+cos(2*pi*41*n+3*cos(2*pi*90*n+0.1*randn(size(n;window=boxcar(l

14、ength(xn;nfft=512;pxx,f=periodogram(xn,window,nfft,fs;subplot(121plot(f,10*log10(pxxxlabel('frequency(hz'ylabel('power spectral density(Db/Hz'title('periodogram psd estimate'order1=50;range='half'magunits='db'subplot(122pburg(xn,order1,nfft,fs,range周期图法求出的功率谱曲

15、线和burg算法求出的AR功率谱曲线(p=50如图所示: 3 常见谱估计法的比较通过实验仿真可以直观地看出以下特性:(1功率谱估计中的相关函数法和周期图法所得到的结果是一致的,其特点是离散性大,曲线粗糙,方差较大,但是分辨率较高。(2平均周期图法和平滑平均周期图法的收敛性较好,曲线平滑,估计的结果方差较小,但是功率谱主瓣较宽,分辨率低。这是由于对随机序列的分段处理引起了长度有限所带来的Gibbs现象而造成的。(3平滑平均周期图法与平均周期图法相比,谱估值比较平滑,但是分辨率较差。其原因是给每一段序列用适当的窗口函数加权后,在得到平滑的估计结果的同时,使功率谱的主瓣变宽,因此分辨率有所下降。经典

16、功率谱估计的分辨率反比于有效信号的长度,但现代谱估计的分辨率可以不受此限制。这是因为对于给定的N点有限长序列x(n,虽然其估计出的自相关函数也是有限长的,但是现代谱估计的一些隐含着数据和自相关函数的外推,使其可能的长度超过给定的长度,不象经典谱估计那样受窗函数的影响。因而现代谱的分别率比较高,而且现代谱线要平滑得多,从上图可以清楚看出。4 结束语近十几年来,随着计算机的广泛应用,各种新算法不断提出,谱估计得到了迅速的发展,出现了许多功率谱估计的实现方法,也有很多具体的算法。Matlab提供的算法函数为功率谱估计提供了一条可行的方便途径,但较为有限。在熟悉了估计原理之后,可以自己动手编写m文件来实现。本文作者创新点:分析了谱估计法的原理、各自特点,实现Matlab仿真,并且附有