经典功率谱估计及其仿真
开篇：润墨网以专业的文秘视角，为您筛选了一篇经典功率谱估计及其仿真范文，如需获取更多写作素材，在线客服老师一对一协助。欢迎您的阅读与分享！
摘 要：介绍了各种经典功率谱估计方法，不仅从理论上对各种方法的谱估计质量进行了分析比较，而且通过Matlab实验仿真验证了理论分析的正确性。着重对使用比较广泛的Welch法进行了深入的研究，给出了窗函数选择的一般要求，通过仿真分析了不同的窗函数对Welch法谱估计质量的影响，比较了他们的优缺点。最后分析了采样点数较少即短数据对Welch法谱估计质量的影响。
关键词：经典谱估计；估计质量；Welch法；窗函数；短数据
中图分类号：TN911.7 文献标识码：B
文章编号：1004-373X(2008)11-159-03
Classical Power Spectrum Density Estimation and Its Simulation
SONG Ning，GUAN Hua
(School of Information & Electric Engineering,Shandong Jianzhu University,Jinan,250101,China)
Abstract：Various classical Power Spectrum Density (PSD) estimation methods are introduced,estimation quality of each method is analyzed and compared in both theory and simulation using the software Matlab.Then further study is made in Welch method which is used most widely.General selecting criterion of window function is presented and estimation quality of Welch method using different window function is compared.Finally,the impact of fewer data on estimation quality of Welch method is analyzed.
Keywords：classical PSD estimation;estimation quality;Welch method;window function;fewer data
1 引 言
信号的频谱分析是研究信号特性的重要手段之一，对于确定性信号，可以用Fourier变换来考察其频谱性质，而对于广义平稳随机信号，由于它一般既不是周期的，又不满足平方可积，严格来说不能进行Fourier变换，通常是求其功率谱来进行频谱分析。
功率谱反映了随机信号各频率成份功率能量的分布情况，可以揭示信号中隐含的周期性及靠得很近的谱峰等有用信息，应用极其广泛，例如，在语音信号识别、雷达杂波分析、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域,发挥了重要作用。
然而，实际应用中的平稳随机信号通常是有限长的，只能根据有限长信号估计原信号的真实功率谱，这就是功率谱估计问题。功率谱估计分为经典谱估计和现代谱估计，本文主要研究经典谱估计。经典谱估计又称非参数模型谱估计,主要方法有：周期图法，相关图法及改进的周期图估计法。
2 周期图法
Schuster于1899年首先提出周期图法，也称直接法，取平稳随机信号X(n)的有限个观察值x(0),x(1),…,x(N－1)对功率谱S(ω)Ы行估计：И
┆(ω)=1N|XN(ejω)|2=1N|∑N－1n=0x(n)e－jωn|2
(1)
И
主要性能指标有：
(1) 估计的均值:
И
E［[AKS^](ω)］=12πS(ω)*W(ω)
И
W(ω)是窗函数ω(n)的Fourier变换。当N∞时，E［[AKS^](ω)］S(ω)，是无偏估计；N是有限值时，是有偏估计，偏差为：
И
bia([AKS^](ω))=12πS(ω)W(ω－λ)dλ－S(ω)
И
(2) 估计的方差:
И
var([AKS^](ω))=|12πN∫π－πS(λ)D0(ω－λ)D0(ω+λ)dλ|2+
E2([AKS^](ω))
И
其中D0(ω)是矩形窗d0(n)=1，|n|≤|N－1|的Fourier变换，可见[AKS^](ω)不是S(ω)У囊恢鹿兰疲
随着N的增大，谱估计起伏增大，N∞时，var([AKS^](ω))S2(ω)。И
(3) 估计的分辨率：数据窗为长度为N的矩形窗时，Re s{S(ω)}=0.892πN，N增大时，分辨率提高，但会使[AKS^]x(ω)У钠鸱加剧，可见方差与分辨率是一对矛盾。
周期图法应用比较广泛，主要是由于它与序列的频谱有对应关系，可以采用FFT快速算法来计算。但是，这种方法需要对无限长的平稳序列进行截断，相当于对其加矩形窗，使之成为有限长数据。同时，这也意味着对自相关函数加三角窗，使功率谱与窗函数卷积，从而产生频谱泄漏，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱的模糊和失真，使得谱分辨率较低。
3 相关图法
相关图法是1958年由Blackman与Tukey首先提出的，理论基础是维纳-辛钦定理，基本思想是通过改善对相关函数的估计方法，对周期图进行平滑处理以改善周期图谱估计的方差性能，亦称BT法、间接法。首先，由信号的有限个观察值x(0),x(1),…,x(N－1)估计出自相关函数[AKR^](m)=1N∑N－1－|m|n=0x(n)x(n+m)，然后在(－M,M)内对[AKR^](m)作Fourier变换，得到功率谱[AKS^](ω)=∑Mm=－M[AKR^](m)ω(m)e－jωn，ω(n)为窗函数。
当M=N－1时，BT法与周期图法估计出的功率谱是一样的；当M
BT法的缺点在于当MN时，[AKR^](m)的方差很大，使谱估计质量下降；由[AKR^](m)得到的[AKS^](ω)Р灰欢ㄎ正值，从而可能失去功率谱的物理意义。
4 经典谱估计方法的改进
4.1 经典谱估计性能分析
平稳随机信号X(n)У淖韵喙睾数为：
И
R(n)=E［X(n)X(n+m)］
(2)
И
功率谱为：
И
S(ω)=∑∞m=－∞R(m)e－jωm
(3)
И
然而，实际应用中求不出总集意义上的自相关函数和功率谱，因而假定X(n)是各态历经的，x(n)是他的一个样本，则：
И
R(m)=limN∞12N+1∑Nn=－Nx(n)x(n+m)
(4)
S(ω)=limN∞E{12N+1|∑Nn=－Nx(n)e－jωn|2}
(5)
И
由以上两式可以看到，尽管自相关函数可以用时间平均来代替总集平均，但功率谱必须用总集平均来计算。这是因为随机过程X(n)的每一次实现x(n)，其Fourier变换仍然是一个随机过程，在每个频率ω处都是一个随机变量，因此，求数学期望是必要的，也就是说，对R(m)ёFourier变换后，S(ω)不具有各态历经性，所以，真实谱S(ω)вυ谧芗意义上求。另外，如果没有求数学期望运算，式(4)的求极限运算在任何统计意义上都不会收敛。而式(1)与式(5)相比，既无求均值计算，也无求极限运算，只能看作是对真实谱S(ω)ё骶值运算时的一个样本，缺少了统计平均，当然也就产生了方差。
为了改进周期图法的估计性能，常用的方法有两种：一是平滑，就是用适当的窗函数对谱进行平滑；二是平均，就是对同一信号做多次周期图估计后再平均，在一定程度上弥补上述所缺的求均值运算。平滑的方法在BT谱估计中已作了介绍，下面介绍两种平均方法。
4.2 Bartlett法
Bartlett谱估计法的主要目标是改善方差，基本原理是：将长度为N的数据分为L段，每段长度为M，先对每段数据用周期图法进行谱估计，然后对L段求平均得到长度为NУ氖据的功率谱。由上述原理可得功率谱为：
И
[AKS^](ω)=1L∑Li=1┆i^(ω)=1ML∑Li=1|∑M－1n=0xim(n)e－jωn|2
(6)
И
(1) 估计的均值：E［[AKS^](ω)］=12πS(ω)*W(ω)，Ш椭芷谕挤ㄒ谎。
(2) 估计的方差：当N∞时，var([AKS^](ω))1LS2(ω)，是周期图法的1/L。И
(3) 估计的分辨率：
Re s{S(ω)}=0.892πM=0.89L2πNИ
由于MNВ故其分辨率比周期图法低，可见，Bartlett法方差的改善是以牺牲分辨率为代价的。
4.3 Welch法
Welch谱估计法是对Bartlett法的改进，旨在保持Bartlett法方差性能的同时，改善其分辨率，又称加权交叠平均法，其基本原理是：对数据分段时，使每一段有部分重叠，然后对每一段数据用一个合适的窗函数进行平滑处理，最后对各段谱求平均。由上述原理可得功率谱为：
[AKS^](ω)=1MUL∑Li=1|∑M－1n=0xim(n)e－jωn|2
(7)
И
其中U=∑nω(n)，ω(n)是窗函数。
(1) 估计的均值:
E［[AKS^](ω)］=12πMUS(ω)*W(ω)ВW(ω)是Е(n)УFourier变换。
(2) 估计的方差： 当N∞时，var([AKS^](ω))9L16NS2(ω)。И
(3) 估计的分辨率：当窗函数为矩形窗、重叠50%时，ИRe s{S(ω)}=1.282πM。И
因为Welch法允许各段数据交叠，所以数据段数L会增加，使方差得到更大的改善，但是数据的交叠有减小了每一段数据的不相关性，使方差的减小不会达到理论程度。另外，采用合适的窗函数可以减小信号的频谱泄漏，同时也可以增加谱峰的宽度，从而提高谱分辨率。
5 Matlab仿真
取平稳随机信号为：
И
x(n)=2sin(2π•200•n+ω1)+
3cos(2π•300•n+ω2)+e(n)
И
其中Е1,ω2为均匀分布的随机相位，e(n)为零均值的白噪声，分别用周期图法、BT法、Bartlett法、Welch法对其进行功率谱估计。
5.1 仿真结果
在各种谱估计方法中，采样频率1 000 Hz，采样点数1 000，FFT点数512，Bartlett法中矩形窗长度为100，Welch法中Hanning窗长度为255，重叠数据为128。仿真结果如图1所示。
图1 仿真结果(一)
5.2 仿真结果分析
由图1可以直观地看出：
(1) 周期图法和BT法的特点是离散性大,曲线粗糙,方差较大,但是分辨率较高；
(2) Bartlett法和Welch法的收敛性较好,曲线平滑,方差较小,但是功率谱主瓣较宽,分辨率低，这是由于对随机序列加窗截断所引起的Gibbs效应造成的；
(3) 与Bartlett法相比,Welch法的估计曲线比较粗糙,但是分辨率较好，原因是Welch法中对数据进行截断时加的是Hanning窗，而在Bartlett法中使用的是矩形窗，相对于矩形窗，Hanning窗的主瓣包含更多的能量，因而使功率谱的主瓣较窄,分辨率较高。
由上述理论分析及仿真实验可知，Welch法采用加窗交叠求功率谱，可以有效减小方差和偏差，一般情况下能接近一致估计的要求，因而得到广泛应用。同时还可以发现，对信号加不同的窗函数，谱估计的质量是不同的，下面就通过进一步的仿真来分析研究Welch法中使用不同的窗函数时谱估计的质量。
5.3 不同窗函数的Welch谱估计
在选择窗函数时，一般有如下要求：
(1) 窗口宽度M要远小于样本序列长度N，б耘懦不可靠的自相关值；
(2) 当平稳信号为实过程时，为保证平滑周期图和真是功率谱也是实偶函数，平滑窗函数必须是实偶对称的；
(3) 平滑窗函数应当在m=0出游峰值，并且m随绝对值增加而单调下降，使可靠的自相关值有较大的权值；
(4) 功率谱是频率的非负函数且周期图是非负的，因而要求窗函数的Fourier变换是非负的。
仍以上述平稳随机信号为例，采样频率、采样点数、FFT点数、窗长度及重叠数据不变，窗函数采用矩形窗、Blackman窗、Hamming窗，仿真结果如图2所示。
图2 仿真结果(二)
由图2可以看出，使用不同的窗函数谱估计的质量是不一样的，矩形窗的主瓣较窄，分辨率较好，但方差较大，噪声水平较高；而Blackman窗和Hamming窗的主瓣较宽，分辨率较低，但方差较小，噪声水平较低。因此，在进行谱分析时选择何种窗函数，要视具体情况而定。如果强调高分辨率，能精确读出主瓣频率，而不关心幅度的精度，例如测量震动物体的自震频率，可以选用主瓣宽度比较窄的矩形窗；对受到强干扰的窄带信号，若干扰靠近信号，则可选用旁瓣幅度较小的窗函数，若离开通带较远，则可选用渐近线衰减速度比较快的窗函数。总之，要针对不同的信号和不同的处理目的来选择合适的窗函数，这样才能得到良好的效果。
5.4 短数据的Welch谱估计
以上研究了不同窗函数的Welch谱估计质量，下面研究采样数据较少时Welch谱估计的质量。仍以上述平稳随机信号为例，采样频率1 000，采样点数100，FFT点数256，采用矩形窗、Blackman窗、Hamming窗，窗长度50，重叠数据25，仿真结果如图3所示。
图3 仿真结果(三)
比较图2和图3可知，在短数据情况下，得到的信号功率谱的谱峰不够尖锐，若信号的功率更小，很容易被淹没在噪声中而无法识别。
6 结 语
在经典谱估计中，无论是周期图法还是其改进的方法，都存在着频率分辨率低、方差性能不好的问题，原因是谱估计时需要对数据加窗截断，用有限个数据或其自相关函数来估计无限个数据的功率谱，这其实是假定了窗以外的数据或自相关函数全为零，这种假定是不符合实际的，正是由于这些不符合实际的假设造成了经典谱估计分辨率较差。另外，经典谱估计的功率谱定义中既无求均值运算又无求极限运算，因而使得谱估计的方差性能较差，当数据很短时，这个问题更为突出。如何选取最佳窗函数、提高频谱分辨率，如何在短数据情况下提高信号谱估计质量，还需要进一步研究。
参 考 文 献
［1］金连文，韦岗.现代数字信号处理简明教程［M］.北京：清华大学出版社，2004.
［2］飞思科技产品研发中心.Matlab7辅助信号处理技术与应用［M］.北京：电子工业出版社，2005.
［3］王晓峰，王炳和.周期图及其改进方法中谱分析率的Matlab分析［J］.武警工程学院学报,2003(6):75-77.
［4］姚武川，姚天任.经典谱估计方法的Matlab分析［J］.华中理工大学学报,2000,28(4)：45-47.
作者简介 宋 宁 男，1983年出生，硕士研究生。研究方向为自动化装置的集成化与智能化。
关 华 男，1963年出生，硕士生导师，博士，副教授。主要从事信号处理、图像处理、DSP等方面的研究。
注：本文中所涉及到的图表、注解、公式等内容请以PDF格式阅读原文。