基于MATLAB实现经典功率谱估计

基于MATLAB实现经典功率谱估计*王春兴(山东师范大学物理与电子科学学院，250014，山东省济南市)摘要:功率谱估计(powerspectrumestimation，PSE)是利用有限长的数据估计信号的功率谱，广泛应用于雷达、声纳、通信、地质勘探、天文、生物医学工程等众多领域．功率谱估计可分为经典谱估计和现代谱估计，文章主要研究经典谱估计中的修正算法．目前经典谱估计中主要采用对2种方法进行估计，分别是周期图法和自相关法，而修正算法就是在上述方法的基础上进行改进得到的．关键词:自相关函数;功率谱估计;窗函数中图分类号:TN911．6文献标识码:A文章编号:1001-5337(2011)02-0059-041引言随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应用功率信号来描述．然而，功率信号不满足傅里叶变换的狄里赫莉绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的．因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱．信号的功率谱密度描述随机信号的功率在频域随频率的分布．利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计．谱估计方法分为2大类:非参数化方法和参数化方法．非参数化方法又叫做经典谱估计，如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低;而参数化谱估计又叫做现代谱估计，如AR模型法、移动平均模型法(简称MA模型法)、自回归移动平均模型法(简称ARMA模型法)等．2经典功率谱估计经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为零，这就相当将数据加一窗函数．根据截取的N个样本数据用估计出其功率谱．其中可以利用相关函数估计功率谱、也可以利用周期图法估计出功率谱．这些方法实质上依赖于傅立叶变换，因此实现较容易，且可以采用FFT等技术使计算量大大降低而收到人们青睐．但由于利用2种方法得到的功率谱方差性能不理想，可以对其进行修正改进．2．1用相关函数获得功率谱［8，10，11，12］根据维纳—辛钦定理，对于离散随机信号有:Sx(ejω)=∑∞m=－∞Rx(m)e－jωm，Rx(m)=12π∫π－πSx(ejω)ejωmdω．上式中RX(m)为离散随机信号的自相关函数，Sx(ejω)为功率谱密度．如果获得随机信号的自相关函数，那么功率谱密度也就不难求出．因此，这种方法的实质是对相关函数的估值求傅立叶变换即为功率谱密度．对于平稳随机信号来说，其相关函数是确定性函数，故其功率谱也是确定的．这样可由平稳随机离散信号的有限个离散值x(0)，x(1)，…，x(N－1)求出自相关函数RX(m)=1N∑N－m－1n=0x(n)x(n+m)．(1)然后在(一M，M)内对RX(m)作Fourier变换，得到功率谱∑Mm=－MRX(m)e－jωm．(2)第37卷第2期2011年4月曲阜师范大学学报JournalofQufuNormalUniversityVol．37No．2Apr．2011*收稿日期:2010-11-18基金项目:济南市高校自主创新项目资助．作者简介:王春兴，男，1962-，博士，副教授;研究方向:智能信号处理，信息安全;E-mail:[e-mail]．上式可用FFT算法求出．当M趋向于N时，由上式求得的相关函数RX(m)的方差较大，也即相关函数的可靠性差，使谱估计质量下降．为了弥补这一缺陷，在RX(m)上加一ω(m)为窗函数，从而使方差大的自相关函数估值加权小，以减少对谱估计的影响．∑Mm=－MRX(m)ω(m)e－jωm，(3)式中窗函数ω(m)=ω(m)，－M≤m≤M，0，其他{．(4)由于功率谱是非负值，那么加窗后的功率谱也应当是非负值．因此窗函数的傅立叶变换必须是非负值，这样窗函数必须是偶对称的，三角形窗是经常用的．ω(m)=1－|m|M，|m|≤M，M＜N，0，其他{．(5)该方法的缺点在于当M趋向于N时，相关函数R(m)的方差很大，使谱估计质量下降;由R(m)得到的S(w)不一定为正值，从而可能失去功率谱的物理意义．该方法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计．当延迟与数据长度相比很小时，可以有良好的估计精度(图1)．图1用相关函数获得的功率谱基于Matlab实现的程序:Fs=500;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;cxn=xcorr(xn，＇unbiased＇);CXk=fft(cxn，nfft);Pxx=abs(CXk);index=0:round(nfft/2－1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));figure(3);plot(k，plot_Pxx);2．2周期图法(Periodogram)Schuster首先提出周期图法．周期图法是根据各态历经的随机过程功率谱的定义进行的普估计．周期图法是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得x(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计．取平稳随机信号X(n)的有限个观察值x(0)，x(1)，…，x(N－1)，求出其傅立叶变换XN(e－jω)=∑N－1n=0x(n)e－jωn，(6)然后进行谱估计S(ω)=1N|XN(e－jω)|2．(7)当M=N－1时，相关函数法与周期图法估计出的功率谱是相同的;当M＜N－1时，相关函数法的偏差大于周期图法，在窗函数满足一定条件时是渐进无偏估计;方差小于周期图的方差;分辨率比周期图法低，与窗函数的选择有关．可见，窗函数的选择对谱估计的修正程度是不一样的．周期图法应用比较广泛，主要是由于它与序列的频谱有直接的对应关系，并且可以采用FFT快速算法来计算．但是，这种方法需要对无限长的平稳随机序列进行截断，相当于对其加矩形窗，使之成为有限长数据．同时，这也意味着对自相关函数加三角窗，使功率谱与窗函数卷积，从而产生频谱泄漏，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱的模糊和失真，使得谱分辨率较低(图2)．该方法基于Matlab实现的程序: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］=periodogram(xn，window，nfft，Fs);figure(1);plot(f，10*log10(Pxx));window=boxcar(length(xn));06曲阜师范大学学报(自然科学版)2011年nfft=1024;［Pxx，f］=periodogram(xn，window，nfft，Fs);figure(2);plot(f，10+log10(Pxx));图2用周期图法获得的功率谱2．3经典谱估计的改进从上面的分析知，周期图法不满足一致估计的条件，必须进行改进，采用的措施主要是将周期图进行平滑，使估计方差减小，从而得到一致谱估计．对于用相关函数法进行谱估计以及修正方法已在前面做了介绍．对于周期图的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进．2种改进的估计法是平均周期图法和平滑平均周期图法，主要目标是改善方差．Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均．将长度为N的数据分为L段，每段长度为M．先对每段数据用周期图法进行谱估计，然后对L段求平均得到长度为N的数据的功率谱．可得功率谱为S(ω)=1ML∑Li=1∑Mn=0xim(n)e－jωn2，(8)该估计的分辨率Res{S(w)}=0．892pM=0．89L2pN，(9)式中由于M远小于N，故其分辨率比周期图法低．可见，Bartlett法方差的改善是以牺牲分辨率为代价的．Welch法谱估计是在上述基础上进行了改进，目的是在保持Bartlett法方差性能的同时，改善其分辨率．其基本原理是先对随机序列分段时，使每一段有部分重叠，然后对每一段数据用一个合适的窗函数进行平滑处理，最后对各段谱求平均．这样可得功率谱S(ω)=1MUL∑Li=1∑M－1n=0xim(n)e－jωn2，(10)其中U=∑M－1n=0ω(n)，ω(n)是窗函数．该估计的分辨率:当窗函数为矩形窗、重叠50%时，Res{S(w)}=1．282pM．(11)因为Welch法允许各段数据交叠，所以数据段数L会增加，使方差得到更大的改善，但是数据的交叠有减小了每一段数据的不相关性，使方差的减小不会达到理论程度．另外，采用合适的窗函数可以减少信号的频谱泄露，同时也可以增加谱峰的宽度，从而提高分辨率(图3)．图3修正后得到的功率谱基于MATLAB实现的程序为: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);%矩形窗windowl=hamming(100);%海明窗window2=blackman(100);%blackman窗16第2期王春兴:基于MATLAB实现经典功率谱估计noverlap=20;%数据无重叠range=＇half＇;%频率间隔为［0Fs/2］，计算一半的频率［Pxx，f］=pwelch(xn，window，noverlap，nfft，Fs，range);［Pxxl，f］=pwelch(xn，windowl，noverlap，nfft，Fs，range);［Pxx2，f］=pwelch(xn，window2，noverlap，nfft，Fs，range);plot_Pxx=10*log10(Pxx);plot_Pxxl=10*log10(Pxxl);plot_Pxx2=10*log10(Pxx2);figure(1)plot(f，plot_Pxx);figure(2)plot(f，plot_Pxxl);figure(3)plot(f，plot_Pxx2);plot_Pxxl=10*logl0(Pxxl);plot_Pxx2=10*logl0(Pxx2);figure(1)plot(f，plot_Pxx);figure(2)plot(f，plot_Pxxl);figure(3)plot(f，plot_Pxx2);3结论对于平稳随机过程来说，功率谱理论上的数值是不可能实现的，只能用有限的观测数据来估计或者逼近真实值，从而比较真实地反映问题的本质．当然估计结果的好坏，与拟合的数学模型及所采用的处理方法等都有很大关系．参考文献:［1］姚武川，姚天任．经典谱估计方法的Matlab分析［J］．华中理工大学学报，2000，28(4):45-47．［2］王晓峰，王炳和．周期图及其改进方法中谱分析率的Mat-lab分析［J］．武警工程学院学报，2003，19(6):75-77．［3］宋宁，关华．经典功率谱估计及其仿真I－J1．现代电子技术，2008，2008，(11):159-164．［4］冯磊．经典功率谱估计与现代功率谱估计的对比［J］．商业文化，2009，2009，(5):223-239．［5］宁长春，陈天禄，索郎桑姆，等．数字信号处理中常用的Matlab工具箱函数简介［J］．西藏科技，2007，(12):75-78．［6］魏鑫，张平．周期图法功率谱估计中的窗函数分析［J］．现代电子技术，2005，(3):14-15．［7］邵玉斌．Matlab/Simulink通信系统建模与仿真实例分析［M］．北京:清华大学出版社，2008．［8］范瑜，邬正义．功率谱估计的Welch方法中的窗函数研