对功率谱估计常用方法的探讨及应用

DSP课程的设计

对功率谱估计常用方法的探讨及应用分析

进行傅里叶变换在频域中研究信号，是研究确定性信号最简单且有效的手段，但在现代信号分析中，对于常见的随机信号，不可能用清楚的数学关系式来描述，其傅里叶变换更不存在，转而可以利用给定的N个样本数据估计一个平稳随机信号的功率谱密度。功率谱估计是数字信号处理的重要研究内容之一。功率谱估计可以分为经典功率谱估计和现代功率谱估计。本文介绍了各种经典功率谱估计方法，不仅从理论上对各种方法的谱估计质量进行了分析比较，而且通过Matlab进行了仿真。在对经典谱估计进行讨论之后，还分析了现代谱估计即参数谱估计方法，通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。现代谱估计的内容极其丰富，设计的学科及应用的领域都相当广泛，至今每年都有大量的科研成果出来。在本文的最后利用现代谱估计的方法讨论了功率谱方法在噪声源信号识别中的应用。文章还给出了常见谱估计方法的比较，便于深刻理解各种方法的特点，从而在实际工作中做出合理的选择。

1.功率谱方法的发展

功率谱估计是随机信号处理的重要内容，其技术渊源很长，而且在过去的40

余年中获得了飞速的发展。涉及到信号与系统、随机信号分析、概率统计、矩阵代数等一系列的基础学科，广泛应用于人民的日常生活及军事、工业、农业活动中，是一个具有强大生命力的研究领域。本文将简要回顾一下功率谱估计的发展历程，对常用的一些方法进行总结。功率谱的估计方法有很多，主要有经典谱估计和现代谱估计。经典谱估计又可以分成两种：一种是BT法，也叫间接法；另一种是直接法又称周期图法。现代谱估计的方法又大致可分为参数模型谱估计和非参数模型谱估计，前者有AR模型、MA模型、ARMA模型、PRONY模型等，后者有最小方差方法、多分量的MUSIC方法等。

1.1功率谱研究的发展过程

功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。下面对谱估计的发展过程做简要回顾:

英国科学家牛顿最早给出了“谱”的概念。后来，1822年，法国工程师傅立叶提出了著名的傅立叶谐波分析理论。该理论至今依然是进行信号分析和信号处理的理论基础。

傅立叶级数提出后，首先在人们观测自然界中的周期现象时得到应用。19世纪末，Schuster提出用傅立叶级数的幅度平方作为函数中功率的度量，并将其命名为“周期图”（periodogram）。这是经典谱估计的最早提法，这种提法至今仍然被沿用，

只不过现在是用快速傅立叶变换（FFT）来计算离散傅立叶变换（DFT），用DFT 的幅度平方作为信号中功率的度量。

1958年，R,Blackman和J.Tukey首先提出BT法，并命名为布莱克曼-杜基谱估计器（简称BT谱估计器）。这种方法是先按照有限个观测数据估计自相关函数，再对其求傅里叶变换得到功率谱。在1965年FFT未出现以前，BT法一直是最常用的方法。

周期图较差的方差性能促使人们研究另外的分析方法。1927年，Yule提出用线性回归方程来模拟一个时间序列。Yule的工作实际上成了现代谱估计中最重要的方法——参数模型法谱估计的基础。

Walker利用Yule的分析方法研究了衰减正弦时间序列，得出Yule-Walker方程，可以说，Yule和Walker都是开拓自回归模型的先锋。

1930年，著名控制理论专家Wiener在他的著作中首次精确定义了一个随机过程的自相关函数及功率谱密度，并把谱分析建立在随机过程统计特征的基础上，即，“功率谱密度是随机过程二阶统计量自相关函数的傅立叶变换”，这就是Wiener—Khintchine定理。该定理把功率谱密度定义为频率的连续函数，而不再像以前定义为离散的谐波频率的函数。

1949年，Tukey根据Wiener—Khintchine定理提出了对有限长数据进行谱估计的自相关法，即利用有限长数据估计自相关函数，再对该自相关函数球傅立叶变换，从而得到谱的估计。1958年，Blackman和Tukey在出版的有关经典谱估计的专著中讨论了自相关谱估计法，所以自相关法又叫BT法。周期图法和自相关法都可用快速傅立叶变换算法来实现，且物理概念明确，因而仍是目前较常用的谱估计方法。

1948年，Bartlett首次提出了用自回归模型系数计算功率谱。自回归模型和线性预测都用到了1911年提出的Toeplitz矩阵结构，Levinson曾根据该矩阵的特点于1947年提出了解Yule-Walker的快速计算方法。这些工作为现代谱估计的发展打下了良好的理论基础。1965年，Cooley和Tukey提出的FFT算法，也促进了谱估计的迅速发展。

现代谱估计的提出主要是针对经典谱估计（周期图和自相关法）的分辨率和方差性能不好的问题。1967 年，Burg 提出的最大嫡谱估计，即是朝着高分辨率谱估计所作的最有意义的努力。虽然，Bartlett 在1948年，Parzem 于1957 年都曾经建议用自回归模型做谱估计，但在Burg 的论文发表之前，都没有引起注意。

现代谱估计的内容极其丰富，涉及的学科及应用领域也相当广泛，至今，每年都有大量的论文出现。目前尚难对现代谱估计的方法作出准确的分类。从现代谱估

计的方法上，大致可以分为参数模型谱估计和非参数模型谱估计，前者有 AR 模型、MA 模型、ARMA 模型、PRONY 模型等；后者有最小方差方法、多分量的 MUSIC 方法等。非参数模型谱估计的特点是其模型不是用有限参数来描述，而直接由相关函数序列得到，这种方法能提高低信噪比时的谱分辨率。参数模型谱估计是先根据过程的先验信息或者一些假定，建立一个数学模型来表示所给定采样数据的过程，或者选择一个较好的近似实际模型，而后利用采样数据序列或者自相关序列，估计该模型的参数，最后把参数代入到该模型对应的理论功率谱表达式，得到所需要的谱估计。目前大量的论文集中在模型参数的求解上，以求得到速度更快、更稳健、统计性能更好的算法。

1.2 功率谱估计方法提出

在通信系统中，往往需要研究具有目中统计特性的随机信号。由于随机信号是一类持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能够应用确定信号的频谱计算方法去分析随机信号的频谱。然而，虽然随机信号的频谱不存在，但其相关函数是可以确定的。如果随机信号是平稳的，那么其相关函数的傅里叶变换就是它的功率谱密度函数，简称功率谱。功率谱反映了单位频带内随机信号的一个样本信号来对该随机过程的功率谱密度函数做出估计。

1.3 功率谱估计应用及用途

功率谱估计有着极其广泛的应用，不仅在认识一个随机信号时，需要估计它的功率谱。它还被广泛地应用于各种信号处理中。在信号处理的许多场所，要求预先知道信号的功率谱密度(或自相关函数)。例如，在最佳线性过滤问题中，要设计一个维纳滤波器就首先要求知道(或估计出)信号与噪声的功率谱密度(或自相关函数)。根据信号与噪声的功率谱(或)(m xx φ)才能设计出能够尽量不失真的重现信号，而把噪声最大限度抑制的维纳滤波器。

常常利用功率谱估计来得到线性系统的参数估计。例如，当我们要了解某一系统的幅频特性)(ωH 时，可用一白色噪声()n ω通过该系统。再从该系统的输出样本y(n)估计功率谱密度)(ωyy P 。由于白色噪声的PSD(用)(ωωωP 表示)为一常数即

2

)(ωωωσω=P ，于是有：

2

2

)()(ωσωωH P yy = (1-1)

故通过估计输出信号的PSD ，可以估计出系统的频率特性)(ωH (模特性)。从宽带噪声中检测窄带信号。这是功率谱估计在信号处理中的一个重要用途。但是这要求功率谱估计有足够好的频率的分辨率，否则就不一定能够清楚地检测出来。所

谓谱估计的分辨率可以粗略地定义为能够分辨出的二个分立的谱分量间的最小频率间隙(距)。提高谱估计的分辨率已成为目前谱估计研究中的一个重要方向

功率谱估计就是通过信号的相关性估计出接受到信号的功率随频率的变化关系，实际用途有滤波，信号识别（分析出信号的频率），信号分离，系统辨识等。谱估计技术是现代信号处理的一个重要部分，还包括空间谱估计，高阶谱估计等。维纳滤波、卡尔曼滤波，可用于自适应滤波，信号波形预测等（火控系统中的飞机航迹预判）。

2 经典谱估计

2.1 周期图法

周期图法又称直接法。它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样.

周期图这一概念早在1899年就提出了，但由于点数Ｎ一般比较大，该方法的计算量过大而在当时无法使用。只是1965年FFT 出现后，此法才变成谱估计的一个常用方法。周期图法包含了下列二条假设：

1．认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段

)(n x N 来估计该随机序列的功率谱。这当然必然带来误差。

2．由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。与相关法相比，相关法在求相关函数)(m R x 时将

)(n x N 以外是数据全都看成零，因此相关法认为除)(n x N 外x(n)是全零序列，这种处

理方法显然与周期图法不一样。

但是，当相关法被引入基于FFT 的快速相关后，相关法和周期图法开始融合。比我们发现：如果相关法中M=N ，不加延迟窗，那么就和充（N-1）个零的周期图法一样了。简单地可以这样说：周期图法是M=N 时相关法的特例。因此相关法和周期图法可结合使用。

2.2 相关法谱估计 （BT ）

这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。它是1958年由Blackman 和Tukey 提出。这种方法的具体步骤是：

第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列)(n x N 第二步：由N 长序列)(n x N 求（2M-1）点的自相关函数)(m R x

序列。即

)()(1

)(1

m n x n x

N

m R N n N N

x +=

∑-=∧

(2-1)

这里，m=-(M-1)…,-1,0,1…,M-1，M N ，)(m R x 是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。。。,M-1的傅里叶变换，另一半也就知道了。

第三步：由相关函数的傅式变换求功率谱。即

jwm

M M m X

jw

x e m R

e S ----=∧∧

∑=

)()(1)

1( (2-2) 以上过程中经历了两次截断，一次是将x(n)截成N 长，称为加数据窗，一次是将想x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。一般取M<

当FFT 问世后，情况有所变化。因为截断后的)(n x N 可视作能量信号，由相关卷积定理可得

)(m R x ∧

)]()([1

m x m x N

N N -*=

(2-3) 这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。我们可对上式两边取（2N-1）点DFT ，则有

2121212)(1)()([1)(K X N

K x k x N m R N N N x ---∧

=*=

(2-4) 于是将时域卷积变为频域乘积，用快速相关求)(m R x ∧

的完整方案如下：

1. 对N 长)(n x N 的充（N-1）个零，成为（2N-1）长的。

2.

求（2N-1）点的FFT ，得∑-=----=2

20

121

212)()(N N mk N N N W n x

K X 。

3.

求

2

12)(1K X N N -。由DFT 性质，)(12n x N -是纯实的，)(12k x N -满足共轭偶对称，而2

12)(1K X N

N -一定是实偶的，且以（2N-1）为周期。

4. 求（2N-1）点的IFFT ：

mk N N N k N x W K X N N m R -----=-∧

∑-=121

)1(2

12)(1121)( (2-5)

这里

2

12)(1K X N N -是实偶的，

m=-(N-1)...0...N-1。本来IFFT 求和范围是0至2N-2，由于2

12)(1K X N

N -的实偶性与周期性，求和范围改为-（N-1）至（N-1）不影响计

算结果。同理可将m 的范围改为-（N-1）至（N-1）。

上述的快速相关中，充零的目的是为了能用圆周卷积代替线性卷积，以便进一步采用快速卷积算法。快速相关输出是-（N-1）至（N-1）的2N-1点，加)(m W M 窗后截取的是-（M-1）至（M-1）的频段，最后作（2M-1）点FFT ，得)(k S x ∧

。我们注意到：如果数据点数与自相关序列点数相同即M=N ，则（2N-1）点的IFFT 后紧跟一个（2N-1）点的FFT ，利用)(m R x ∧

的对称性，FT 运算框的计算式变为

=

)(K S X ∑---=∧

1)

1()(N N m X

m R

1212---N mk N X W (2-6)

由于N=M 并假设窗形状是矩形的，第二次()m W M 的截断就不需要了。比较式

（3-5）和式（3-6）， ](m)R FFT[k S x x ∧

=）

（,])(1[)( 2

12K X N

IFFT m R N x -∧

=正反傅氏变换可以抵消，直接得

）

（k S x =2

12)(1K X N

N - (2-7) 为了实行基2FFT ，也可将（2N-1）点换成2N 点，这样做不影响结果的正确性。

2.3 巴特利特(Bartlett)平均周期图的方法

首先让我们来看一下为什么周期图经过某种平均(或平滑)后会使它的方差当∞→N 时趋于零，达到一致估计的目的。如果L x x x , , ,21 是不相关的随机变量，每一个具有期望值μ，方差2σ，则可以证明它们的数学平均L x x x x l /)...(21+++=的期望值等于μ，数学平均的方差等于L /2σ，即：

[][]μμ=?=+++=

L L

x x x E L x E L 1

121 [][]

[]222)(][))((x E x E x E x E x Var -=-=

[]

22212)(1

μ-+++=

L x x x E L

()[]

2112

22212

1μ-??

???

??

???++++=∑∑=≠=L j L j i i j

i L x x E x x x E L

[][][]∑∑∑∑=≠==≠=?=L j L

j

i i i

j

L j j

i i j

i x E x E x x E 1

111

222)1(μμμμL L L L -=-=

所以

[][][]

??

????-=-??????-+=∑∑==21222

22122211μμμμL x E L L L x E L x Var L i i L i i

[][][][][][]{}

2

2222221212

])[(])[(])[(1L L x E x E x E x E x E x E L

-++-+-=

L L L

2221σσ== (2-8)

由(2-8)可见，L 个平均的方差比每个随机变量的单独方差小L 倍。当

[]0 →∞→x Var L ，可达到一致谱估计的目的。因而降低估计量的方差的一种有效方法是将若干个独立估计值进行平均。把这种方法应用于谱估计通常归功于Bartlett 。

Bartlett 平均周期图的方法是将序列)10( )(-≤≤N n n x 分段求周期图再平均。设将)(n x 分成L 段，每段有M 个样本，因而LM N =，第i 段样本序列可写成

L i M n M iM n x n x i ≤≤-≤≤-+=1 ,10 )()(

第i 段的周期图为

2

10

)(1)(∑-=-=

M n n j j

i

M e n x

M

I ωω

如果)( ,m M m xx φ>很小，则可假定各段的周期图)(ωi

M I 是互相独立的。对功率

谱密度的概念的讨论，谱估计可定义为L 段周期图的平均，即

∑==L i i

M xx I L P 1

)(1)(?ωω (2-9)

于是它的期望值为

[]

[][]

∑===L i i M i M xx I E I E L P E 1

)()(1)(?ωωω

[]

?-

-=π

π

θωθθπ

ωd e W P P E j B xx xx

)()(21)(?)(

()[()[?-??

??

??--=

π

πθθωθωθπd M P M

xx 2

2/sin 2sin )(21 (2-10) 这里L N M /=，因此Bartlett 估计的期望值是真实谱)(ωxx P 与三角窗函数的卷积。由于三角窗函数不等于δ函数，所以Bartlett 估计也是有偏估计即0≠Bias ，但当∞→N 时，0→Bias 。

由于我们假定各段周期图是相互独立的，所以可按式(3-8)得到下式：

[]

[])(1)(?ωωM

xx I Var L

P Var = ???

???????? ??+≈22

sin )sin(1)(1ωωωM M P L xx (2-11)

由此可见，随着L 的增加[]

)(?ωxx P Var 是下降的，当∞→L 时，[]

0)(?→ωxx

P Var 。因此Bartlett 估计是一致估计。

比较式(2-10)的[]

)(?ωxx

P E 与式(2-1)的[])(ωN I E 可见在二种情况的估计量的期望值都是真值)(ωxx P 与窗口函数)(ωj B e W 的卷积形式，但后者将前者WB 中N 改为M ，N L N M <=/。

因而使)(ωj B e W 主瓣的宽度增大。由于主瓣的宽度愈窄愈接近δ函数，偏倚愈小。今式(2-10)中)(ωj B e W 的主瓣宽度大于后式中的主瓣宽度，因而

[]

[])()(?ωωN

xx I Bias P Bias >，而主瓣愈宽显然使分辨率愈差。因此Bias 可用来说明谱的分辨率，Bias 愈大说明谱分辨率愈差。一个固定的记录长度N ，周期图分段的数目L 愈大将使方差愈小，但M 也愈小，因而使Bias 愈大，谱分辨率变得愈差。因此Bartlett 方法中Bias 或谱分辨率和估计量的方差间是有互换关系的。M 和N 的选择一般是由对所研究的信号的预先了解来指导的。例如，如果我们知道谱有一个窄峰，同时如果分辨出这个峰是重要的，那么我们必须选择M 足够大。又从方差的表达式我们可以确定谱估计的可接受的方差所要求的记录长度N=(LM)。

由此可见Bartlett 法使谱估计的方差减小是用增加Bias 以及降低谱分辨率的代价换来的。实际上，当N 一定时，Bias 与Var 的互换性是谱估计的一个固有特性。

为了说明经平均后的周期图作为功率谱估计的实际效果，设有一零均值高斯分布的随机过程，其功率谱密度为

1.01 )

1)(1(1

)(1

-=--=

-a az az z P xx

这一功率谱密度是由一零值高斯分布单位方差的噪声序列通过一个其

)1/(1)(1--=az z H 的滤波器形成的。为了简单设选用8=N 的矩形窗函数。其周期图的

期望值(用)]([8w I E 表示)与真值)(w P xx 均表示在图2.4中。它说明8=N 的周期图可以得到的Bias 的情况。图2-2表示M=8分4段与16段二种情况平均后的周期图。显然L=4的方差比L=16的大。将L=16的曲线与图2.4的曲线比较可见在这种估计中误差的大部分起因于Bias 而不是方差(因为二种情况均为8=M ，Bias 相仿，误差也相仿)。M=16，L=2及8的周期图表示在图2-3。图2-3中L 不同造成的影响也是明显的。但是这二个曲线的起伏都很大，因此有理由认为误差主要起因于方差。比较

8=M 与

M=16的周期图可见， M 愈大(即§2.3中的N 愈大)，将使周期图的起伏愈

增快的结论，在这里也同样成立。

比较图2-2与图2-3发现在这个例子中最好的选择是应用L=16，M=8的估计而不是L=8、M=16的估计，即宁可减小方差，牺牲Bias 。在实际中，当然功率谱密度的真值是不知道的。但是谱的窗口函数以及关于功率谱密度的某些信息往往是预

先知道的。通过改变M 和L 以及利用预先知道的情况，通常可以很好地进行选择。

平均周期图的方法特别适合于应用FFT 算法。因此在FFT 出现以后这个方法比下面将要讨论的利用窗口函数的处理法用得更多。而在FFT 出现以前主要是用窗口函数处理法平滑周期图。

图2-1 )(ωx P 与)]([8ωI E 的特性 图2-2 平滑后的周期图（每段取8个数据）

图2-3 平均后的周期图(每段取16个数据)

2.4 Welch 法

Welch 提出了对Bartlett 法的修正使更适合于用FFT 进行计算。他主要提出二方面的修正，其一是选择适当的窗函数)(n ω，并在周期图计算前直接加进法，这样得到的每一段的周期图为

)()(1)(21

)

(∑-=-=

N n n

j i i M e

n n x MU

I ωωω (2-12)

这里∑-==

1

2

)(1

N n n W

M

U 为归一化因子，而Bartlett 法每段的周期图为

21

)

()(1)(∑-=-=

N n n

j i i M e

n x M

I ωω (2-13)

这样加窗函数的优点是无论什么样的窗函数均可使谱估计非负。其二是在分段时，可使各段之间有重迭，这样将会使方差减小(当N 与M 一定时)。他建议可以重迭50％。

3 经典谱估计方法编程与分析

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));

图3-1 周期图法中N=1000功率谱

图3-2 周期图中N=100000功率谱

图3-3 周期图当N=100功率谱

随着采样点数的增加，该股集是渐进无骗的。从图中可以看出，采用周期突发估计得出的功率谱很不平滑，相应的估计协方差比较大。而且采用增加采样点的办法也不能吃周期图变得更加平滑，这是周期图法的缺点。周期图法得出的估计谱方差特性不好：当数据长度N 太大时，扑线的起伏加剧；N 太小时谱的分辨率又不好。对其改进的主要方法有二种，即平均和平滑，平均就是将截取的数据段)(n x N 再分成L 个小段，分别计算功率谱后取功率谱的平均，这种方法使估计的方差减少，但偏差加大，分辨率下降。平滑是用一个适当的窗函数)(jw e W 与算出的功率谱

)(jw e S X

进行卷积，使谱线平滑。这种方法得出的谱估计是无偏的，方差也小，但

分辨率下降。

3.2 间接法

间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 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);

图3-4 直接发功率谱图

当M=N-1时,BT 法与周期图法估计出的功率谱是一样的;当 M < N-1 时,BT 法的偏差大于周期图法,在窗函数满足一定条件时是渐进无偏估计;方差小于周期图的方差;分辨率比周期图法低,与窗函数的选择有关。

BT 法的缺点在于当 M →N 时, )(m R ∧

的方差很大,使谱估计质量下降;由)(m R ∧

得到的)(w S ∧

不一定为正值,从而可能失去功率谱的物理意义。

3.3 改进的直接法：

对于直接法的功率谱估计，当数据长度N 太大时，谱曲线起伏加剧，若N 太小，谱的分辨率又不好，因此需要改进。

3.3.1 Bartlett 法

对周期图法改进的思想是将信号分段进行估计，然后再将这些估计结果进行平均，从而减小估计的协方差，是功率谱图变得比直接法更平滑。增加分段数可以进一步减低估计的协方差，然而若每段中的数据点数太少，就会使估计的频率分辨率

下降很多。从样本信号序列总点数一定的条件下，可以采用使分段相互重叠的方法来增加分段数，从而保持每段信号点数不变，这样就在保证频率分辨率的前提下进一步降低估计协方差。而一般的Bartlett法是通过降低分辨率来降低其方差的。

Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。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_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

图3-5 bartlett法估计功率谱图

3.3.2Welch法

现在比较常用的改进方法是Welch 法，又叫加权交叠平均法，简记为WOSA 法，这种方法以加窗（加权）求取平滑，以分段重叠求得平均，因此集平均与平滑的优点于一体，同时也不可避免带有两者的缺点，因此归根到底是一种折中。其主要步骤是：

（1）将N 长的数据段分成L 个小段，每小段M 点，相邻小段间交叠M/2点（即2:1分段）。因为L （M/2）+M/2=N,所以段数

2

/2

/M M N L -=

(3-1)

（2）对各小段加同样的平滑窗

后起傅氏变换

,,...1,)()()(10

L i e n ax n x e X jwn M n i jw

i ==--=∑ （3-2）

（3）用下式求各小段功率谱的平均

21)(1

1)(jw i l i jw

i e X MU L e S ∑=∧

= （3-3）

这里，代表窗函数平均功率，MU 是M 长窗函数的能量。

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';

[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);

pause;

figure(2)

plot(f,plot_Pxx1);

pause;

figure(3)

plot(f,plot_Pxx2);

图3-6 矩形窗

图3-7 海明窗

图3-8 blackman窗

Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。

不同窗函数的Welch 谱估计在选择窗函数时，一般有如下要求：

(1)窗口宽度M要远小于样本序列长度N ，以排除不可靠的自相关值；

(2)当平稳信号为实过程时，为保证平滑周期图和真是功率谱也是实偶函数，平滑窗函数必须是实偶对称的；

(3)平滑窗函数应当在m=0出游峰值，并且m随绝对值增加而单调下降，使可靠的自相关值有较大的权值；

(4)功率谱是频率的非负函数且周期图是非负的，因而要求窗函数的Fourier变换是非负的。

仍以上述平稳随机信号为例,采样频率、采样点数、FFT点数、窗长度及重叠数据不变,窗函数采用矩形窗、Blackman 窗、Hamming窗，仿真结果如图所示。

由仿真结果可以看出：使用不同的窗函数谱估计的质量是不一样的，矩形窗的主瓣较窄，分辨率较好，但方差较大，噪声水平较高；而Blackman 窗和Hamming 窗的主瓣较宽，分辨率较低，但方差较小，噪声水平较低。因此，在进行谱分析时选择何种窗函数，要视具体情况而定。如果强调高分辨率，能精确读出主瓣频率，而不关心幅度的精度，例如测量震动物体的自震频率，可以选用主瓣宽度比较窄的矩形窗；对受到强干扰的窄带信号若干扰靠近信号，则可选用旁瓣幅度较小的窗函数，若离开通带较远，则可选用渐近线衰减速度比较快的窗函数。总之，要针对不同的信号和不同的处理目的来选择合适的窗函数，这样才能得到良好的效果。

3.4经典谱估计方法的比较总结

（1）周期图法和BT法的特点是离散性大，曲线粗糙，方差较大，但是分辨率较高；

（2）Bartlett法和Welch法的收敛性较好，曲线平滑，方差较小，但是功率谱主瓣较宽，分辨率低，这是由于对随机序列加窗截断所引起的Gibbs效应造成的；

（3）与Bartlett 法相比，Welch 法的估计曲线比较粗糙，但是分辨率较好，原因是 Welch 法中对数据进行截断时加的是Hanning 窗，而在Bartlett 法中使用的是矩形窗，相对于矩形窗，Hanning 窗的主瓣包含更多的能量，因而使功率谱的主瓣较窄，分辨率较高。

由上述理论分析及仿真实验可知，Welch 法采用加窗交叠求功率谱，可以有效减小方差和偏差，一般情况下能接近一致估计的要求，因而得到广泛应用。同时还可以发现，对信号加不同的窗函数，谱估计的质量是不同的。

在经典谱估计中,无论是周期 图法还是其改进的方法,都存在着频率分辨率低、方差性能不好的问题，原因是谱估计时需要对数据加窗截断，用有限个数据或其自相关函数来估计无限个数据的功率谱，这其实是假定了窗以外的数据或自相关函数全为零，这种假定是不符合实际的，正是由于这些不符合实际的假设造成了经典谱估计分辨率较差。另外，经典谱估计的功率谱定义中既无求均值运算又无求极限运算，因而使得谱估计的方差性能较差，当数据很短时，这个问题更为突出。如何选取最佳窗函数、提高频谱分辨率， 如何在短数据情况下提高信号谱估计质量 ，还需要进一步研究。

4.现代谱估计

现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱,主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。常用模型有 ARMA 模型、 AR 模型、 MA 模型,AR 模型的正则方程是一组线性方程，而MA 和ARMA 模型是非线性方程。由于AR 模型具有一系列良好的性能，因此被研究最多也得到最广泛的应用。本节将较为详细的讨论AR 模型，并对MA 和ARMA 模型谱估计方法做简要的讨论。

4. 1 AR 模型

AR 模型又称为自回归模型,是一个全极点的模型,可用如下差分方程来表示:

)

()()(1n w k n x n x p

k k a +--=∑=

（4-1）

其中，p 是 A R 模型的阶数， {a k }=l ， 2 ， …， p 是p 阶 A R 模型的参数．将该模型记为AR( p) ，它的系统转移函数为

∑=-+==

p

k k

k z

a z W z X z H 111)

()

()( (4-2)

在功率谱估计中，若观测的数据 x(n) 是平稳随机过程，则该系统的输入w( n ) 也可认为是平稳的，因 而根据线性系统对平稳随机信号的响应理论可得观测数据的功率谱为

2

1

22

2

|

1||)(|)(e

a e P jkw p

k k w

jw

w x

H w -=∑+==σ

σ (4-3)

由式可知，利用A R 模型进行功率谱估计的实质是求解模型系数 {a k } 和σ2

w 的问 题．

将式 ( 1 )两端乘以x ( n-m ) 求平均 ( 数学期望) ，可以求得观测数据的A R( p) 模型参数与自相关函数的关系式为

000)()()

()(21

1

<=>????

?

????

-+----=∑∑==m m m m k m k m m R R a R a R xx

w

xx p k k xx p

k k xx σ (4-4) 可见， 阶 AR 模型输出的相关函数具有递推的性质， 因而选用 AR 模型进行谱估计只需较少的观测数据将式 ( 4 )写成矩阵形式得

??????????????=

??????????????????????????--001)0()1()()1()0()1()()1()

0(21 σw p a a R p R p R p R R R p R R R (4-5) 上式就是著名的Yule-Walker( Y —W)方程．它表明，只要已知观测数据的自相关数，就能求出AR 模型参数{a k } 和σ2

w ，进而按式 ( 3 )求得信号功率谱的估值。

另外，从 AR 模型的差分方程式可知，该模型的现在输出值是它本身过去值的回归，这与预测器存在着一定的相似性，它们之间有着非常密切的关系，即它们的系统函数互为倒数，也就是说预测误差滤波器

)(z A p

就是 AR 信号模型H(z)

的逆滤波器．因此通过预测误差滤波器优化设计使预测均方误差最小就可求得A R 信号模型的最优参数，即P 阶线性预测器的预测系数{)(k a p }等于p 阶A R 模型的系数{a k }，其最小均方预测误差E p 等于自噪声方差σ2

w 。

因此， 根据上述的Y-W 方程以及 AR 模型与预测误差滤波器之间的关系，就可提