经典功率谱估计课件

数字信号处理电气信息工程学院蔡超峰数字信号处理电气信息工程学院引言对各态遍历随机信号X(n)，自相关函数和功率谱密度均可用时间平均来定义：

维纳-辛钦定理：引言对各态遍历随机信号X(n)，自相关函数和功率谱密度均可第十三章经典功率谱估计周期图法（直接法）间接法直接法和间接法的关系直接法和间接法估计的质量、直接法的改进经典功率谱估计总结短时傅里叶变换第十三章经典功率谱估计周期图法（直接法）1.周期图法（直接法）周期图法是把随机信号X(n)的N点观察数据xN(n)视为一能量有限信号，直接取xN(n)的DTFT得到XN(ejω)，然后再取其幅值的平方，并除以N，作为对真实功率谱P(ejω)的估计：表示用周期图法估计出的功率谱。因为功率谱密度直接由傅里叶变换得到，所以周期图法又称直接法。自从1965年FFT出现后，该方法就成了谱估计中的一个常用方法。将ω在单位圆上等间隔取值，得由于XN(k)可以用FFT快速计算，所以可以方便地求出。1.周期图法（直接法）周期图法是把随机信号X(n)的1.周期图法（直接法）比较以下两种计算方法：易知，直接法包含了下述假设及步骤：①把平稳随机信号X(n)视为各态遍历的，用其一个样本x(n)来代替X(n)，并且仅利用x(n)的N个观察值xN(n)来估计功率谱P(ejω)。②从记录到一个连续信号x(t)到估计出，还包括了对x(t)的离散化、必要的预处理（如除去均值和趋势项、滤波等）。1.周期图法（直接法）比较以下两种计算方法：1.周期图法（直接法）一个实际的例子（fs=250Hz）：1.周期图法（直接法）一个实际的例子（fs=250Hz2.间接法间接法的理论基础是维纳-辛钦定理。1958年Blackman和Tukey给出了这一方法的具体实现，即先由xN(n)估计出自相关函数，然后求自相关函数的傅里叶变换得到的功率谱，记之为，并以此作为对P(ejω)的估计，即因为这种方法求出的功率谱是通过自相关函数间接得到的，所以称为间接法，又称自相关法或BT法。当M

较小时，上式计算量不是很大，因此该方法是FFT问世之前常用的谱估计方法。与维纳-辛钦定理相比较：2.间接法间接法的理论基础是维纳-辛钦定理。1958年Bl2.间接法如果X(n)是各态遍历随机信号，x(n)是其一个样本函数，则自相关函数可定义如下：实际中的信号大多是因果信号，所以上式可以表示为：本章所涉及的都是自相关函数，因此将rx(m)简写为r(m)。如果观察值的个数为有限值，则求r(m)的一种方法为：由于x(n)只有N个观察值，因此对于每一个固定的延迟m，可以2.间接法如果X(n)是各态遍历随机信号，x(n)是2.间接法利用的数据只有N-1-|m|个，且在0~N-1的范围内，xN(n)=x(n)，所以实际计算时，上式变为：的长度为2N-1，它是以m=0为偶对称的。由偏差的定义可知：2.间接法利用的数据只有N-1-|m|个，且在0~N2.间接法可以看出：①对于一个固定的延迟|m|，当N→∞时，，因此

是对r(m)的渐进无偏估计；②对于一个固定的N，只有当|m|<<N时，的均值才接近于真值r(m)，即当|m|越接近于N时，估计的偏差越大；③的均值是真值r(m)和一三角窗函数的乘积，w(m)的长度是2N-1。该窗函数对r(m)加权，致使产生了偏差。2.间接法可以看出：2.间接法三角窗w(m)：当我们对一个信号做自然截短时，就不可避免地对该数据施加了一个矩形窗，由此矩形窗就产生了加在自相关函数上的三角窗，该三角窗影响自相关函数的估计质量。2.间接法三角窗w(m)：2.间接法由方差的定义可知：当N→∞时，，又因为，所以，对固定的延迟|m|，是r(m)的渐进一致估计。2.间接法由方差的定义可知：2.间接法计算时，如果N和m

都比较大，则需要的乘法次数很多。可以利用FFT实现对的快速计算。上式也可以写为：求的离散时间傅里叶变换，得：2.间接法计算时，如果N和m都比较2.间接法把xN(n)补N个零，得x2N(n)，即：记x2N(n)的傅里叶变换为X2N(ejω)，则有其中X2N(ejω)为有限长信号x2N(n)的能量谱，除以N以后即为功率谱。这说明自相关函数的估计值和x2N(n)的功率谱是一对傅里叶变换。2.间接法把xN(n)补N个零，得x2N(n)2.间接法利用FFT计算自相关函数的步骤：①对xN(n)补N个零，得x2N(n)，对x2N(n)做DFT得X2N(k)，k=0,1,…,2N-1；②求X2N(k)的幅平方，然后除以N，得；③对做逆变换，得。将中的部分向右平移2N点后形成的序列即为。2.间接法利用FFT计算自相关函数的步骤：3.直接法和间接法的关系直接法：间接法：其中自相关函数与x2N(n)的功率谱是一对傅里叶变换：因此有令M=N-13.直接法和间接法的关系直接法：令M=N-13.直接法和间接法的关系由此可知，直接法可以看作是间接法的一个特例，即当间接法中使用的自相关函数的最大延迟M=N-1时，二者是相等的。前面已经指出：这就意味着，当M较大，特别是接近于N-1时，对r(m)的估计偏差变大，此时估计出的功率谱的质量也必然下降。因此，在使用间接法时，都是取M<<N-1，很显然此时令M<<N-1，这意味着对最大长度为2N-1的自相关函数的截短，亦即施加了一个窗函数，记之为v(m)，得：3.直接法和间接法的关系由此可知，直接法可以看作是间接法的3.直接法和间接法的关系的均值等于真实的自相关函数r(m)乘以三角窗w(m)，这是第一次加窗。该三角窗是由数据截短而产生的，其宽度为2N-1。v(m)是对自相关函数r(m)的第二次加窗，宽度为2M-1，M<<N-1。因为v(m)的宽度远小于w(m)，所以v(m)的频谱V(ejω)主瓣的宽度远大于w(m)的频谱W(ejω)主瓣的宽度。这时，对r(m)施加v(m)的作用等效于在频域做和V(ejω)的卷积，这样就起到了对周期图平滑的作用。直接法和间接法往往结合起来使用，步骤如下：①对xN(n)补N个零，求；②做的傅里叶逆变换得，这时|m|≤M=N-1；③对加窗函数v(m)，这时|m|≤M<<N-1，得；④求的傅里叶变换即可得到。3.直接法和间接法的关系的均值等于真实的自3.直接法和间接法的关系直接法和间接法的关系可归纳如下：x(n)

截短乘矩形窗d0(n)求线性

相关函数DFTIDFTxN(n)

3.直接法和间接法的关系直接法和间接法的关系可归纳如下：x(3.直接法和间接法的关系一个实际的例子（fs=250Hz）：直接法间接法M=1003.直接法和间接法的关系一个实际的例子（fs=250H4.直接法和间接法估计的质量当M=N-1时，直接法和间接法估计出的结果和是相同的，因此可以把这两个估计方法的质量一起讨论。为了书写的方便，下面把和都简写为。由偏差的定义可知：

其中4.直接法和间接法估计的质量当M=N-1时，直接法和间4.直接法和间接法估计的质量三角窗w(n)是由矩形窗d0(n)做自相关得到的。记W(ejω)和D0(ejω)分别是w(n)和d0(n)的傅里叶变换，则有：因此有：当时，d0(n)趋于无限宽，W(ejω)和D0(ejω)都趋于δ

函数，这时因此，对于固定的数据长度N，是有偏估计。当时，它的期望值等于真值P(ejω)，所以它又是渐进无偏的。4.直接法和间接法估计的质量三角窗w(n)是由矩形窗d04.直接法和间接法估计的质量由协方差的定义可知：当ω=ω1=

ω2时，依据4.直接法和间接法估计的质量由协方差的定义可知：4.直接法和间接法估计的质量可以得到估计的方差为：当时，上式右边第一项趋近于零，第二项趋近于。由此可知，是真实功率谱P(ejω)的渐进无偏估计，但却不是一致估计。不管N如何选取，估计值的方差总大于等于估计值均值的平方。我们知道，是r(m)的一致估计，但把做傅里叶变换得到的功率谱却不是P(ejω)的一致估计，所以功率谱的估计要比相关函数的估计复杂的多。4.直接法和间接法估计的质量可以得到估计的方差为：如果选择一个好的窗函数，使其频谱在主瓣以外的部分基本为零（这样的窗函数是不存在的），如左下图所示，其中B1是主瓣的宽度。如果限定B1/2<ω<(π-B1/2)，则有D0(ω-λ)

D0(ω+λ)=0,如右下图所示。很显然，此时估计的方差为：4.直接法和间接法估计的质量D0(λ)

λ

π

0

λ

π

0

D0(ω+λ)

D0(ω-λ)

如果选择一个好的窗函数，使其频谱在主瓣以外的部分基本为零（这如果限定ω1和ω2在0<ω<(π-B1/2)内取值，且|ω2-ω1|>B1，如下图所示。此时估计的协方差为：可见在ω1和

ω2上是不相关的，因此呈现较大的起伏。4.直接法和间接法估计的质量λ

π

0

λ

π

0

D0(ω1-λ)

D0(λ-ω2)

D0(ω2+λ)

D0(ω1-λ)

ω1

-ω2

ω2

ω1

如果限定ω1和ω2在0<ω<(π-B1/2)内取值4.直接法和间接法估计的质量回忆平稳随机信号X(n)的自相关函数和功率谱的定义：通常求不出集总平均意义下的自相关函数和功率谱，因而假定X(n)是各态遍历的，取其一个样本x(n)，于是有:尽管自相关函数可以用时间平均来代替集总平均，但功率谱必须保留集总平均。这是因为对随机过程X(n)的每一次实现x(n)，其傅里叶变换仍然是一个随机过程，在每一个频率处，它都是一个随机变量，因此集总平均是必须的。这也说明，P(ejω)并不具备各态遍历性。4.直接法和间接法估计的质量回忆平稳随机信号X(n)的自4.直接法和间接法估计的质量对功率谱的估计：既无求均值运算，也无求极限运算，它只能看作是对真实功率谱P(ejω)做均值运算时的一个样本。缺少了统计平均，当然也就产生了较大的方差。根据也可以看出来，由单个样本x(n)估计出的功率谱不会收敛到真实功率谱。4.直接法和间接法估计的质量对功率谱的估计：4.直接法和间接法估计的质量的频率分辨率是指保持真实功率谱P(ejω)中两个靠的很近的谱峰仍被分辨出来的能力。决定分辨率的主要因素是所使用的数据的长度，也即数据窗d0(n)的宽度。若数据的长度为tp，采样频率为fs，采样后的点数为N，即tp=N/fs，那么，估计谱的分辨率正比于fs/N或2π/N。长度为N的各种窗函数，其主瓣宽度为2kπ/N，所以的分辨率也正比于2kπ/N。

若P(ejω)中有两个相距为BW的谱峰，为了要区分他们，则要求这样，数据的长度N应该满足：为了保证的分辨率，希望N要大。但增大N时，又使

起伏加剧。4.直接法和间接法估计的质量的频率分辨4.直接法和间接法估计的质量当M<N-1时，直接法和间接法估计出的结果和不相等，是对的平滑。此时：间接法也是一种有偏估计。由于在上施加了一个较短的窗口，使得间接法估计的偏差大于直接法，但方差小于直接法。谱的平滑（方差的减小）是以频率分辨率为代价的。由于V(ejω)的主瓣比W(ejω)的主瓣宽，因而致使其分辨率下降。4.直接法和间接法估计的质量当M<N-1时，直接法和间4.直接法和间接法估计的质量对，在0<ω

<(π-B1/2)的范围内，当|ω2-ω1|>B1

时，

和是不相关的，这时主瓣的宽度B1=4π/N。对，

,也可认为在上述频率范围内，当|ω2-ω1|>B1时，和是不相关的。不过此时的B1=4π/M，因为M<<N-1，所以B1增大。因此临近频率上的估计值变得较为相关。从这一角度也可以解释对平滑的原因。4.直接法和间接法估计的质量对，5.直接法估计的改进直接法的估计结果性能不好。当数据长度N太大时，谱线起伏加剧，N太小时，谱的分辨率又不好。因此需要改进。这里所说的改进，主要是改进其方差特性。间接法是对直接法的一种改进，又称为周期图的平滑。对直接法进行改进的另外一种方法是所谓的平均法，其指导思想是把长度为N的数据xN(n)分成L段，分别求每一段的功率谱，然后加以平均，以达到减小方差的目的。几种主要的改进方法：Bartlett法Welch法Nuttall法5.直接法估计的改进直接法的估计结果5.直接法估计的改进由概率论可知，对具有相同均值μ

和方差σ2的独立随机变量X1,X2,…,XL,新随机变量X=（X1+X2+…+XL）/L的均值也是μ，但方差为σ2

/L

。此即Bartlett法的指导思想。将数据分成L段，每段的长度都是M，即N=LM，第i段数据为：然后分别计算每一段的功率谱：求平均，得到平均周期图：5.直接法估计的改进由概率论可知，对具有相同均值μ和方的均值为：式中D1(ejω)是矩形窗d1(n)的频谱，W1(ejω)是d1(n)做自相关所得到的三角窗w1(n)的频谱，w1(n)的长度是2M-1。因为W1(ejω)主瓣的宽度远大于W(ejω)，所以平均后偏差增大，分辨率下降。M的选择主要取决于所需的分辨率。因为W1(ejω)的主瓣宽度是4π/M，若P(ejω)中有两个相距为BW的谱峰，为了要分辨它们，需要4π/M<BW，即

M>4π/BW。如果数据长度N

已确定，根据所需的M，L也就自然被确定。如果N可以变化，则应根据方差的要求确定L，然后再确定要记录的数据长度N。M5.直接法估计的改进的均值为：M5.直接法估计的改5.直接法估计的改进Welch法是对Bartlett法的改进。改进之一是在对xN(n)分段时，允许每一段的数据有部分的交叠。例如，若每一段数据重合一半，这时的段数：式中M仍然是每段的长度。改进之二是，每一段的数据窗口可以不是矩形窗，例如汉宁窗或汉明窗，记之为d2(n)。按Bartlett法求每一段的功率谱，记之为，即其中U是归一化因子，使用它是为了保证所得到的功率谱是渐进无偏估计。5.直接法估计的改进Welch法是对Bartlett5.直接法估计的改进如果d2(n)是一矩形窗，则平均后的功率谱为：其均值为：记D2(ejω)是d2(n)的频谱，及则有：（证明过程略）5.直接法估计的改进如果d2(n)是一矩形窗，则平均后5.直接法估计的改进Welch法中各段允许交叠，因而段数增大，这样方差就可以得到更大程度地改善。但是数据的交叠减小了每一段数据的不相关性，使方差的减小不会达到理想的程度。Welch法允许分段时交叠，这样就增加了段数，当然也就增加了做FFT的次数，如果采用的数据窗非矩形窗，这又大大增加了做乘法的次数，因此Welch法的计算量较大。Nuttall等人提出的算法步骤如下：①②与Bartlett法相同，即对xN(n)自然分段（加矩形窗），且不交叠，求得平均后的功率谱。③由做反变换，得到对应的自相关函数，其最大宽度是2M-1，M=N/L。

5.直接法估计的改进Welch法中各段允许交叠，因而段数5.直接法估计的改进④此步如同间接法，对加延迟窗w2(m)，w2(m)的最大单边宽度为M1，这样得到，即：⑤由做正变换得到对x(n)功率谱的估计：很显然，该方法是把间接法和间接法结合了起来，同时也把平滑和平均也结合了起来。这样就保持了平滑和平均减小方差的优点，同时计算量也小于Welch法。5.直接法估计的改进④此步如同间接法，对5.直接法的改进上述三种改进方法可以归纳为：x(n)

截短乘矩形窗d0(n)不交叠分段d1(n)为矩形窗求平均功率谱交叠分段对每一段加权，

d2(n)可以不是矩形窗平均做逆变换加权w2(n)正

变换BartlettWelchNuttall5.直接法的改进上述三种改进方法可以归纳为：x(n)截短5.直接法的改进仍然考虑前面的例子：Bartlett,M=100Welch,M=100Welch,M=250M1=505.直接法的改进仍然考虑前面的例子：Bartlett,M6.经典功率谱估计方法总结①不论是直接法、间接法还是直接法的改进方法，都可以利用FFT实现快速计算，且物理概念明确，仍是目前较常用的谱估计方法。②谱的分辨率较低，它正比于2π/N，N是所使用数据的长度。③由于窗函数的影响不可避免，使得功率谱P(ejω)在窗口主瓣内的功率向边瓣部分泄露，降低了分辨率。较大的边瓣有可能掩盖P(ejω)中较弱的成分，或是产生假的峰值。④方差性能不好，不是P(ejω)的一致估计，且N增大时谱线起伏加剧。⑤周期图的平滑和平均与窗函数的使用密不可分。平滑和平均主要是用来改善周期图的方差性能，但往往又减小了分辨率和增大了偏差。6.经典功率谱估计方法总结①不论是直接法、间接法还是直接法7.短时傅里叶变换对于平稳信号，前述的经典功率谱估计都是建立在传统的傅里叶变换的基础之上。其实，傅里叶变换在信号的分析中自身就存在着不足。回忆福利叶变换的表达式：显然，如果我们想知道在某一个特定时间所对应的频率是多少，或对某一个特定频率所对应的时间是多少，那么傅里叶变换则无能为力。傅里叶变换的这一缺点对统计特征随时间变化的非平稳随机信号来说，使用起来更加困难。7.短时傅里叶变换对于平稳信号，前述的经典功率谱估计都是建7.短时傅里叶变换设信号x(n)由三个不同频率的正弦首尾相接所组成，即其波形和幅频特性分别为：ω1ω2ω3nx(n)0π7.短时傅里叶变换设信号x(n)由三个不同频率的正弦首7.短时傅里叶变换短时傅里叶变换（short-timeFouriertransform,STFT）于1946年由Gabor提出，其定义为：上式中窗函数w(τ)的作用可以从不同的角度来解释：①当窗函数w(τ)沿着t轴移动时，相当于不断地截取一小段又一小段