经典功率谱估计方法实现问题的研究

PAGE231随机信号的经典谱估计方法 估计功率谱密度的平滑周期图是一种计算简单的经典方法。它的主要特点是与任何模型参数无关，是一类非参数化方法[4]。它的主要问题是：由于假定信号的自相关函数在数据观测区以外等于零，因此估计出来的功率谱很难与信号的真实功率谱相匹配。在一般情况下，周期图的渐进性能无法给出实际功率谱的一个满意的近似，因而是一种低分辨率的谱估计方法。本章主要介绍了周期图法、相关法谱估计（BT）、巴特利特(Bartlett)平均周期图的方法和Welch法这四种方法。2.1周期图法周期图法又称直接法。它是从随机信号x(n)中截取N长的一段，把它视为能量有限x(n)真实功率谱的估计的抽样.周期图这一概念早在1899年就提出了，但由于点数Ｎ一般比较大，该方法的计算量过大而在当时无法使用。只是1965年FFT出现后，此法才变成谱估计的一个常用方法。周期图法[5]包含了下列两条假设：1．认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段来估计该随机序列的功率谱。这当然必然带来误差。2．由于对采用DFT，就默认在时域是周期的，以及在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段的周期延拓，这也就是周期图法这个名字的来历。与相关法相比，相关法在求相关函数时将以外是数据全都看成零，因此相关法认为除外x(n)是全零序列，这种处理方法显然与周期图法不一样。但是，当相关法被引入基于FFT的快速相关后，相关法和周期图法开始融合。通过比较我们发现：如果相关法中M=N，不加延迟窗，那么就和补充（N-1）个零的周期图法一样了。简单地可以这样说：周期图法是M=N时相关法的特例。因此相关法和周期图法可结合使用。2.2相关法谱估计（BT）法这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。它是1958年由Blackman和Tukey提出。这种方法的具体步骤是：第一步：从无限长随机序列x(n)中截取长度N的有限长序列列第二步：由N长序列求（2M-1）点的自相关函数序列。即(2-1)这里，m=-(M-1)…,-1,0,1…,M-1，MN，是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。。。,M-1的傅里叶变换，另一半也就知道了。第三步：由相关函数的傅式变换求功率谱。即(2-2)以上过程中经历了两次截断，一次是将x(n)截成N长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的代表估值。一般取M<<N，因为只有当M较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。因此，在FFT问世之前，相关法是最常用的谱估计方法。当FFT问世后，情况有所变化。因为截断后的可视作能量信号，由相关卷积定理可得(2-3)这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。我们可对上式两边取（2N-1）点DFT，则有(2-4)于是将时域卷积变为频域乘积，用快速相关求的完整方案如下：对N长的补充（N-1）个零，成为（2N-1）长的。求（2N-1）点的FFT，得。求。由DFT性质，是纯实的，满足共轭偶对称，而一定是实偶的，且以（2N-1）为周期。求（2N-1）点的IFFT：(2-5)这里是实偶的，m=-(N-1)...0...N-1。本来IFFT求和范围是0至2N-2，由于的实偶性与周期性，求和范围改为-（N-1）至（N-1）不影响计算结果。同理可将m的范围改为-（N-1）至（N-1）。上述的快速相关中，补充零的目的是为了能用圆周卷积代替线性卷积，以便进一步采用快速卷积算法。快速相关输出是-（N-1）至（N-1）的2N-1点，加窗后截取的是-（M-1）至（M-1）的频段，最后作（2M-1）点FFT，得。我们注意到：如果数据点数与自相关序列点数相同即M=N，则（2N-1）点的IFFT后紧跟一个（2N-1）点的FFT，利用的对称性，FT运算框的计算式变为(2-6)由于N=M并假设窗形状是矩形的，第二次的截断就不需要了。比较式（2-5）和式（2-6），,正反傅氏变换可以抵消，直接得(2-7)为了实行基2FFT，也可将（2N-1）点换成2N点，这样做不影响结果的正确性。2.3巴特利特(Bartlett)平均周期图法首先让我们来看一下为什么周期图经过某种平均(或平滑)后会使它的方差当时趋于零，达到一致估计的目的。如果是不相关的随机变量，每一个具有期望值，方差，则可以证明它们的数学平均的期望值等于，数学平均的方差等于，即：所以 (2-8)由(2-8)可见，L个平均的方差比每个随机变量的单独方差小L倍。当，可达到一致谱估计的目的。因而降低估计量的方差的一种有效方法是将若干个独立估计值进行平均。把这种方法应用于谱估计通常归功于Bartlett。Bartlett平均周期图的方法是将序列分段求周期图再平均。设将分成L段，每段有M个样本，因而，第i段样本序列可写成第i段的周期图为如果很小，则可假定各段的周期图是互相独立的。对功率谱密度的概念的讨论，谱估计可定义为L段周期图的平均，即(2-9)于是它的期望值为(2-10)这里，因此Bartlett估计的期望值是真实谱与三角窗函数的卷积。由于三角窗函数不等于函数，所以Bartlett估计也是有偏估计即，但当时，。由于我们假定各段周期图是相互独立的，所以可按式(2-8)得到下式：(2-11)由此可见，随着L的增加是下降的，当时，。因此Bartlett估计是一致估计。比较式(2-10)的与式(2-1)的可见在二种情况的估计量的期望值都是真值与窗口函数的卷积形式，但后者将前者WB中N改为M，。因而使主瓣的宽度增大。由于主瓣的宽度愈窄愈接近函数，偏倚愈小。今式(3-10)中的主瓣宽度大于后式中的主瓣宽度，因而，而主瓣愈宽分辨率就愈差。因此Bias可用来说明谱的分辨率，Bias愈大说明谱分辨率愈差。一个固定的记录长度N，周期图分段的数目L愈大将使方差愈小，但M也愈小，因而使Bias愈大，谱分辨率变得愈差。因此Bartlett方法中Bias或谱分辨率和估计量的方差间是有互换关系的。M和N的选择一般是由对所研究的信号的预先了解来指导的。例如，如果我们知道谱有一个窄峰，同时如果分辨出这个峰是重要的，那么我们必须选择M足够大。又从方差的表达式我们可以确定谱估计的可接受的方差所要求的记录长度N=(LM)。由此可见Bartlett法使谱估计的方差减小是用增加Bias以及降低谱分辨率的代价换来的。实际上，当N一定时，Bias与Var的互换性是谱估计的一个固有特性。[例如]为了说明经平均后的周期图作为功率谱估计的实际效果，设有一零均值高斯分布的随机过程，其功率谱密度为这一功率谱密度是由一零值高斯分布单位方差的噪声序列通过一个其的滤波器形成的[6]。为了简单设选用的矩形窗函数。说明的周期图可以得到的Bias的情况。图2表示M=8分4段与16段二种情况平均后的周期图。显然L=4的方差比L=16的大。M=16，L=2及8的周期图表示在图3。图3中L不同造成的影响也是明显的。但是这二个曲线的起伏都很大，因此有理由认为误差主要起因于方差。比较与M=16的周期图可见，M愈大，将使周期图的起伏愈增快的结论，在这里也同样成立。比较图2与图3发现在这个例子中最好的选择是应用L=16，M=8的估计而不是L=8、M=16的估计，即宁可减小方差，牺牲Bias。在实际中，当然功率谱密度的真值是不知道的。但是谱的窗口函数以及关于功率谱密度的某些信息往往是预先知道的。通过改变M和L以及利用预先知道的情况，通常可以很好地进行选择。平均周期图的方法特别适合于应用FFT算法。因此在FFT出现以后这个方法比下面将要讨论的利用窗口函数的处理法用得更多。而在FFT出现以前主要是用窗口函数处理法平滑周期图。图1与的特性[6]图2平滑后的周期图[6]（每段取8个数据）图3平均后的周期图[6](每段取16个数据)2.4Welch法Welch提出了对Bartlett法的修正使更适合于用FFT进行计算。他主要提出二方面的修正，其一是选择适当的窗函数，并在周期图计算前直接加进法，这样得到的每一段的周期图为(2-12)这里为归一化因子，而Bartlett法每段的周期图为(2-13)这样加窗函数的优点是无论什么样的窗函数均可使谱估计非负。其二是在分段时，可使各段之间有重迭，这样将会使方差减小(当N与M一定时)，重迭可以达到50％。3经典谱估计方法的仿真和分析要对谱估计方法进行分析和比较，需要用到功能强大的matlab编程软件[7]对这些方法编程仿真。MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。3.1周期图法仿真周期图法又称直接法，它是把随机序列x(n)的N个观测数据视为一个能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。Matlab代码示例：clear;Fs=1000;%采样频率n=0:1/Fs:1;%产生含有噪声的序列改变n的取值范围观察图形的变换xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));window=boxcar(length(xn));%矩形窗nfft=1024;[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法plot(f,10*log10(Pxx));xlabel('频率')ylabel('功率/DB')图4周期图当N=100功率谱图5周期图法中N=1000功率谱图6周期图中N=100000功率谱通过matlab程序仿真从图4至图6可以看出随着采样点数的增加，该估计是渐进无偏的。从图5和图6可以看出，采用周期突发估计得出的功率谱很不平滑，相应的估计协方差比较大。而且采用增加采样点的办法也不能吃周期图变得更加平滑，这是周期图法的缺点。周期图法得出的估计谱方差特性不好：当数据长度N太大时，扑线的起伏加剧；N太小时谱的分辨率又不好。对其改进的主要方法有二种，即平均和平滑，平均就是将截取的数据段再分成L个小段，分别计算功率谱后取功率谱的平均，这种方法使估计的方差减少，但偏差加大，分辨率下降。平滑是用一个适当的窗函数与算出的功率谱进行卷积，使谱线平滑。这种方法得出的谱估计是无偏的，方差也小，但分辨率下降。3.2相关法谱估计法仿真间接法是先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计的方法。功率谱图如图7所示。Matlab代码示例：clear;Fs=1000;%采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));nfft=1024;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));plot(k,plot_Pxx);xlabel('频率')ylabel('功率/DB')图7直接法功率谱图当M=N-1时,BT法与周期图法估计出的功率谱是一样的;当M<N-1时,BT法的偏差大于周期图法,在窗函数满足一定条件时是渐进无偏估计;方差小于周期图的方差;分辨率比周期图法低,与窗函数的选择有关。BT法的缺点在于当M→N时,的方差很大,使谱估计质量下降;由得到的不一定为正值,从而可能失去功率谱的物理意义。3.3改进的周期图法仿真对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。对周期图法改进的思想是将信号分段进行估计，然后再将这些估计结果进行平均，从而减小估计的协方差，是功率谱图变得比直接法更平滑。增加分段数可以进一步降低估计的协方差，然而若每段中的数据点数太少，就会使估计的频率分辨率下降很多[8]。从样本信号序列总点数一定的条件下，可以采用使分段相互重叠的方法来增加分段数，从而保持每段信号点数不变，这样就在保证频率分辨率的前提下进一步降低估计协方差。主要的改进方法有Bartlett法和Welch法3.3.1Bartlett法Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。功率谱图如图8所示。

Matlab代码示例：clear；Fs=1000;n=0:1/Fs:1;xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));nfft=1024;window=boxcar(length(n));%矩形窗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(1)plot(k,plot_Pxx);pause;figure(2)plot(k,[plot_Pxxplot_Pxx-plot_Pxxcplot_Pxx+plot_Pxxc]);xlabel('频率')ylabel('功率/DB')

图8bartlett法估计功率谱图3.3.2Welch法现在比较常用的改进方法是Welch法，又叫加权交叠平均法[9]，简记为WOSA法，这种方法以加窗（加权）求取平滑，以分段重叠求得平均，因此集平均与平滑的优点于一体，同时也不可避免带有两者的缺点，因此归根到底是一种折中。其主要步骤是：（1）将N长的数据段分成L个小段，每小段M点，相邻小段间交叠M/2点（即2:1分段）。因为L（M/2）+M/2=N,所以段数(3-1)（2）对各小段加同样的平滑窗后起傅氏变换（3-2）（3）用下式求各小段功率谱的平均（3-3）这里，代表窗函数平均功率，MU是M长窗函数的能量。仍以上述平稳随机信号为例,采样频率、采样点数、FFT点数、窗长度及重叠数据不变,窗函数采用矩形窗、Blackman窗、Hamming窗，仿真结果如图所示。图9是加矩形窗得到的功率谱图；图10是加海明窗所得的功率谱图；图11是加blackman窗所得到的功率谱图。Matlab代码示例：clear; Fs=1000;n=0:1/Fs:1;xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));nfft=1024;window=boxcar(100);%矩形窗window1=hamming(100);%海明窗window2=blackman(100);%blackman窗noverlap=20;%数据无重叠range='half';%频率间隔为[0Fs/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(Pxx1);plot_Pxx2=10*log10(Pxx2);figure(1)plot(f,plot_Pxx);xlabel('频率')ylabel('功率/DB')pause;figure(2)plot(f,plot_Pxx1);xlabel('频率')ylabel('功率/DB')pause;figure(3)plot(f,plot_Pxx2);xlabel('频率')ylabel('功率/DB')图9矩形窗图10海明窗图11blackman窗Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。不同窗函数的Welch谱估计在选择窗函数时[10]，一般有如下要求：(1)窗口宽度M要远小于样本序列长度N，以排除不可靠的自相关值；(2)当平稳信号为实过程时，为保证平滑周期图和真是功率谱也是实偶函数，平滑窗函数必须是实偶对称的；(3)平滑窗函数应当在m=0出游峰值，并且m随绝对值增加而单调下降，使可靠的自相关值有较大的权值；(4)功率谱是频率的非负函数且周期图是非负的，因而要求窗函数的Fourier变换是非负的。3.4各种方法的比较与分析(1)在采样点相同的时，周期图法和BT法的特点是离散性大，曲线粗糙，方差较大，但是分辨率较高；(2)Bartlett法和Welch法的收敛性较好，曲线平滑，方差较小，但是功率谱主瓣较宽，分辨率低，这是由于对随机序列加窗截断所引起的Gibbs效应造成的；(3)由仿真结果可以看出：使用不同的窗函数谱估计的质量是不一样的，从图9、10、11可以观察到矩形窗的主瓣较窄，分辨率较好，但方差较大，噪声水平较高；而Blackman窗和Hamming窗的主瓣较宽，分辨率较低，但方差较小，噪声水平较低。因此，在进行谱分析时选择何种窗函数，要视具体情况而定。如果强调高分辨率，能精确读出主瓣频率，而不关心幅度的精度，例如测量震动物体的自震频率，可以选用主瓣宽度比较窄的矩形窗；对受到强干扰的窄带信号若干扰靠近信号，则可选用旁瓣幅度较小的窗函数，若离开通带较远，则可选用渐近线衰减速度比较快的窗函数。总之，要针对不同的信号和不同的处理目的来选择合适的窗函数，这样才能得到良好的效果。所以，与Bartlett法相比，Welch法的估计曲线比较粗糙，但是分辨率较好。原因是Welch法中对数据进行截断时加的是Hanning窗，相对于矩形窗，Hanning窗的主瓣包含更多的能量，因而使功率谱的主瓣较窄，分辨率较高。由上述理论分析及仿真实验可知，经过对几种方法的比较我们可以看出Welch法采用加窗交叠求功率谱，可以有效减小方差和偏差，一般情况下能接近一致估计的要求，因而得到广泛应用。同时还可以发现，对信号加不同的窗函数，谱估计的质量是不同的。4应用举例现取具体信号加随机信号对随机振动信号进行功率谱估计，取周期信号和随机信号叠加，即:随机信号，用经典Welch估计法进行谱估计,程序如下：clear; Fs=1000;n=0:1/Fs:1;xn=sin(2*pi*80*n)+2*sin(2*pi*140*n)+randn(size(n));nfft=1024;window=boxcar(100);%矩形窗window1=hamming(100);%海明窗window2=blackman(100);%blackman窗noverlap=20;%数据无重叠range='half';%频率间隔为[0Fs/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(Pxx1);plot_Pxx2=10*log10(Pxx2);figure(1)plot(f,plot_Pxx);xlabel('频率')ylabel('功率/DB')pause;figure(2)plot(f,plot_Pxx1);xlabel('频率')ylabel('功率/DB')pause;figure(3)plot(f,plot_Pxx2);xlabel('频率')ylabel('功率/DB')功率谱图如图12：图12随机振动功率谱图5结论1、各种方法优缺点的总结由前一章的仿真结果和分析可得:（1）在采样点相同的时，周期图法和BT法有离散性大，曲线粗糙，方差较大的缺点；优点是分辨率较高；(2)Bartlett法和Welch法的优点是收敛性较好，曲线平滑，方差较小，缺点是功率谱主瓣较宽，分辨率低；(3)对于Welch方法的三种窗函数，矩形窗的优点是：主瓣较窄，分辨率较好，缺点是方差较大，噪声水平较高；而Blackman窗和Hamming窗的缺点是主瓣较宽，分辨率较低；优点是方差较小，噪声水平较低。因此，在进行谱分析时选择何种窗函数，要视具体情况而定。如果强调高分辨率，能精确读出主瓣频率，而不关心幅度的精度，例如测量震动物体的自震频率，可以选用主瓣宽度比较窄的矩形窗；对受到强干扰的窄带信号若干扰靠近信号，则可选用旁瓣幅度较小的窗函数，若离开通带较远，则可选用渐近线衰减速度比较快的窗函数。总之，要针对不同的信号和不同的处理目的来选择合适的窗函数，这样才能得到良好的效果。通过对功率谱估计概念的理解和对经典谱估计的几种方法的认识以及仿真