现代功率谱估计淮北师范大学物理与电子信息学院 235000摘要功率谱估计就是基于有限的数据寻找信号、随机过程或系统的频率成分。
它是随机信号处理的重要内容，广泛应用于人民的日常生活及军事、工业、农业活动中。
其实现方法主要可分为经典谱估计和现代谱估计。
经典谱估计方法由于其种种缺点，迫使人们大力研究现代谱估计方法。
现代谱估计法是以参数模型为基础的方法，大致可以分为参数模型谱估计和非参数模型谱估计，前者有AR模型、MA模型、ARMA模型、PRONY模型等；后者有最小方差方法、多分量的MUSIC 方法等。
本文将着眼于现代谱估计的各种方法，首先简要介绍随机信号功率谱估计的相关基础知识，然后从经典法入手，探讨现代谱估计的理论基础，分析各种方法的优劣性及适用范围，并且给出对应的Matlab仿真结果，从而深刻理解各种方法的特点，从而在实际工作中做出合理的选择。
关键词功率谱估计现代信号处理 Matlab引言功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。
英国科学家牛顿最早给出了“谱”的概念。
后来，1822年，法国工程师傅立叶提出了著名的傅立叶谐波分析理论。
该理论至今依然是进行信号分析和信号处理的理论基础。
傅立叶级数提出后，19世纪末，Schuster提出用傅立叶级数的幅度平方作为函数中功率的度量，并将其命名为“周期图”（periodogram）。
这是经典谱估计的最早提法，这种提法至今仍然被沿用。
周期图较差的方差性能促使人们研究另外的分析方法。
1927年，Yule提出用线性回归方程来模拟一个时间序列。
Yule的工作实际上成了现代谱估计中最重要的方法——参数模型法谱估计的基础。
Walker利用Yule的分析方法研究了衰减正弦时间序列，得出Yule-Walker方程，可以说，Yule和Walker都是开拓自回归模型的先锋。
1948年，Bartlett首次提出了用自回归模型系数计算功率谱。
自回归模型和线性预测都用到了1911年提出的Toeplitz矩阵结构，Levinson曾根据该矩阵的特点于1947年提出了解Yule-Walker的快速计算方法。
这些工作为现代谱估计的发展打下了良好的理论基础。
1965年，Cooley和Tukey提出的FFT算法，也促进了谱估计的迅速发展。
现代谱估计的提出主要是针对经典谱估计（周期图和自相关法）的分辨率和方差性能不好的问题。
1967 年，Burg 提出的最大熵谱估计，即是朝着高分辨率谱估计所作的最有意义的努力。
由于随机信号是一类持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能够应用确定信号的频谱计算方法去分析随机信号的频谱。
然而，虽然随机信号的频谱不存在，但其相关函数是可以确定的。
如果随机信号是平稳的，那么其相关函数的傅里叶变换就是它的功率谱密度函数，简称功率谱。
功率谱反映了单位频带内随机信号的一个样本信号来对该随机过程的功率谱密度函数做出估计。
本文将着眼于现代谱估计的各种方法，首先简要介绍随机信号功率谱估计的相关基础知识，然后从经典法入手，探讨现代谱估计的理论基础，分析各种方法的优劣性及适用范围，并且给出对应的Matlab仿真结果，从而深刻理解各种方法的特点，从而在实际工作中做出合理的选择。
1 从经典谱估计到现代谱估计首先，给出功率谱的两个最基本的定义如下：S xx(e jω)=∑R xx(k)e−jωk∞k=−∞(1)P x(e jω)=limM→∞E{12M+1|∑x(n)e−jωnMn=−M|2} (2)可以证明，这两个定义是等效的。
⋯无论是建立在第一个还是第二个公式上的定义，在实际中都几乎是不可能实现的（除非x(n)可以用解析法精确的表示），因此，只能用所得的有限次记录（往往仅一次）的有限长数据来予以估计，这就产生了功率谱估计这一极其活跃，同时也极其重要的研究领域。
1 经典谱估计之周期图法在实际应用中，通常观测到的是信号的有限个(例如N个)取样值，用y N(n)表示。
可以认为它是分段平稳随机信号的一段，也可将它看成是从平稳随机信号中截取的一段数据。
对于平稳随机信号，无论从何时开始任取一段长为N的数据，所计算出来的均值或自相关函数都是相同的。
信号y N(n)可以看成是用一个宽为N的数据窗w(n)从平稳随机信号y N(n)中截取出来的，即y N(n)=y(n)w(n) (3)根据遍历性，用时间平均代替集合平均，若已知N个数据为y N(n)=(y0,y1,…,y N−1) (4)则用时间平均来近似计算的自相关函数为R̂yy(k)=1N∑y n+k y kN−1−|k|n=0,|k|≤N−1 (5)称之为取样自相关。
它可以看成是有限长序列y N(n)和y N(−n)的卷积运算结果除以N，即R̂yy(k)=1N(y N(n)∗y N(−n)) (6)取样自相关函数的双边Z 变换叫做周期图，它是功率谱的一种估计，用 Ŝyy (z ) 表示， S ̂yy (z )=∑R̂yy (k )N−1k=−(N−1)z −k (7) 联系式(9)，由上式得到S ̂yy (z )=1NY (z )Y (z −1) (11) 这里Y (z )是y N (n )的Z 变换。
式(10)和(11)是计算周期图的两种基本方法，前者称为间接法，后者成为直接法。
令z =e jω,由式(11)得到Ŝyy (z )=1N |Y(ω)|2=1N |∑y N (n)e −jωn N−1n=0|2 (12) 该式很适合用FFT 计算。
改进周期图的有4种办法：修正周期图法，平均周期图法，加床平滑法，Welch 法。
2 经典谱估计之自相关法根据Wiener-Khintchine 定理，平稳离散随机信号x(n)的自相关函数R xx (m )=E[x ∗(n )x (m +n )] (3)与功率谱S xx (ω)之间构成一对傅里叶变换关系，即S xx (ω)=∑R xx (m)e −jωm ∞m=−∞(4)R xx (m )=12π∫S xx (ω)e jωm dωπ−π(5) 这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。
它是1958年由Blackman 和Tukey 提出。
这种方法的具体步骤是：第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列)(n x N第二步：由N 长序列)(n x N 求（2M-1）点的自相关函数)(m R x序列。
即)()(1)(10m n x n xN m R N n N N x +=∑-=∧ (6)这里，m=-(M-1)…,-1,0,1…,M-1，M N ，)(m R x 是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,…,M-1的傅里叶变换，另一半也就知道了。
第三步：由相关函数的傅式变换求功率谱。
即jwm M M m Xjw x e m R e S ----=∧∧∑=)()(1)1( (7)以上过程中经历了两次截断，一次是将x(n)截成N 长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。
因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。
一般取M<<N ，因为只有当M 较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。
因此，在FFT 问世之前，相关法是最常用的谱估计方法。
当相关法被引入基于FFT 的快速相关后，相关法和周期图法开始融合。
简单地可以这样说：周期图法是M=N 时相关法的特例。
因此相关法和周期图法可结合使用。
周期图和自相关法及它们的改进方法称为谱估计的经典方法。
然而，传统方法并不是功率谱的良好估计。
事实上，随着记录长度的增加，这两种估计的随机起伏反而更加严重。
此外，他们还有着频率分辨率不高和旁瓣泄露两个难以克服的缺点。
这就促进了现代谱估计方法研究的展开。
2 现代谱估计的原理及方法现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱,主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。
常用模型有 ARMA 模型、 AR 模型、 MA 模型。
由于AR 模型具有一系列良好的性能，因此被研究最多也得到最广泛的应用。
本节将较为详细的讨论AR 模型，并对MA 和ARMA 模型谱估计方法做简要的讨论。
以参数模型为基础的谱估计方法一般按下列3步进行：(1)选择合适的信号模型。
(2)根据已知的有限个观测数据，或者它的有限个自相关函数估计值，估计模型的参数；(3)用估得的模型参数计算输出功率谱。
下面解释参数模型法的理论基础。
1 ARMA模型如图1所示，输入激励u(n)是均值为零，方差为σ2的白噪声序列。
线性系统传输函数为H(z)=B(z)A(z)=∑b k z−kqk=0∑a k z−kpk=0式中，b k是前馈支路的系数，称为MA系数；a k是反馈支路的系数，称为AR系数。
系统的输出序列是被建模的离散随机信号。
u(n(n)图1 离散随机信号x(n)的有理传输函数模型该模型的输出功率谱和输入功率谱之间存在下列关系：S xx(z)=σ2 H(z)H∗(1z∗)=σ2B(z)B∗(1z∗)A∗(1z∗)或者S xx(e jω)=σ2 |H(e jω)|2=σ2|B(e jω) A(e jω)|2设a0=1和b0=1，其余所有的系数不全为零，这种模型称为ARMA(p,q)模型。
2 AR模型若除b0=1外所有其他的MA系数都等于零，则称为p阶自回归模型或简称为AR(p)模型。
AR模型的传输函数为H AR(z)=1A(z)=11+∑a k z−kpk=1模型输出功率谱为S xx (z )=σ2A (z )A (z −1)或S xx (e jω)=σ2|A (e jω)|2=σ2|1+∑a ke −jωk p k=1|2这是一个全极点模型。
AR 模型的Yule-Walker 方程如下：R xx (m )={−∑a k R xx (m −k )+σ2p k=1 , m =0−∑a k R xx(m −k )p k=1 , m >0 自相关函数的头p+1个值是{R (0),R (1),…,R (p )}，因此，式__表示成下列矩阵形式：[ R (0)R (1)R (2) … R (p )R (1)R (0)R (1) … R (p −1)R (2) R (1) R (0) … R (p −2)⋮ ⋮ ⋮ ⋱ ⋮R (p )R (p −1)R (p −2)… R (0)] [ 1a 1a 2⋮a p ] = [ σ200⋮0]这就是AR(p)模型的Yule-Walker 方程。
![功率谱计算[解说]](https://img.taocdn.com/s1/m/3cec2a18f02d2af90242a8956bec0975f465a478.png)