《2022年第六章功率谱估计电子教案》由会员分享，可在线阅读，更多相关《2022年第六章功率谱估计电子教案（24页珍藏版）》请在金锄头文库上搜索。
1、名师精编精品教案第六章功率谱估计的经典方法 6.1 引言从第二章的讨论中， 我们已经知道一个随机信号在各时间点上的值是不能先验确定的，它的每个实现 (样本 )往往是不同的，因此无法象确定信号那样可以用数学表达式或图表精确地表示它，而只能用它的各种统计平均量来表征它。其中，自相关量作为时移的函数是最能较完整地表征它的特定统计平均量值。而一个随机信号的功率谱密度(函数)，则是自相关函数的傅立叶变换。对于一个随机信号来讲，其能量通常为无限大，它本身的傅立叶变换是不存在的，常常需要研究其功率在频域上的分布。因此，功率谱密度是一个随机信号的一种最重要的表征形式。我们要在统计意义下了解一个随机信号，就要求
2、知道 (或估计 )它的功率谱密度。如果我们用()xxRm 表示随机信号)(nx的自相关函数，)(xxP表示它的功率谱密度 (以下简写成 PSD)，则有：()()jmxxxxmPRm e(6.1) 而其中()( )()xxRmE x n xnm(6.2) 即为滞后积的数学期望。根据各态历经假设， 零均值广义平稳随机过程的集合的平均可以用一个样本序列的时间的平均代替，于是上式可写成1()lim( )()21NxxNnNRmx n xnmN(6.3) 实际上，首先不可能获得样本序列的所有数据，即无数个x(n)，其次通常检测到的数据都是含有噪声或者干扰。因此，只能根据有限个含有噪声的检测数据估计随机信
3、号的自相关序列，进而估计功率谱。将式(6.3)代入式 (6.1)得1()lim( )()21Nj mxxNmnNPx n xnm eN1l i m()()21NjmNnNmx nxnmeN令mnl，则mln，上式可写成精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 24 页名师精编精品教案1()lim( )( )21NjnjlxxNnNlPx n ex l eN2)(121limnjNNnNenxN(6.4) 式(6.4)在N的极限情况下是不可能收敛的， 这是因为对于无限时域的随机信号，它的傅氏变换是不存在的。实际上只有将式(6.4)求
4、平均，成为2)(121lim)(njNNnNxxenxNEP(6.5) 才有意义。以后我们还将会看到，只有将式(6.4)经过求平均 (或平滑 )，即只有式 (6.5)才能满足一个正确的估计必须满足的一致估计的条件。由于实际得到的随机信号只能是它的一个实现或一个样本序列的片段，因此问题是如何根据它的有限个样本序列来估计信号的自相关函数或功率谱密度。这是本章要讨论的中心内容。当我们用一个样本的记录的有限个数据)1(,),1 (),0(Nxxx来估计自相关函数和功率谱密度时，有101()( ) ()NxxnRmx n x nmNnNNmnxnxN)()(1(6.6) ?( )( )jmxxxxmPR
5、m e(6.7) 这里( )( ) ( )()() ()NNNNxnwn x nxnmwnm x nm(6.8) wN为矩形函数，1, 00,0NnNwnnN及或按式 (6.4) 2210)(1)(1)(?NnjNnxxXNenxNP(6.9) 这里)(NX是有限长序列)1, 1, 0()(NnnXN的傅氏变换。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 24 页名师精编精品教案一个好的估计应该是无偏估计，最小方差估计。如果我们用表示某个随机变量的真值，?表示它的估计值，则希望满足：(1) 无偏估计无偏估计，即?的偏差 (Bias)
6、为零，所谓偏差 (用 B 表示)定义为?BBiasE无偏估计即 B0，?E的估计。图 6.1 中的估计 1 和估计 2都属于无偏估计。(2) 最小方差估计最小方差估计，即方差2? VarEE为最小的估计。图6.1 中，估计 2 较之估计 1 方差小。图 6.1 二种估计的概率密度分布但是常常发生这种情况，一种估计的偏倚较小，而方差较大；另一种估计偏倚较大而方差较小；此时很难定哪一种估计好。因此也常常用均方误差的大小来衡量估计的优劣。在第二章中我们已经讨论到均方误差定义为22?EeE不难证明2?22BeE均方误差为2eE与偏差和方差均有关，要2eE最小就要求 B2与2?之和最小。由于，当N时式(
7、6.6)就成为式 (6.3)。因此应有xxxxNRR这就是说，当观察到的样本的数据有无限多个时，按照无穷多个这样的样本数据估计到的自相关函数应该就是自相关函数的真值(各态历经假设 )。换句话说，一个正确的估计精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 24 页名师精编精品教案应有0?ar0?VNBiasN(6.10) 满足式 (6.10)的估计称为一致估计。一个正确的估计应该满足一致估计的条件(这是正确估计的必要条件，不是充分条件)。反之，如果某种估计方法不能满足式(6.10)一致估计的条件，则这种估计方法一定是不正确的。下面我们在
8、讨论各种估计方法时，常常以此作为估计正确与否的主要准则之一。功率谱估计有着极其广泛的应用，不仅在认识一个随机信号时，需要估计它的功率谱。它还被广泛地应用于各种信号处理中。下面我们举三个应用的例子。 在信号处理的许多场所， 要求预先知道信号的功率谱密度(或自相关函数 )。 例如，在最佳线性过滤问题中，要设计一个维纳滤波器就首先要求知道(或估计出 )信号与噪声的功率谱密度 (或自相关函数 )。根据信号与噪声的功率谱(或()xxRm )才能设计出能够尽量不失真的重现信号，而把噪声最大限度抑制的维纳滤波器(见第二章 )。 常常利用功率谱估计来得到线性系统的参数估计。例如，当我们要了解某一系统的幅频特性
9、)(H时，可用一白色噪声n 通过该系统。再从该系统的输出样本y(n)估计功率谱密度)(yyP。由于白色噪声的PSD(用)(P表示)为一常数即2)(P，于是有22)()(HPyy故通过估计输出信号的PSD，可以估计出系统的频率特性)(H(模特性 )。在第七章将要讨论到用自回归模型法估计PSD 的一节中，我们将要具体讨论系统参数估计与 PSD 估计间的关系。 从宽带噪声中检测窄带信号。这是功率谱估计在信号处理中的一个重要用途。但是这要求功率谱估计有足够好的频率的分辨率，否则就不一定能够清楚地检测出来。所谓谱估计的分辨率可以粗略地定义为能够分辨出的二个分立的谱分量间的最小频率间隙(距)。谱估计问题无
10、论从认识一个随机信号或从其它应用方面来讲都是重要的。因此对谱估计方法的研究引起了国内外学者的广泛注意与重视。它是当前在信号处理中的一个十分活跃的课题。功率谱估计总的来讲可以分为经典谱估计方法和现代谱估计方法。经典谱估计为线性、非参数化方法，需采用经典的傅里叶变换及窗口截断。经典谱估计方法包括周期图法，相关图法等，对长序列有良好估计。现代谱估计为非线性、参数化方法，包括最大似然估计，最大熵法， AR 模型法，预测滤波器法， ARMA 模型等。对短序列的估计精及精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 24 页名师精编精品教案度高，与
11、经典法相互补充。 现代谱估计是融合经典变换理论、统计估计理论、 系统辨识、信息论、时间序列分析及计算方法等理论与技术的一门新学科。目前应用广泛，发展迅速。功率谱估计的经典方法以傅立叶变换为基础，很早就被提出， 在 FFT 算法出现后得到了广泛的应用。本章内容为经典谱估计方法。 6.2 自相关序列的估计功率谱的估计要求计算自相关序列，下面讨论自相关序列的估计方法。 6.2.1 自相关序列的无偏估计设观察到 N 个样本序列nx的值：1,),1(),0(Nxxx。现在要由此N 个数据来估计自相关函数()xxRm 。由于nx 只能观察到10Nn的 N 个值，而10Nnn与时的nx 值是不知道的，那么滞
12、后积序列*( )( )()mynx nx nm 是一个长为N|m|的序列，因此，式 (6.6)成为10101()11( )NmxxnnmnNmmnRmx xmNNmynNm(6.11) 式中 m 取绝对值是因为()()xxxxRmRm ，m 为负值时上式仍适用。式 (6.11)规定的求和上下限的原则是保持充分利用全部(N 个)数据。这种估计方法的效果如何，我们首先需要看它的偏差与方差是否满足一致估计的条件。由式(6.11)，得110011()( ) ()()NmNmxxxxnnE RmE x n x nmRmNmNm() 1xsRmmN上面 Rxs(m)是自相关函数的真实值。所以()()( )
13、0 1xxxxxxBias RmRmE RmmN故这种估计，当1Nm时，属于无偏估计。现在来求()xxVarRm，按定义22()()()xxxxxxVarRmE RmERm精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 24 页名师精编精品教案22()()1xxxxERmRmmN(6.12) 又按式 (6.11) 122001()( ) () ( ) ()()NmxxnkE RmE x n x nm x k x kmNm(6.13) 当随机序列nx是零均值，实、高斯序列时，有3241423143214321xxExxExxExxExxE
14、xxExxxxE所以)()()()(mkxkxmnxnxE)()()()()()()()(mkxmnxEkxnxEmkxkxEmnxnxE)()()()(kxmnxEmkxnxE22()()()()xxxxxxxxRmRknRkmn Rkmn代入式 (6.13)，得12222001( )()()()()()NmxxxxxxxxxxnkE RmRmRknRkmn RkmnNm代入式 (6.12)，得122001()()()()()NmxxxxxxxxnkVar RmRknRkmn RkmnNm令nkr，显然 r 的最小值为1mN，最大值为1mN，且0r)(nk即的情况将出现mN次，1r的情况将出
15、现1mN次以此类推，对不同 r 值的情况，出现的次数将为rmN，于是上式可写成(1)22(1)1()()()()()rNmxxxxxxxxrNmVar RmNmrRRrm RrmNm122(1 )11()()()()NmxxxxxxrNmmrRrRrm RrmNmN122(1 )1()()()()NmxxxxxxrNmRrRrm RrmNm(6.14) 当 Nm 时，上式以 1/ N 趋于零，即lim()0xxNVar Rm。故( )xxRm 的方差满足一致估计的条件。如果nx不是高斯过程，在上式中需要再加一项，但此项往往是可以忽精选学习资料 - - - - - - - - - 名师归纳总结
16、- - - - - - -第 6 页，共 24 页名师精编精品教案略的，因此，式 (6.14)仍近似适用。式(6.11)这种估计自相关函数的方法，虽然当m很小于 N 时能得到一致估计，但当N 一定，|m|接近于 N 时，即Nm时，()xxVar Rm就变得很大，因而不能得到有用的估计。 6.2.1 自相关序列的有偏估计因此很多学者如 Jenkins-Watts和 Parzen等都主张按下式估计()xxRm ， 我们用()xxRm表示101()( ) ()NmxxnRmx n x nmN() ,(1 )xxNmRmmNN(6.15) 同时()()xxxxNmE RmRmN这相当于它的均值等于真值
17、()xxRm 用三角窗函数加权。故()xxRm 是有偏的，其偏差为()()()xxxxxxBias RmRmE Rm()()xxxxNmRmRmN()xxmRmN因此，()xxRm 是()xxRm 的渐进无偏估计。同时2()()xxxxNmVar RmVarRmN（()xxVar Rm）122(1 )11()()()NmxxxxxxrNmmrRrRrm RrmNN事实上，将xxR 用三角窗函数加权后，不仅使方差减小，而且有利于钝化()xxRm 的截断边界，从而改进对 P的估计。讨论：虽然()xxRm 的 Bias和 Var都不等于零， 但当N时0,0 VarBias，说明()xxRm 是()x
18、xRm 的渐进无偏估计和有效估计，因此，()xxRm 是( )xxRm 的一致估计。且精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 24 页名师精编精品教案()( )xxxxVar RmVar Rm，同时可以证明( )xxRm 的均方误差小于()xxRm的，以及可改进对 P的估计，所以今后我们还是用( )xxRm 作为自相关函数的估计，即101()( ) ()NmxxnRmx n x nmN(6.16) 此外， m 越大，分辨率越高，但是自相关序列估计的偏差也相应增大，通常取10 5mNN，效果较好。通过将自相关函数的估计进行傅氏变换
19、求得功率谱估计的方法即为BT PSD 法。 6.3 周期图作为功率谱的估计功率谱估计的经典法实质上就是传统的傅立叶分析法，它又可分成二种。一种是间接法，它先通过式 (6.6)对自相关函数进行估计 (一般都需要窗函数将自相关值加权，以减小自相关序列截断的影响)，然后再通过式 (6.7)作傅立叶变换得功率谱估计值。 这种方法是 1958 年由 Blackman 与 Tukey 提出的，简称 BT PSD 估计法，也称为相关图法。另一种是直接法，它是将观察到的有限个样本数据)1(,),1(),0(Nxxx利用 FFT 算法作傅立叶变换直接按式 (6.9)进行功率谱估计 (不通过自相关函数的估计)，这
20、种方法称为周期图法。本章的经典法中主要讨论周期图法。用周期图 (包括平滑后的周期图 )作为功率谱估计的方法可利用FFT 进行计算，因而有计算效率高的优点，在谱分辨率要求不高的地方常用这种周期图法进行谱估计。它的一个主要缺点是频率分辨率低。 这是由于周期图法在计算中把观察到的有限长的N 个数据以外的数据认为是零。 这显然与事实不符合。 把观察不到的值认为是零， 相当于将 x(n)在时域里乘上了一个矩形窗口函数，在频域里相当于与一个sinc 函数卷积，由于sinc函数与函数比较有二方面的差别，其一是其主瓣不是无限窄，其二是有旁瓣，因此卷积的结果必然造成失真。由于主瓣不是无限窄，如果原来真实的功率谱
21、是窄的，与主瓣卷积后将使功率向附近频域扩散，使信号模糊，降低了分辨率，主瓣愈宽分辨率愈差。由于存在旁瓣，又将产生两个后果，一是PSD 主瓣内的能量“泄漏”到旁瓣将使谱估计的方差增大，二是与旁瓣卷积得到的信号功率谱完全属于干扰。在严重的情况下，强信号与旁瓣的卷积可以大于弱信号与主瓣的卷积，使弱信号淹没在干扰的强信号中，而无法检测出来。这正是周期图作为功率谱估计的二个主要缺点。对于 BT PSD 估计法，由于它是按式 (6.7)将功率谱用有限个自相关函数值的傅氏变换代表按式 (6.1)用无限个自相关函数值的傅氏变换求得的真实功率谱。这相当于将序列()sxRm 乘了一个矩形窗函数，因此也同样存在上述
22、二个缺点。许多学者想尽办法，企图选择适当窗口函数的形状来提高经典法的谱分辨率，但是发现，所有能降低旁瓣的窗口函数都是以主瓣的增宽为代价的；反之亦然。这两个缺点只能互换而无法同时改善。因精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 24 页名师精编精品教案此，近几年来， 在提高功率谱估计的分辨率方面提出了很多新的方法，其中包括以 1967年 Burg 提出的最大熵谱分析法为代表的现代谱估计法，现代谱估计将在下一章阐述。 6.3.1 周期图将式(6.16)代入(6.7)得?()()jmxxxxmPRm e101()()Nmjmmnx n
23、x nm eN1()()jmNNmnxn xnm eN这里)()(mnxnxNN与的下标 N 表示它们是有限长 (长为 N)的序列，令 l=n+m，有1?( )( )( )jnjlxxNNnlPxn exl eN2)(1)()(1XNeXeXNjNjN其中的)(NX为nNnnjnjNjNNenxenxeXX10)()()()(即)(NX是有限长序列x(n)的傅氏变换。显然)(NX是周期性的。直接将)(NX的模的平方除以 N 求得的功率谱的估计称为周期图，并用)(NI表示。于是有2)(1)()(?NNxxXNIP(6.17) 如果我们观察到x(n)的 N 个值：)1(,),1 (),0(Nxxx
24、，可以通过FFT 直接求得)(KXN频率离散化的)(jeX。 然后按照式 (6.17)直接求得)(xxP(不必先通过估计自相关函数 )。这种将周期图作为功率谱估计的方法的主要优点是计算方便，它可以直接用FFT 算法从 x(n)得到)(NX。从而得到)(xxP。正是由于它的这个优点，使这种方法成为一种十分通用的方法。 6.3.2 周期图的估计性能为了了解周期图作为功率谱估计的估计效果，让我们来讨论它的偏差和方差。为此首先需要求周期图的期望值。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 24 页名师精编精品教案按式(6.16)，得1()
25、( )()xxNNnE RmE xn xnmN1( ) ( )() ()NNnE wn x n wnm x nmN1( )()( ) ()NNnwn wnmE x n x nmN(6.18) 这里 wN代表矩形序列。令11()( )()()()NNnmwn wnmw mwmNN(6.19) 由于)(m是二个矩形函数的卷积，因此它必是一个三角函数，常常称它为Bartlett 窗函数，用)(mB表示，不难证明0,1)(其它，NmNmmB而它的傅氏变换为222/sin2/sin1)(1)(NNRNWNB(6.20.0) 图 6.1 所示为 Bartlett 函数及其傅立叶变换示意图。图 6.2 Ba
26、rtlett 函数及其傅立叶变换示意图又( ) ()()xxE x n x nmRm自相关函数的真值将上式与式 (6.19)代入式 (6.18)，并求其傅氏变换，得精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 24 页名师精编精品教案1()()()2NBxxE IWPdPWxxB)()(21(6.20) 由式(6.20)可见，除非)(BW为函数，)(NIE将不等于)(xxP，故周期图作为)(xP的估计是有偏的：0)()()(NxxNIEPIBias当N时NBNmm11)(故)()(NBW)()(limxxNNPIE因此，周期图作为功
27、率谱的估计，当N时是无偏的。( )perS是渐进无偏估计。其次求周期图的方差。为了得到周期图的方差，首先假设序列)10()(Nnnx是一个实，白色，零均值过程的样本，具有高斯概率分布函数。按方差的定义应有)()()(22NNNIEIEIVar(6.21) 按式(6.17)，周期图可表示为2)(1)(jNNeXNI1()()jnjkNNnkxn exk eN()1( )( )jn lNNnxn xk eN将( )( ) ( )NNxnwn x n 和( )( ) ( )NNxkwk x k代入上式，并求期望值：()1()( )( )( ) ( )jn kNNNnkE Iwn wk E x n x
28、 keN(6.22) 其中( ) ( )()xxE x n x kRnk精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 24 页名师精编精品教案由于我们假设)(nx是白色的，所以2()()xxxRnknk 代入上式得22221()( )( )xNNxNnnE IwnwnNN22( )xNnwnN为了求得)(NIVar，按式(6.21)除了需要求得)(NIE以外，还需要求得)(2NIE，为此我们先求在频率为的及21)(NI的协方差，参照式 (6.22)可得12()()NNE II()()1221( )( )( )( )( ) ( ) (
29、) ( )jjn kpqNNNNnkpqnk wp wqE x n x k x p x qeeN(6.23) 对于白色高斯过程可以证明)()()()(qxpxkxnxE( ) ( )( ) ( )( ) ( )( ) ( )( ) ( )( ) ( )E x n x kE x p x qE x n x pE x k x qE x n x qE x k x p矩分解定理考虑到对于白色过程有2()( ) ()()xxxRmE x n x nmm 所以其它，及或及或及当，0)()()()(4pkqnqkpnqpknqxpxkxnxEx代入式 (6.23)得4()()222212122()()( )(
30、 )( )( )j n kxNNNNNNppnkE IIwnwpwn wk eN() ()2221( )( )j n kNNnkwn wk e(6.24) 上式当 n=k 及 p=q 时只有第一项存在 (其它二项为 0)； 当 n=p 及 k=q 时只有第二项存在；当 n=q，k=p 时，只有第三项存在。又因为22222( )( )()NNNnkpwnwkwpN故式(6.24)成为精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 24 页名师精编精品教案41111()()()()212121220000()()NNNNj nkj n k
31、xNNnknkE IINeeN(6.25) 而其中1010)()(1010)(212121NnNkjkjnNnNkknjeee)()()()(212121211111jjNjjNeeee2/ )(2/ )(2/ )(2/)(2/ )(2/ )(2/)(2/ )(2/ )(2/ )(2/ )()(212121212121212121212121jjjjNjNjNjjjjNjNjNeeeeeeeeeeee22/)(2/ )(2/ )(2/ )(21212121jjjNjNeeee21212sin()2sin() 2N同理可得211()()12120012sin()2sin() 2NNj nknkN
32、e代入式 (6.25)得2212122121224212/)sin(2)(sin2/)sin(2)(sin)()(NNNNIIExNN221212212142/ )sin()(sin2/ )sin()(sin1NNNNx当21时得1sinsin1)(242NNIEXN因为上式方括号中只有一项不为0，故有)()()(22NNNIEIEIVar4x由公式可见，当N，0)(4xNIVar。这说明周期图不满足一致估计的条件。不论 N 取多长，)(NIVar都有4x的量级。因此周期图不是对功率谱的好的估计。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13
33、 页，共 24 页名师精编精品教案实际上当N时式(6.17)成为式 (6.4)。正像我们在讨论式 (6.4)时所指出的，对于无限能量的随机序列， 它的傅氏变换是不存在的， 因此式 (6.4)在N的极限情况下是不可能收敛的。 故也就不能期望当N时)(NI会等于它的真值)(xxP，而满足一致估计的条件。为了使周期图作为功率谱估计满足一致估计的条件，我们必须将周期图进行平滑 (或平均 )处理。如果我们需要求周期图的协方差，则有)()()(),()(),(212121NNNNNNIEIEIIEIICOV4221212212142/ )sin(2/)(sin2/ )sin(2/ )(sin1xxNNNN
34、221212212142/ )sin(2/)(sin2/ )sin(2/)(sinNNNNx如果Nk /21以及Nl /22(k 和 l 为整数 )则上式成为22421/ )(sin)(sin/)(sin)(sin)(),(NlkNlkNlkNlkIICOVxNN上式在lk以及lklk与不是 N 的整倍数时等于零。 因此以N/2的整倍数为频率间距的周期图的值是不相关的，当N 增大时，协方差为零的功率谱样本之间的间距减小。因此可以预料到 N 增加将会使周期图的起伏加快。图6.3 说明了周期图的这种特性。图(a)、(b)、(c)分别表示 N 为 16，32 和 64 时的高斯白噪声序列的周期图。图
35、 6.3 高斯白噪声序列的周期图精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 14 页，共 24 页名师精编精品教案由于22sin)sin(1)()(NNPIVarxxN故有N)()(2xxNPIVar虽然本节推导的结果是以假设高斯概率密度为根据的，但其定性结果在一个相当宽广的范围内成立。总之，周期图法不是功率谱的一致估计， 而且周期图都在真实功率谱附近随机起伏。这主要是由于时域序列加窗后， 频域发生泄漏现象， 包括主瓣平滑作用和旁瓣泄漏作用。所以人们作了许多改进，加各种形式的窗，还有分段求估计再求平均的平均周期图法，平滑周期图法等。然而，由于将
36、测得数据以外的数据均看作零（实质相当于乘上了一个矩形窗口），因此频率分辨率低，频谱能量向旁瓣泄漏，且不是功率谱的一致估计。即使进行改进，也不能从根本上改善周期图的性能，只适用于数据较长和频率分辨率要求不高的场合同时，如果增加信号的数据点数，数据量增大，相当于滞后窗加宽，相应的傅立叶变换的主瓣变窄，即功率扩散的频率范围变窄。同时，方差不会降低，随机起伏现象仍然存在。且数据点越多，随机起伏越密集。经典谱估计的优点是计算简单。 6.4 平滑后的周期图作为PSD 的估计上面我们已证明周期图作为功率谱的估计不满足一致估计的条件。因为当N，方差0)(NIVar，因此必须作一些修正。本节将表明如果将周期图进
37、行平滑(平均是一种主要的平滑方法)将会使方差减小，得到一致的谱估计。 平滑的主要方法有二种。 在 FFT 出现并广泛应用以前主要用窗口处理法进行平滑，即选择适当的窗口函数作为加权函数进行加权平均来加快收敛速度。另一种平滑方法是平均周期图的方法，即先将数据分段，再求各段周期图的平均值。后一种方法又称为 Bartlett 方法，是当前用得最多的一种平滑方法。Welch 又对 Bartlett 方法进行了改进并提出了用FFT 计算的具体方法。6.4.1 巴特利特 (Bartlett)平均周期图的方法首先让我们来看一下为什么周期图经过某种平均(或平滑 )后会使它的方差当N时趋于零，达到一致估计的目的。
38、如果Lxxx,21是不相关的随机变量，每一个具有期望值，方差2，则可以证明它们的数学平均LxxxxL/ )(21的期望值等于，数学平均的方差等于L/2，现证明如下：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 15 页，共 24 页名师精编精品教案LLxxxELxEL1121222)()(xExExExExVar22212)(1LxxxEL2112222121LjLjiijiLxxExxxELLjLjiiijLjjiijixExExxE1111222)1(LLLL所以2122222122211LxELLLxELxVarLiiLii22222221
39、212)()()(1LLxExExExExExELLLL2221(6.26) 由(6.26)可见，L 个平均的方差比每个随机变量的单独方差小L 倍。 当0xVarL，可达到一致谱估计的目的。 因而降低估计量的方差的一种有效方法是将若干个独立估计值进行平均。把这种方法应用于谱估计应归功于Bartlett。Bartlett 平均周期图的方法是将序列)10()(Nnnx分段求周期图再平均。设将)(nx分成 L 段，每段有 M 个样本，因而LMN，第 i 段样本序列可写成LiMnMiMnxnxi1, 10)()(第 i 段的周期图为210)(1)(MnnjjiMenxMI假定各段的周期图)(iMI是互
40、相独立的，则谱估计可定义为L 段周期图的平均，即LiiMxxILP1)(1)(?(6.27) 于是它的期望值为精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 16 页，共 24 页名师精编精品教案LiiMiMxxIEIELPE1)()(1)(?将式(6.20.0)与式(6.20)代入上式得deWPPEjBxxxx)()(21)(?)(dMPMxx22/sin2sin)(21(6.28) 这里LNM/，因此 Bartlett 估计的期望值是真实谱)(xxP与三角窗函数的卷积。由于三角窗函数不等于函数，所以 Bartlett 估计也是有偏估计即0Bia
41、s，但当N时，0Bias。由于我们假定各段周期图是相互独立的，所以可按式(6.26)得到下式：)(1)(?MxxIVarLPVar22s i n)si n (1)(1MMPLxx(6.29) 由此可见，随着 L 的增加)(?xxPVar是下降的， 当L时，0)(?xxPVar。因此 Bartlett估计是一致估计。 M 和 N如何选择？比较式 (6.28)的)(?xxPE与式 (6.20)的)(NIE可见在二种情况的估计量的期望值都是真值)(xxP与窗口函数)(jBeW的卷积形式，但后者将前者WB中 N 改为M，NLNM/。因而使)(jBeW主瓣的宽度增大。 由于主瓣的宽度愈窄， 愈接近函数，
42、偏 差 愈 小 。 式 (6.28)中)(jBeW的 主 瓣 宽 度 大 于 式 (6.20) 中 的 主 瓣 宽 度 ， 因 而)()(?NxxIBiasPBias，而主瓣愈宽显然使分辨率愈差。因此Bias可用来说明谱的分辨率，Bias 愈大说明谱分辨率愈差。一个固定的记录长度N，周期图分段的数目L 愈大将使方差愈小，但 M 也愈小，因而使 Bias愈大，谱分辨率变得愈差。因此Bartlett 方法中 Bias 或谱分辨率和估计量的方差间是有相互制约关系的。M 和 N 的选择一般是要根据所研究的信号由用户选定。例如，如果我们知道谱有一个窄峰，同时如果分辨出这个峰是重要的，那么我们必须选择M
43、足够大。又从方差的表达式我们可以确定谱估计的可接受的方差所要求的记录长度)(LMN。由此可见 Bartlett 法使谱估计的方差减小是以增加Bias以及降低谱分辨率为代价的。实际上，当 N 一定时， Bias与 Var的相互制约性是谱估计的一个固有特性。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 17 页，共 24 页名师精编精品教案例如 为了说明经平均后的周期图作为功率谱估计的实际效果，设有一零均值高斯分布的随机过程，其功率谱密度为1.01)1)(1(1)(1aazazzPxx这一功率谱密度是由一零值高斯分布单位方差的噪声序列通过一个其)1/
44、(1)(1azzH的滤波器形成的。 为了简便，设选用8N的矩形窗函数。 其周期图的期望值 (用)(8IE表示)与真值)(xxP均表示在图 6.4 中。它说明8N的周期图可以得到的Bias的情况。图6.5 表示8M分 4 段与 16 段二种情况平均后的周期图。 显然 L=4 的方差比 L=16的大。将 L=16 的曲线与图 6.4 的曲线比较可见在这种估计中误差的大部分起因于Bias 而不是方差(因为二种情况均为8M，Bias 相仿，误差也相仿 )。M=16，L=2 及 8 的周期图表示在图 6.6。图 6.6 中 L 不同造成的影响也是明显的。但是这二个曲线的起伏都很大，因此有理由认为误差主要
45、起因于方差。比较8M与 M=16 的周期图可见，在 6.3 中得到的 M 愈大(即 6.3 中的 N 愈大)， 将使周期图的起伏愈增快的结论， 在这里也同样成立。比较图 6.5 与图 6.6 发现在这个例子中最好的选择是应用L=16，8M的估计而不是 L=8、M=16 的估计，即宁可减小方差，牺牲Bias。在实际中，当然功率谱密度的真值是不知道的。但是谱的窗口函数以及关于功率谱密度的某些信息往往是预先知道的。通过改变 M 和 L 以及利用预先知道的情况，通常可以很好地进行选择。平均周期图的方法特别适合于应用FFT 算法。 因此在 FFT 出现以后这个方法比下面将要讨论的利用窗口函数的处理法用得
46、更多。而在 FFT 出现以前主要是用窗口函数处理法平滑周期图。图 6.4 )(xP与)(8IE的特性图 6.5 平滑后的周期图 （每段取 8 个数据）精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 18 页，共 24 页名师精编精品教案图 6.6 平均后的周期图 (每段取 16 个数据 )6.4.2 加窗平滑周期图滞后窗的加窗效应源自无限长序列取样时加的数据窗。周期图的期望值被主瓣平滑的程度和功率向旁瓣泄露的多少，取决于加在信号x 上的数据窗 w(n)的类型。当 w(n)为矩形窗，其傅里叶变换的主瓣比较窄，故周期图被平滑的程度最小，但是由于矩形窗的
47、傅里叶变换的旁瓣比其他窗的旁瓣要高，因而旁瓣泄漏也最严重。如果可以选择合适的窗， 比如海明窗，以减弱旁瓣泄漏， 使弱窄带信号也显现出来，从而提高频率分辨率，但是旁瓣的下降是以主瓣的变宽为代价的。因此，平滑周期图的另一种常用的方法，下面讨论加窗平滑后周期图的偏差。平滑周期图，实际上就是用一个合适的窗口谱函数)(jeW与周期图卷积，即deWIPjNxx)()(21)(?)(6.30) 或1(1)?( )()()Mj mxxxxmMPRmm e(6.31) 这里?()xxRm与)(m分别是)(NI与)(jeW的傅氏反变换， 并设序列)(m长 2M-1。在第一章中我们已讲到)(?xxP是一个实，偶，非
48、负函数，)(m应是一个偶序列，并且满足条件：当0)(jeW(6.32) 例如前面用到的三角或Bartlett 窗函数是满足此条件的，而海宁(Hanning)与汉明(Hamming)窗函数就不满足式 (6.32)的条件虽然这后二种窗函数能够提供较好的频率分辨率以及较低的旁瓣，但会产生负的功率谱估计。如果我们将时域与频域对换， 由式(6.30)的卷积形式，可以看成)(?xxP是)(NI通过一个具有“单位样本响应”)(jeW的滤波器产生的。)(jeW的“频率响应”)(m具有低通特性，因此可以平滑)(NI，对)(NI中快变化成分有滤除作用， 因此我们称)(?xxP为平滑后的周期图。精选学习资料 - -
49、 - - - - - - - 名师归纳总结 - - - - - - -第 19 页，共 24 页名师精编精品教案因式(6.30)可得deWIEPEjNxx)()(21)(?)(又按式 (6.20)，有deWPIEjBxxN)()(21)()(因此，有)()()()(?jjBxxxxeWeWPPE1(1)()()()MjmxxBmMRmmm e这里NmNmmB,1)(现在如果 M 比 N 小得多，则1,BNm，于是式 (5.48)可近似写成deWPPEjxxxx)()(21)(?)(6.33) 将式(6.33)与式(6.20)比较可见，如将式 (6.20)右边的)()( jBeW换成)()( j
50、eW，则式 (6.20)就全同于式 (6.33)。事实上，由于周期图法是将观察到的x(n)的 N 个以外的值均认为是零，这实际上已把 x(n)乘上了一个矩形窗口 (即使我们并没有作窗口处理)。式(6.20)中的三角窗函数)(jBeW正是这个矩形窗口造成的，加窗口处理后将式(6.20)改变成式 (6.33)。等于将式 (6.20)中的窗函数)(jBeW改变成)(jeW。如果)(jeW的主瓣宽度大于)(jBeW的主瓣宽度，则可进一步平滑谱估计，减小谱估计的方差，但是同样带来 Bias的增加和谱分辨率的变差。 至于选用怎样的窗函数才能有较少的 Bais 的增加与分辨率的降低来换取较多的谱估计的方差的
51、改善以及窗口处理对谱估计方差的影响的具体研究，这里不再讨论。常用窗函数例如：汉宁窗函数：2( )cos (/)0.50.5cos(2/)/ 2w nnNnNnN海明窗函数：2( )0.50.5cos(2/)sin (/)/ 2w nnNnNnN精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 20 页，共 24 页名师精编精品教案例题 1：1212( )sin()cos()( ),1,0.05 2 ,0.12 2 ,(x nnnw nnNw n ） 为高斯白噪声，用平均周期图法和Bartlett 法求信号的功率谱密度估计。N=1024; Ns=256
52、; pxxo=10*log10(abs(fft(xn).2)/(N-1);%周期图无重叠分段平均法pxx1=abs(fft(xn(1:256),Ns).2)/Ns; pxx2=abs(fft(xn(257:512),Ns).2)/Ns; pxx3=abs(fft(xn(513:768),Ns).2)/Ns; pxx4=abs(fft(xn(769:1024),Ns).2)/Ns; pxxb=10*log10(pxx1+pxx2+pxx3+pxx4)/4); 结果分析：从频谱看出， 与周期图谱估计比较， 无重叠分段平均周期图法 Bartlett 法估计曲线均较为平坦，但是波峰加宽，频率分辨率下降
53、。6.4.3 Welch法Welch 提出对 Bartlett 法的修正，使其更适合于用FFT 进行计算。他主要提出二方面的修正，一方面是在分段时，可使各子序列数据段之间有重迭，这样将会使方差减小(当 N 与 M 一定时 )。他建议以重迭5070%为宜；另一方面是是选择适当的窗函数)(n，并在周期图计算前直接加进每个子序列。设子序列为 xi(n) 长度为 L,相邻两个子序列有LD 点重叠，于是第i 个子序列为：( )(),0,1,1;0,1,1ix nx niDnLiKK 个子序列总长度为N=L+D( K-1)。若 D=L，则 N=KL，即为 Bartlett 法。这样得到的每一段的周期图为k
54、iennxMUINnnjiiM)()(1)(210)(6.34) 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 21 页，共 24 页名师精编精品教案这里102)(1NnnWMU为归一化因子，而 Bartlett 法每段的周期图为210)()(1)(NnnjiiMenxMI(6.35) 这样， Welch 周期图表示式为：211101()() ()KLjjnWinSeN x nid eKLU或1( )01()( )KjiwMiS eIK(6.36)这样加窗函数的优点是无论什么样的窗函数均可使谱估计非负。例题 2: 对例题 1 的信号分别分析分段加
55、窗平滑周期图和重叠分段加窗平滑周期图。pxxo=10*log10(abs(fft(xn).2)/(N-1);%周期图无重叠分段加窗平均周期图w=hanning(256); pxx1=abs(fft(w.*xn(1:256),Ns).2)/norm(w)2; pxx2=abs(fft(w.*xn(257:512),Ns).2)/norm(w)2; pxx3=abs(fft(w.*xn(513:768),Ns).2)/norm(w)2; pxx4=abs(fft(w.*xn(769:1024),Ns).2)/norm(w)2; pxxp=10*log10(pxx1+pxx2+pxx3+pxx4)/
56、4); 重叠分段加窗平均周期图法，即Welch 法pxx1=abs(fft(w.*xn(1:256),Ns).2)/Ns; pxx2=abs(fft(w.*xn(129:384),Ns).2)/norm(w)2; pxx3=abs(fft(w.*xn(257:512),Ns).2)/norm(w)2; pxx4=abs(fft(w.*xn(385:640),Ns).2)/norm(w)2; pxx5=abs(fft(w.*xn(513:768),Ns).2)/norm(w)2; pxx6=abs(fft(w.*xn(641:896),Ns).2)/norm(w)2; pxx7=abs(fft(
57、w.*xn(769:1024),Ns).2)/norm(w)2; pxxw=10*log10(pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7); 频谱图为：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 22 页，共 24 页名师精编精品教案结果分析： 从上述例子看到： 与 Bartlett 法相比，Welch法与它的周期图谱估计频率分辨率近似，但是方差减小了。说明：对于给定数据量N，增加子序列重叠点数，将增加子序列之间相关性，加大运算量。一般情况下， Welch法以重叠 50%70%为宜。周期图性能评价：描述谱估计技术
58、的性能参数有两个：1、变异性：2()=()jxxjxxVar SeESe2、品质因数：=，它是变异性与分辨率之积。经典谱估计方法的每种技术的品质因数都是近似的，都与记录数据长度成反比，但每种方法的分辨率和方差都不相同（见下页表格数据）。然而，它们的总的性能指标受到数据记录长度的限制。所以，经典谱估计以傅立叶变换为基础，具有计算效率高的优点，但是由于将未观测数据认为 0 和数据加窗，而具有频率分辨率低、旁瓣泄漏等严重的缺陷。下表为三种经典谱估计方法的性能参数比较：变异性分辨率（w）品质因数周期图1 0.89(2/)N0.89(2/)NBartlett 法1/K 0.89(2/)KN0.89(2/
59、)NWelch 法9/8K 1.28(2/)N0.72(2/)N复习题：1、功率谱如何定义？谱估计的任务是什么？经典谱估计和现代谱估计的主要区别是什么？经典谱估计方法有哪些优缺点？精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 23 页，共 24 页名师精编精品教案2、如何估计随机信号的自相关序列？自相关序列的无偏估计和渐近无偏估计有何区别？3、什么是周期图？什么是计算周期图的直接方法和间接方法？4、滞后窗对周期图的期望有何影响？什么是谱估计的频率分辨率和旁瓣泄漏？5、周期图的偏差和方差如何计算？周期图的分辨率如何估算？6、改善周期图质量有哪些方法？请简要说明。7、Welch 周期图法和 Bartlett 法周期图法的主要区别是什么？8、比较各种周期图的性能。9、在使用 Welch 周期图法时，如何考虑序列分段长度和相邻分段间的重叠部分的长度的选取？精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 24 页，共 24 页
