现代数字信号处理课件：功率谱估计

功率谱估计4.1经典功率谱估计4.2AR模型功率谱估计的方法和性质4.3最大熵谱估计方法4.4最大似然谱估计4.5互协方差估计与互谱估计4.6特征分解法谱估计4.1经典功率谱估计

4.1.1BT法

根据维纳-辛钦定理，1958年Blackman和Tukey给出了这一方法的具体实现，即先由N个观察值xN(n)估计出自相关函数r

x(m)，求其傅立叶变换，以此变换结果作为对功率谱

Px(ω)的估计。

如果我们得到的是x(n)的N个观察值x(0)，x(1)，…，x(N－1)，令

xN(n)=a(n)·x(n)

(4.1.1)

其中a(n)是数据窗，对于矩形窗

计算rx(m)的估计值的一种方法是(4.1.2)

的均值为若a(n)是矩形窗，则所以，偏差(4.1.4)由此可以看出:

(1)这种自相关函数的估计是一个有偏估计，且估计的偏差是，当N→+∞时,。因此是rx(m)的渐近无偏估计。

(2)对于一个固定的N，当|m|越接近于N时，估计的偏差越大。

(3)由式(4.1.3)可看出，是真值rx(m)和三角窗函数

,,的乘积，q(m)的长度是2N－1，它是由矩形数据窗ar(n)的自相关所产生的。的方差是(4.1.5)而(4.1.6)假定x(n)是零均值的高斯随机信号，有

E[x(n)x(k)x(n+m)]=r2x(n－k)+rx(n－k－m)rx(n－k+m)

所以

将上式和式(4.1.3)代入式(4.1.5)得(4.1.7)令n－k=i，式(4.1.7)可写成

(4.1.8)在大多数情况下，rx(m)是平方可求和的，所以当N→+∞时，var[rx(m)]→0，又因为limbias[rx(m)]=0，所以对于固定的延迟|m|，rx(m)是rx(m)的一致估计。对由式(4.1.2)得到的自相关函数估计rx(m)进行傅立叶变换:^^^^(4.1.9)其中，v(m)是平滑窗，其宽度为2M+1，以此作为功率谱估计，即为BT谱估计。因为用这种方法求出的功率谱是通过自相关函数的估计间接得到的，所以此法也称为间接法。4.1.2周期图法

周期图法是把随机信号的N个观察值xN(n)直接进行傅立叶变换，得到XN(ejω)，然后取其幅值的平方，再除以N，作为对x(n)真实功率谱Px(ω)的估计。以Pper(ω)表示周期图法估计的功率谱，则^(4.1.10)其中(4.1.11)a(n)为所加的数据窗，若a(n)为矩形窗，则因为这种功率谱估计的方法是直接通过观察数据的傅立叶变换求得的，所以习惯上又称之为直接法。周期图法功率谱估计的均值为(4.1.12)(4.1.13)其中4.1.3周期图法与BT法的关系

式(4.1.9)中取M为其最大值N－1，且平滑窗v(m)为矩形窗，则

令n+m=l，上式可变成所以(4.1.14)由此可见，周期图法功率谱估计是BT法功率谱估计的一个特例，当间接法中使用的自相关函数延迟M=N－1时，二者是相同的。

1.Welch法

假定观察数据是x(n)，n=0，1，…，N－1，现将其分段，每段长度为M，段与段之间的重叠为M－…k。如图4.1所示，第i个数据段经加窗后可表示为

xiM(n)=a(n)x(n+ik)，i=0,1,…,L－1,n=0,1,…,M－1

其中，k为一整数，L为分段数，它们之间满足如下关系:

(L－1)k+M≤N

该数据段的周期图为

(4.1.15)其中(4.1.16)U为归一化因子，使用它是为了保证所得到的谱是真正功率谱的渐进无偏估计。由此得到平均周期图图4.1数据分段方法(4.1.17)如果x(n)是一个平稳随机过程，每个独立的周期图的期望值是相等的，根据式(4.1.15)和式(4.1.16)有(4.1.18)其中(4.1.19)A(ω)是对应M个点数据窗a(n)的傅立叶变换，若M值较大，则Q(ω)主瓣宽度较窄，如果Px(ω)是一慢变的谱，那么认为Px(ω)在Q(ω)的主瓣内为常数，这样式(4.1.17)可以写成(4.1.20)为了保证Welch法估计的谱是渐进无偏的，必须保证

(4.1.21)或(4.1.22)根据Parseval定理，式(4.1.22)可写成所以归一化因子U应取成(4.1.23)的方差表达式为

(4.1.24)如果x(n)是一个平稳随机过程，上式的协方差仅仅取决于i－l=r，令(4.1.25)式(4.1.24)可写成单求和表示式(4.1.26)其中表示某一数据段的周期图方差，即而是与的归一化协方差，如果

各个数据段的周期图之间的相关性很小，那么式(4.1.26)可近似写成(4.1.27)这也就是说,平均周期图的方差减小为单数据段图方差的1/L。

2.Bartlett

对应Welch法，如果段与段之间互不重叠，且数据窗选用的是矩形窗，此时得到的周期图求平均的方法即为Bartlett法。可以从上面讨论的Welch法得到Bartlett法有关计算公式，第i个数据段可表示为

xiM(M)=x(n+iM)，i=0，1,…,L－1,n=0,1,…,M－1

(4.1.28)

其中，LM≤N。该数据段的周期图为

(4.1.29)其中(4.1.30)平均周期图为(4.1.31)其数学期望为(4.1.32)其中(4.1.33)将式(4.1.33)与式(4.1.13)相比，取平均情况下A(ω)的主瓣宽度是不取平均情况下A(ω)的主瓣宽度的N/M。由此可知，取平均以后，由式(4.1.32)与式(4.1.33)计算的平均周期图偏差要比式(4.1.12)与式(4.1.13)计算的平均周期图偏差大，同时分辨率也下降。而平均周期图的方差仍可应用式(4.1.26)计算，由于数据段非重叠，各数据段的相关性比Welch法各数据段的相关性要小，因此平均周期图的方差更趋向于式(4.1.27)的理论结果，但要注意，在N一定的情况下，此时所能分的段数比Welch法有重叠情况下所能分的段数L小，因此总的来说，Welch法的计算结果要比Bartlett法好。4.1.5经典功率谱估计性能比较

为了对几种经典功率谱估计方法的性能进行比较，我们采用参考文献27中的数据。信号表示为

其中包含有四个复正弦，其归一化频率分别是f1=0.15，f2=0.16，f3=0.252，f4=－0.16,Ak对应不同的系数，可得到不同的信噪比，本数据在f1处的信噪比为64dB，在f2处的信噪比为54dB，在f3处的信噪比为2dB，在f4处的信噪比为30dB。y(n)是一个复值的噪声序列，其功率谱为其中，σ2=0.01，b(k)是模型系数。

x(n)的真实功率谱曲线如图4.2(a)所示，注意其频率范围是从－0.5~0.5，即－π~π。令归一化频率f1和f2相差0.01，目的是检验算法的分辨能力；f3的信噪比很低，目的是检验算法对弱信号的检测能力。

现取N＝128，图4.2(b)示出了该数据段直接求出的周期图，所用数据窗为矩形窗，由于主瓣过零点宽度B=2/128=0.015625>0.01，所以f1和f2不能完全分开，只是在波形的顶部能看出两个频率分量。

图4.2(c)是利用Welch平均法求出的周期图，共分四段，每段32点，没有重叠，使用Hamming窗，这时谱变得较平滑，但分辨率降低。图4.2功率谱估计方法比较图4.2(d)也是用Welch平均法求出的周期图，共分七段，每段32点，重叠16点，使用Hamming窗，谱变得更加平滑，分辨能力和图4.2(c)大体一致。

图4.2(e)是用BT法求出的功率谱曲线，M=32，没有加窗；图4.2(f)也是用BT法求出的欧尼功率谱曲线，M=16，使用了Hamming窗。

4.2AR模型功率谱估计的方法和性质

4.2.1AR模型功率谱估计的引出

由第2章讨论可知，p阶AR模型满足如下差分方程:

xAn+a1xAn－1+…+apxAn－p=εn

其中，a1，a2，…，ap为实常数，且ap≠0；εn是均值为0、方差为σ2ε的白噪声序列，也就是说，随机信号xAn可以看成是白噪声εn通过一个系统的输出，如图4.3所示。图4.3AR模型信号图4.3中

(4.2.1)而A(z)=1+a1z－1+…+apz－p

(4.2.2)已经证明(4.2.3)其中，rA(m)是AR模型的自相关函数，尤其对于0≤m≤p，由式(4.2.3)可写出矩阵方程为

(4.2.4)这就是AR模型的正则方程，又称Yule-Walker方程。对于一个p－1阶预测器，预测值为(4.2.5)其中，h(k)=－a(k)，预测误差为(4.2.6)其中，a(0)=1。p阶预测误差滤波器Ap(z)如图4.4所示。图4.4p阶预测误差滤波器图中，

Ap(z)=1+a(1)z－1+…+a(p)z－p

当E[e2(n)]达到其最小值E[e2(n)]min时，必满足Yule-Walker方程

(4.2.7)当x(n)就是图4.3所产生的p阶AR过程xAn，也即xn=xAn或rx(m)=rA(m)时，m=0，1，…，p，必满足关系式(4.2.8)或

A(z)=Ap(z)

(4.2.9)

此时，预测误差滤波器Ap(z)就是AR模型H(z)的逆滤波器，实际上也就是一个白化滤波器，而且它的输出en得到了完全的白化，也即en是一个方差为σ2ε的白噪声εn。

通过上面的分析可以看出，对于一个p阶的AR过程xn，如果首先建立阶数等于p或大于p的预测误差滤波器Ap(z)，然后以1/Ap(z)构成一个AR模型，那么以方差为E[e2(n)]min的白噪声εn通过此线性系统，其输出功率谱必定与待估计的随机信号的功率谱完全相同，因此，模型H(z)可以完全表示出AR(p)过程的xn功率谱，它们的关系即是

(4.2.10)当然，实际待估计的随机信号xn可能是一个阶数大于p的AR过程，也可能根本就不是一个AR过程，但仍可以采用上述方法建立一个p阶的AR模型，作为对随机信号xn的功率谱估计，此功率谱估计可作为xn真实功率谱的一个近似，其步骤是:

(1)对此随机信号xn建立p阶的线性预测误差滤波器，求得系数a(1)，a(2)，…，a(p)和E[e2(n)]min；

(2)令A(z)=1+a1z－1+…+apz－p，其中a1=a(1)，…，ap=a(p)，并构成一线性系统(4.2.11)那么将一方差为σ2ε的白噪声εn通过该系统，其输出的功率谱可作为待估计随机信号xn的功率谱估计(4.2.12)其中，σ2ε=E[e2(n)]min。有下述几点需要注意:

(1)预测误差滤波器Ap(z)实际上就是一个白化滤波器，但一般情况下，它不能将随机信号xn进行完全白化，所以图4.4中预测误差滤波器的输出en并不是白噪声εn,因此E[e2(n)]min表示的也是非白噪声的方差；只有当被估计的随机信号xn本身就是一个AR(p)过程,且当预测误差滤波器的阶数大于或等于p时，其输出en才是一个白噪声εn，E[e2(n)]min

即表示此白噪声的方差。

(2)对比式(4.2.4)和式(4.2.7)，当ak=a(k)，k=1，2，…，p和E[e2(n)]min时，用这样的方法所建立的p阶AR模型，其自相关函数rA(m)和待估计随机信号xn的自相关函数rx(m)必有如下关系:

(4.2.13)由式(4.2.13)可见，用AR(p)模型来对任意的随机信号建模，一定具有自相关函数的匹配性质，随着模型阶数p的增大，匹配的程度会越来越好；当p→+∞时，AR模型的自相关函数和随机信号xn的自相关函数可达到完全的匹配，即(4.2.14)而当p→+∞时，p阶FIR预测误差滤波器退化为一个IIR预测误差滤波器，它一定可将随机信号xn进行完全的白化。为了说明这个问题，我们又重画出预测误差滤波器(如图4.5所示)，图中G(z)是随机信号xn的建立模型。图4.5AR模型与预测误差滤波器当p→+∞时因此，预测误差滤波器(4.2.15)这证明了上述预测误差滤波器可将随机信号完全白化的结论。4.2.2AR模型谱估计的性质

1.隐含的自相关函数延拓的特性

在前面讨论的经典BT法功率谱估计中，假定由给定的数据

xN(n)，n=0，1，…，N－1，可估计出自相关函数rx(m)，m=－(N－1)~(N－1)，在这个区间以外，用补零的方法将其外推，对此求其傅立叶变换

^(4.2.16)就可得到BT法的功率谱估计PBT(ω)，此PBT(ω)的分辨率显然是随着信号长度N的增加而提高的。^^而在AR模型谱估计中，上述限制不再存在。虽然给定的数据xN(n)，n=0，1，…，N－1,是有限长度，但现代谱估计的一些方法，包括AR模型谱估计法，隐含着数据和自相关函数的外推，使其可能的长度超过给定的长度。前面讨论的AR模型的建立，用到了单步预测的概念，预测值为这样x(n)可能达到的长度是1~(N－1+p)，如果在递推的过程中，用x(n)代替x(n)，那么还可继续不断地外推。同样，从AR模型的建立过程看，AR过程的自相关函数必满足^^由此式可见，AR模型的自相关函数在0~p范围内与rx(m)完全匹配，而在这区间外，可用递推的方法求得。自相关函数

rA(m)实际上就是被估计信号xn的自相关函数的估计

将其进行傅立叶变换，就可得到随机信号xn的功率谱(PSD)估计，即(4.2.17)比较式(4.2.16)和(4.2.17)可见，AR模型法避免了窗函数的影响，因此它可得到高的谱分辨率，同时它所得出的功率谱估计Px(ω)与真实的功率谱Px(ω)偏差较小。图4.6示出了AR模型谱估计和BT谱估计法的比较。^^图4.6AR模型谱估计和BT谱估计性能比较2.AR模型的稳定性

AR模型稳定的充分必要条件是H(z)=1/A(z)的极点(即A(z)的零点)必须在单位圆内，而且这一条件也是保证x(n)是一个广义平稳过程所必需的。这很容易证明，如果H(z)有一个极点在单位圆外，那么x(n)的方差将趋于无穷，则x(n)是非平稳的。

因为AR模型的系数是由Yule-Walker方程解得的，可以证明如果(p+1)×(p－1)的全相关矩阵

是正定的，AR模型的系数具有非零解，此时预测误差滤波器A(z)一定具有最小相位的特性；换句话说，AR模型H(z)一定是一个稳定的全极点滤波器。使用AR模型对纯正弦信号建模是不合适的，因为此时Rp可能出现奇异或非正定的情况，但在信号处理中经常要用正弦信号作为试验信号以检验某个算法或系统的性能。为克服自相关矩阵奇异的情况，最常用的方法是加上白噪声εn，这样det(Rp)就不会等于零。

3.谱的平坦度

前面的讨论已经指出，AR(p)的系数ak就是预测误差功率最小时的p阶线性预测误差滤波器的系数，由于预测误差滤波器是一个白化滤波器，它的作用是去掉随机信号x(n)的相关性，在自己的输出端得到白噪声εn，因此在这一节中，把白化的概念加以推广，表明AR参数也可以使用预测误差滤波器Ap(z)的输出过程具有最大的谱平坦度的方法得到。利用谱平坦度的概念可以把AR谱估计得到的结果看成是最佳白化处理的结果。

功率谱密度的谱平坦度可定义为

(4.2.18)它是Px(ω)的几何均值与算术均值之比，可以证明

0≤εx≤1

如果Px(ω)有很多峰(也就是它的动态范围很大)，例如，在由p个复正弦所组成的随机信号的式子中，Ak、ωk为常量，jk是在(－π，+π)范围内均匀分布的随机变量，则此种信号的功率谱具有最大的动态范围。将上式代入式(4.2.18)，分子显见为零，因此εx=0。但如果Px(ω)是一个常数(也就是它的动态范围为零)，也即相当于xn是一个白噪声，则由式(4.2.18)显见εx=1，由此可见，谱平坦度εx直接度量了谱的平坦程度。

现设预测误差滤波器

为最小相位，输入时间序列x(n)是任意的(不一定是AR过程)，按照使输出误差序列e(n)的谱平坦度最大的准则来确定预测系数，为此，引入下述结果:如果Ap(z)是最小相位，则

(4.2.19)计算预测误差滤波器输出过程e(n)的平坦度。因为对上式两端取指数并除以，得由于对于随机信号x(n)来说，ξx和rx(0)均是固定的，因此要使ξe最大，必须使re(0)最小，因为re(0)=E[e2(n)]，因此使预测误差e(n)的功率谱平坦度最大和使p阶预测误差滤波器输出的误差功率最小是等效的，亦即条件maxξe和条件a完全等效。如果x(n)本身就是一个AR(p)过程，也即其中，A(ω)=1+a1e－jω+…+ape－jpω。现使其通过一个p阶预测误差滤波器

Ap(ω)=1+a(1)e－jω+…+a(p)e－jpω

在满足maxξe的条件下，一定有A(ω)=Ap(ω)，也即aak=a(k),

k=1，2，…，p此时预测误差滤波器输出的误差序列e(n)一定是一个白噪声序列。反之，如果将AR(p)通过一个k阶预测误差滤波器k＜p，同样在满足maxξe的条件下，误差序列e(n)不可能是一个白噪声序列，这一结果与前面讨论的AR模型谱估计的引出中所得的结果是完全一致的。a在这里有一个重要的概念需强调，即预测误差滤波器的输入、输出功率谱总满足关系式在满足E[e2(n)]min的条件下，根据求得的Ap(ω)建立AR模型:其中，σ2ε=E[e2(n)]min，是一个常数。比较这两式可见，用PA(ω)作为Px(ω)的一个估计，其估计的好坏完全取决于Pe(ω)与一个常量相逼近的程度，换句话说，在建立AR模型时，正是由于用σ2ε代替了Pe(ω)，才使得建立的AR模型功率谱PA(ω)中丧失了很多Pe(ω)的重要细节，而只有当误差序列e(n)是一个白噪声序列，Pe(ω)是一个常量，并且就等于σ2ε时，才能得到Px(ω)=PA(ω)=Px(ω)。^4.2.3AR模型参数提取方法

由前面所讨论的AR模型建立过程可知，AR模型参数必满足Yule-Walker方程:对于某随机过程来说，如果已知rx(0)、rx(1)、…、rx(p)，就可由上式求解出(σ2ε，a1，…，ap)，但实际上，对于一个所要估计的随机过程{xn}来说，通常已知的仅仅是有限的观察数据，如x(0)，x(1)，…，x(N－1)，而其自相关函数值通常是未知的，因此通常是首先求得自相关函数的估计值，然后再通过某种算法，得到AR模型参数的估计值。下面讨论常用的几种AR模型参数的提取方法。

1.自相关法

假定观察到的数据是x(0)，x(1)，…，x(N－1)，而对于无法观察到的区间(即n<0和n>N－1)，x(n)的样本假定为零，预测误差功率的表示式为

(4.2.20)式(4.2.20)的求和限又可以写为[0，N－1+p]，这是因为在此区间外误差表示式总为0的缘故。为了得到预测误差功率的极小值，将式(4.2.20)对a(l)求微分并令其等于0，得

(4.2.21)令(4.2.22)代入式(4.2.21)得(4.2.23)写成矩阵形式，这组方程为(4.2.24)求出白噪声方差σ2ε的估计值σ2ε，它是

^将式(4.2.21)代入上式，得将式(4.2.25)与式(4.2.23)合并，得(4.2.26)(4.2.25)或写成(4.2.27)例4.1

利用有限长加窗数据其中，a(n)为矩形窗，x∞(n)=rnu(n)，|r|<1，采用自相关法对其用一阶的AR模型建模。解:

首先计算rx(0)和rx(1):^^

由此得而a(1)的真正解应该是－r，这说明当观察数据为有限值N时，其解总是有偏的，直到N→+∞时，a(1)才收敛于－r。^^

2.协方差法

假定观察到的数据是x(0)，x(1)，…，x(N－1)，写出预测误差功率的表示式为(4.2.28)比较式(4.2.28)和式(4.2.20)可见，协方差法与自相关法的区别主要在于预测误差功率求和式的上下限取的不同。由于协方差法对于观察区间[0，N－1]外的x(n)样本并未假定为零，这就要求式(4.2.28)中x(n－k)总是落在观察区间[0，N－1]中，为此预测误差功率的求和上下限必须取在[p，N－1]之间。将式(4.2.28)对a(l)求微分，井令其等于0，以得到预测误差功率的极小值，得

(4.2.29)令将其代入式(4.2.29)，得写成矩阵形式，这组方程为

(4.2.30)求出白噪声方程σ2ε的估计值σ2ε:^将式(4.2.29)代入上式，得将上式与式(4.2.30)合并，得

或写成(4.2.31)比较式(4.2.31)和式(4.2.27)，表面上看，协方差法和自相关法最终所得到的正则方程具有相同的形式，但实际上两者具有完全不同的内容，主要表现在以下两个方面。式(4.2.31)中自相关矩阵Rp是对称的半正定矩阵，且不具有Toeplitz性质，故式(4.2.31)的求解不能采用Levinson递归算法来进行，虽然现在已经提出一系列协方差法的求解算法，但总的来说，其算法求解的复杂性仍远远超过自相关法。另外由于Rp的非正定性，它也不能保证所得到的预测误差滤波器具有最小相位特性，换句话说，利用协方差法估计出的AR模型极点不能保证在单位圆内。举个例子，当p=1和N＝2时，有^^而因此

显然a(1)的幅度也可能大于或等于1。采用协方差法对信号进行建模，能够较好地反映出信号真正的模型。例4.2

采用与自相关法讨论中相同的例子，根据x(n)，n=0，1，…，N－1的观察数据采用协方差法对其用一阶的AR模型建模。首先计算rx(1，0)和rx(1，1)，即^^^由此得

由此可见，只要N≥2，所建立的模型就与信号的真正模型严格相同。同理可算出由此求得这说明一阶AR模型和数据的真正模型完全一致。

3.修正的协方差法

假定观察到的数据是x(0)，x(1)，…，x(N－1)，分别写出前向预测与后向预测的表示式

和其中a(k)是AR模型系数，然后分别写出前向预测误差功率和后向预测误差功率的表示式则前向预测误差功率ρf

和后向预测误差功率ρb

的平均值ρ为

^^^(4.2.32)修正的协方差法采用前向和后向预测误差功率的平均值ρ为极小的方法估计AR参数。为使式(4.2.32)极小化，可以将ρ

相对于a(l)求微分，并令其等于0，得^^经整理得

令式(4.2.34)可写成(4.2.35)(4.2.36)l=1，2，…，p将式(4.2.36)写成矩阵形式，即求出白噪声方差σ2ε的估计值σ2ε:^其中利用了式(4.2.36)的结果，最后得

将上式与式(4.2.36)合并得(4.2.37)或写成修正的协方差法的一个优点是，它的误差功率的计算(见式(4.2.32))是在相对于协方差法多一倍的数据点上进行，这在观察数据长度很短的情况下，是非常有利的，但是这要求信号在正反两个方向上呈现相同的特性。例如，噪声中叠加上多个正弦信号就具有这种特性，因为正弦信号在正反两个方向上看起来是相同的，而噪声的自相关函数也是与方向无关的。另一方面，若对于一个在两个方向呈现不同特性的信号，如一个指数衰减的信号，从相反的方向看，就是一个指数增长的信号，对这种信号如仍采用修正的协方差法来进行AR模型的估计，可能就会得到不好的结果。例4.3

采用与自相关法和协方差法讨论中相同的例子。根据x(n)，n=0，1，…，N－1的观察数据，采用修正的协方差法对其用一阶的AR模型建模。

解:

首先计算rx(1，0)和rx(1，1)，即

^^和由此得因为模型的真正解是a(1)<－r，这说明此模型参量的估计是有偏的，且与数据的长度N无关。因为r<1，所以|a(1)|>|a(1)|，因此一阶AR模型的极点比真正的极点更接近于单位圆。换句话说，用修正协方差法对一随机信号进行建模，可以得到一“锐化”的谱分析结果，因此用修正协方差法对随机信号建模，可以很好地用来进行高分辨率谱估计。

式(4.2.37)中自相关矩阵Rp是一个非Toeplitz矩阵，另外，虽然Rp是正定的(纯正弦信号情况除外)，但它也不能保证所得到的预测误差滤波器具有最小相位特性，即利用修正协方差法所估计出的极点与协方差法一样，不能保证在单位圆内。但是有一个例外，在一阶的情况下，^^^p=1时用修正协方差法所建立的AR模型又总是稳定的，换句话说，所估计出的极点总是在单位圆内。这是因为在p=1时，由式(4.2.35)可知和由此得因为在观察数据x(n)不恒为0的情况下，上式的分母总是大于分子，因此|a(1)|<1。^4.2.4AR模型阶次的选择

在AR谱估计中模型阶次的选择是一个关键问题，正如前面所讨论的：阶次太低将会导致平滑的谱估计结果，如图4.7所示；而阶次太高，将会产生虚假谱峰，并且估计的方差也会增大。

一种简单而直观的模型阶次估计方法是基于预测误差功率，因为对于所有讨论过的AR模型参数估计方法，预测误差功率都是随模型阶次的增加而减小，但是我们不能简单地只把监测预测误差功率的减小作为确定模型阶次的方法，还必须考虑模型阶次的增加会导致谱估计方差的增大。图4.7AR谱估计值与阶次的关系

1.最终预测误差准则

最终预测误差(FPE)准则是通过使下式达到最小来估计模型阶次的，即(4.2.39)其中，ρk是k阶AR模型的预测误差功率估计值，N是观察数据的长度。可以看出，尽管ρk随k的增大而减小，但(N+k)/(N－k)却随k的增大而增大，它表示了由于AR模型参数估计的不精确而作出的对预测误差功率的一个修正，所以FPE(k)将有一个最小值，这个FPE(k)的最小值所对应的阶便是最后确定的阶。该准则的实际应用表明，虽然对于AR过程来说效果很好，但在处理地球物理数据时，一般都认为这一准则确定的阶次偏低。^^

2.阿凯克信息准则

阿凯克(Akaike)信息准则(AIC)是通过使下式达到最小来估计模型阶次的，即

AIC(k)=Nlnρk+2k

(4.2.40)

可以证明AIC表示的是AR模型估计的PDF和数据的真实PDF之间的库尔贝克-利布勒(Kullback-Leibler)距离的估计值，这个准则不仅可用来确定AR模型的阶次，还可用于MA模型和ARMA模型阶次的确定。

AIC和FPE的性能是相近的，建议对于短观察数据使用AIC。对于长观察数据(N→+∞),这两种估计方法将得到相同的模型阶次估计。^

3.自回归传递函数准则

自回归传递函数准则(CAT)是通过使下式达到最小来估计模型阶次的，即(4.2.41)这里其中，ρi是i阶AR模型的预测误差功率估计值。CAT准则选择AR模型阶次，使得由该模型估计出的预测误差滤波器最接近于最佳无限长滤波器。^

4.3最大熵谱估计方法

给定一个宽平稳过程的自相关rx(k)在|k|≤p时的值，要解决的问题是如何外推出k>p

时的rx(k)值。用re(k)表示外推的值，显然应对re(k)加一些约束。例如，若(4.3.1)则Px(ejw)应对应于一个合法的功率谱，即Px(ejw)对所有w都是实值的且非负的。但一般来说,只约束Px(ejw)是实值和非负的还不能保证获得唯一的外推结果，因此必须对可允许的外推再加上一些约束。Burg就提出一个这样的约束，就是使随机过程的熵达到最大化。由于熵是随机性或不确定性的一个测度，最大熵外推等价于要找出使x(n)尽可能“白”(随机)的自相关序列re(k)。在某种意义上，这种外推对x(n)加了尽可能少的约束或最少量的结构要求。对功率谱而言，这对应于约束Px(ejw)“尽可能平坦”。如图4.8所示，不同的外推方法会得到不同的功率谱形状，其中最上面的外推方法有最“平坦”的谱。图4.8自相关序列的不同外推方法及其对应功率谱对一个功率谱为Px(ejw)的高斯随机过程，其熵是

(4.3.2)因此对于高斯过程，若已知部分自相关序列rx(k)(|k|≤p)，其最大熵功率谱是使式(4.3.2)最大化,同时约束Px(ejw)的逆DTFT在|k|≤p时等于给定的自相关值，即(4.3.3)为求使熵最大化的re(k)值，取H(x)对r＊e(k)的导数并使之等于0，即(4.3.4)由式(4.3.1)可得，再代入上式就有

(4.3.5)定义Qx(ejw)=1/Px(ejw)，则式(4.3.5)表明Qx(ejw)的逆DTFT应是一个有限长序列qx(k),在|k|>p时值为0，即因此且高斯过程的最大熵功率谱Pmem(ejw)是一个全极点功率谱，为

^(4.3.6)利用谱因子化定理，式(4.3.6)可以表达为若用矢量ap=[1，ap(1)，…，ap(1)]和e=[1，ejw，…，ejpw]T来表示，则MEM谱可写为(4.3.7)确定了MEM谱的形式后，剩下的是求系数ap(k)和b(0)。由于式(4.3.3)给定的约束，必须选择这些系数使得Pmem(ejw)的逆DTFT产生的自相关序列与给定的rx(k)在|k|≤p时的值相匹配。在4.2节已看到，若系数ap(k)是如下自相关正则方程的解:^(4.3.8)且若则所获得的自相关将满足式(4.3.6)的约束。因此，MEM谱估计是(4.3.9)其中，ap是方程(4.3.8)的解。显然，这时的MEM谱估计就是AR模型谱估计中自相关法的结果，只是表达为上式的形式。当然，若对信号的类型或熵的表述采用不同的定义，则最大熵谱估计还有其他的结果，这里不再深入讨论。例4.4(噪声中复指数的MEM估计)考虑一个白噪声中含有复指数的随机过程

x(n)=A1ejnw1+w(n)

式中A1=|A1|ejj，j是在[－π，π]区间均匀分布的随机变量，白噪声w(n)的方差为σ2w，求其p阶MEM谱。

解:先解自相关正则方程以获得AR系数ap(k)。x(n)的(p+1)×(p+1)阶自相关阵是

Rx=P1e1eH1+σ2wI

其中

e1=[1，ejw1，…，ejpw1]T

P1=|A1|2Rx的逆阵是

因此(4.3.11)其中，u1=[1，0，…，0]T。

MEM谱估计应为由于

其中，WR(ejw)是矩形窗的DTFT，因此MEM谱变成为确定εp值，注意由于ap(0)=1，因此由式(4.3.11)得解上式得MEM谱又变成

与最小方差法的结果一样，MEM谱的峰值发生在频率w=w1处，且若信噪比很大，即P1>>σ2w，则MEM谱估计在w=w1处的值近似为因此MEM谱的峰值正比于复指数分量的平方。

4.4最大似然谱估计

当信号形式为y(t)=Aejω0t时，第k个传感器的接收信号为

zk(t)=yk(t+τk)+vk(t)，1≤k≤p

(4.4.1)

其中，τk是在第k个传感器上的传输延迟，vk(t)为在第k个传感器上的白噪声过程。传输延迟为其中，·表示向量的点积，zk为第k个传感器的位置，k0是被传输的信号的方向，c是光速。对于一个均匀线性传感器阵列(见图4.9)，有

zk=kdlx式中，lx为沿x轴的单位向量，d是传感器之间的距离。图4.9传感器阵列由图4.9可知

zk·k0=kdsinθ0

其中，θ0为接收信号的到达角。该传感器的输出延迟τm秒后，再乘以权系数am，将乘积项求和，可得如下形式的输出:

(4.4.2)第m个传感器的输出为zm(t－τm)=y(t)+vm(t－τm)

(4.4.3)且有若vm(t)，m=(1，…，p)是互不相关的均值为0、方差为σ2v的随机过程，则均方功率为(σ2y为信号功率)

E[x2(t)]=p2σ2y+p2σ2v

(4.4.4)

因此在理想情况下，输出信噪比是任一传感器上的信噪比的p倍。波束形成的目的之一是通过调节延迟来收集来自该信号源的信号能量，而拒绝所有其他信号(包括噪声在内)。

若定义A=[a1，…，ap]T，Z=[zn，…，zn－p]T，则其采样输出为

其方差为式中，R=E[ZZT]，且Δt为采样间隔。

在高斯环境下最大似然估计与最小方差估计一致，所以对x(t)的估计即为

(4.4.5)

由于为常数，故应使，即:

使方差在下列约束条件下为最小:

式中，ED=[1，D1，D2，…，DP]，D=ejωΔt。

最小化可减小来自k0以外的其他方向的能量，约束条件使接收器的增益为1。对于正弦输入，频率ω0处的正弦输出信号的增益为1。引入拉格朗日乘子，构造目标函数:

(4.4.6)其中，λ是拉格朗日乘子。令则可求出其中，EHD=(E*D)T。代入方差表达式得当向量ED中的频率ω等于信号频率时，上式就代表信号功率。当频率ω的值遍历定义域时，上式就可以看做是过程x(t)的功率谱密度的估值。这样就得到了最大似然谱估计的值，为(4.4.7)可见，自相关函数矩阵R→ΦML(ω)。

可以证明:最大似然谱估计ΦML(ω)和最大熵谱估计

ΦME(ω)之间有密切联系，即:最大似然谱估计的倒数等于1~p阶的最大熵谱估计的倒数和。^^^式中，为m阶AR模型(最大熵)谱估计。通常，先用自回归AR模型法计算，m=1，2，…，p,再用上式计算最大似然谱估计。一般来说，最大似然谱估计的分辨率要比最大熵谱估计的分辨率低。

4.5互协方差估计与互谱估计

假设x(n)和y(n)都是零均值的随机信号，因此γxy(m)=

φxy(m)。与自协方差的估计相同，互协方差估计为

(4.5.1)当y(n)=x(n)时，上式就完全转化为自协方差估计式。式(4.5.1)的期望值为其中，jxy(m)是随机信号x(n)和y(n)的互相关。同理

合并上述二式，得(4.5.2)可见，γxy(m)是协方差jxy(m)的一个渐近无偏估计，并且估计的方差与N成反比。为了估计互功率谱，取γxy(m)的傅立叶变换得到谱估计^^(4.5.3)若y(n)=x(n)，上式显然转化为周期图。因为γxy(m)一般没有对称特性，所以Γxy(ω)一般是复函数。容易证明^^显然从上式可知，Γxy(ω)是互功率谱的渐近无偏估计，但是完全与周期图中的情况一样，随着N的增加，Γxy(ω)的方差并不趋向于零。为了把方差减小并把估计平滑，必须对根据较短的记录段做出的估计进行平均或加窗处理。^^

4.6特征分解法谱估计

前面我们介绍了基于模型的宽平稳随机过程的功率谱估计问题，信号可以建模为白噪声激励的线性移不变滤波器的输出。另一个重要的模型是x(n)为白噪声与多个复指数分量之和，即

(4.6.1)这里假设幅度Ai是复数，有其中,fi是互不相关的且在区间[-π,π]内均匀分布的随机变量。4.6.1自相关阵的特征分解

为理解自相关阵的特征分解可作为频率估计的一种方法，先来看一个一阶谐波过程:

x(n)=A1ejnw1+w(n)

它表示白噪声中有一个复指数分量，复指数的幅度是Ai=|Ai|ejf1，f1是均匀分布的随机变量，w(n)是方差为σ2w的白噪声。前面已介绍过，该x(n)的自相关序列是

rx(k)=P1ejkw1+σ2wδ(k)

其中，P1=|A1|2是复指数分量的功率。因此，x(n)的M×M阶自相关阵是信号的自相关阵Rs和噪声的自相关阵Rn的和，即

Rx=Rs+Rn

(4.6.2)其中信号自相关阵的具体形式是(4.6.3)其秩为1。噪声自相关阵是个对角阵Rn=σ2wI，它是满秩的。注意，若我们定义

e1=[1，ejw1，ej2w1，…，ej(M－1)w1]T

(4.6.4)则Rs可用e1表示为Rs=P1e1eH1由于Rs的秩为1，因此Rs只有一个非零特征值。又由

Rse1=P1(e1eH1)e1=P1e1(eH1e1)=MP1

得Rs的非零特征值应等于MP1，且e1就是对应的特征矢量。另外，由于Rs是共轭对称阵，因此剩下的特征矢量v2,v3，…，vM应与e1正交(共轭对称阵的不同特征值对应的特征矢量是相互正交的)，即

e1Hvi=0，i=2，3，…，M

(4.6.5)

最后，若用λsi表示Rs的特征值，则有

Rsvi=(Rs+σ2wI)vi=λsivi+σ2wvi=(λsi+σ2w)vi

(4.6.6)

这说明Rx的特征矢量和Rs的特征矢量是相同的，只是Rx的特征值是

λi=λsi+σ2w其中，Rx的最大特征值是λmax=MP1+σ2w，剩下的Rx的特征值都等于σ2w。这样，就可以由Rx的特征值和特征矢量获得所有的有关x(n)的参数，具体方法如下:

(1)对自相关阵Rx进行特征分解，其最大特征值应等于MP1+σ2w，剩下的特征值都是σ2w。

(2)利用Rx的特征值求信号功率P1和噪声方差σ2w得

(3)用对应于最大特征值的特征矢量vmax确定频率w1，例如，利用vmax的第二个系数vmax(1)，应有

例4.5

设x(n)是一个一阶谐波过程，为白噪声加一个复指数信号x(n)=A1ejnw1+w(n),设其自相关阵为,求x(n)中复指数信号的频率。

解:

由于x(n)的自相关阵Rx的特征值为

特征矢量为因此白噪声的方差是复指数分量的功率是

最后得复指数分量的频率为这里要指出的是，对一阶谐波过程，实际上可以更容易地求出该过程的参数。由于

，所以立即可得，w1=π/4，而有了P1，白噪声的方差就可以确定如下:实际上，以上介绍的方法在频率估计中的实用价值是有限的，因为它要求精确地已知自相关矩阵。虽然也可以用估计的自相关阵代替Rx，但这样做后，其最大特征值只是近似等于MP1+σ2w，对应的特征矢量也只是近似为e1，而特征值和特征矢量可能对rx(k)的很小的估计误差就非常敏感。因此一般不用一个特征矢量来估计复指数的频率，而利用的是它们的加权平均。设vi表示Rx的噪声特征矢量，即特征值为

σ2w的特征矢量，vi(k)表示vi的第k个分量。如果计算vi各系数的DTFT为(4.6.7)则式(4.6.5)所给定的正交性条件就意味着Vi(ejw)在w=w1处(即复指数的频率处)将等于0，因此，若我们构造的频率估计函数如下:

(4.6.8)则Pi(ejw)在w=w1处的值将很大(理论上是无穷大)。这样，通过对频率估计函数的峰值进行定位，就可以估计出复指数的频率。但是，式(4.6.8)只利用了一个特征矢量，因此可能对Rx的估计误差很敏感，所以可以考虑利用所有噪声特征矢量的加权平均，即^(4.6.9)其中，ai是适当选择的常数。

现在来考虑白噪声中有两个复指数分量的情况，即

x(n)=A1ejnw1+A2ejnw2+w(n)

其中，Ai=|Ai|ejj1(i=1，2)，是复指数分量的幅值，w1和w2是其频率，且w1≠w2。若w(n)是方差为σ2w的白噪声，则x(n)的自相关序列是

其中，P1=|A1|2，P2=|A2|2。显然其自相关阵可写成如下几项的求和:

其中，Rs=P1e1eH1+P2e2eH2代表的是Rx中的信号分量，且Rs的秩为2；而Rn=σ2wI是Rx中的噪声分量，Rn是一个对角阵。Rx的更简洁的分解形式是

Rx=EPEH+σ2wI

(4.6.10)

其中，E=[e1，e2]是M×2的矩阵，且由两个信号矢量e1和e2组成，P=diag[P1，P2]是信号功率值组成的对角阵。

除了像式(4.6.10)那样将Rx分解为两个自相关阵之和外，也可以对Rx进行特征分解。设vi和λi分别表示Rx的特征矢量和特征值，且特征值按降序排列为

λ1≥λ2≥…≥λM

由于Rx=Rs+σ2wI，因而有

λi=λsi+σ2w

其中，λsi是Rs的特征值。由于Rs的秩等于2，所以Rs只有两个非零特征值，且都大于零(Rs是非负定的)，因此Rx的前两个特征值都大于σ2w，其他剩下的特征值等于σ2w。这样Rx的特征值和特征矢量可以分为两组:第一组是两个特征值大于σ2w的特征矢量，称为信号特征矢量，张成一个二维的子空间，称为信号子空间；第二组是那些特征值等于σ2w的特征矢量，称为噪声特征矢量，张成一个M－2维的子空间，称为噪声子空间。注意Rx是共轭对称矩阵，所以其特征矢量vi形成一个正交集合，这样，信号和噪声子空间就应是相互正交的，也就是说，对信号子空间中的任意矢量u和噪声子空间的任意矢量v，有uHv=0。与单个复指数分量时的情况不同，对两个复指数分量加上白噪声的情况，信号特征矢量一般不等于e1和e2(因为信号矢量一般不是正交的，而特征矢量总是相互正交的)，但e1和e2将位于信号特征矢量v1和v2所张成的子空间中，且由于信号子空间和噪声子空间相互正交，所以e1和e2与噪声特征矢量vi是正交的，即

(4.6.11)这样，与单个复指数分量一样，复指数的频率w

1和w

2可以用如下形式的频率估计函数来进行估计:最后来看更一般的情况，即p个不同频率的复指数加上白噪声构成一个宽平稳过程，其自相关序列应是其中，Pi=|Ai|2是第i个分量的功率。它的M×M的自相关阵是(4.6.12)其中，ei=[1，ejwi，ej2wi，…，ej(M－1)wi]T(i=1，2，…，p)是p个线性独立的矢量。就像前面一样，式(4.6.12)可以更简洁地表达为

Rx=EPEH+σ2wI

(4.6.13)其中，E=[e1…ep]是p个信号矢量ei组成的M×p的矩阵，P=diag[P1，…，Pp]是信号功率值组成的对角阵。由于Rx=

Rs+σ2wI(λsi是Rs的特征值)，且Rs的秩是p，因此Rx的前p个特征值将大于σ2w，其他剩下的特征值等于σ2w。这样Rx的特征值和特征矢量可以分为两组:第一组是特征值大于σ2w的信号特征矢量v1，…，vp，另一组是特征值等于σ2w的噪声特征矢量vp+1，…，vM。假设特征矢量都归一化为单位模，则Rx可以分解为若写成矩阵形式，则有

Rx=VsVssVHs+VnVnnVHn

(4.6.14)

其中，Vs=[v1…vp]是信号特征矢量组成的M×p矩阵，Vn=

[vp+1…vM]是噪声特征矢量组成的M×(M－p)矩阵，Vss和

Vnn分别是特征值λi=λsi+σ2w和λi=σ2w组成的对角阵。

后面我们将要把一个矢量投影到信号子空间或噪声子空间，完成这种投影的投影矩阵Ps和Pn应是

Ps=VsVHs，Pn=VnVHn

(4.6.15)

与前面一样，信号子空间和噪声子空间的正交性可用于估计复指数的频率。具体来说，由于每个信号矢量e1…ep都属于信号子空间，所有的ei与每个噪声特征矢量都是正交的，即

eHivk=0，i=1，2，…，p，k=p+1，p+2，…，M

结果频率可用如下形式的频率估计函数进行估计:(4.6.16)下面各节将基于上式给出几种不同类型的频率估计算法。首先介绍Pisarenko谐波分解，它就采用式(4.6.16)的频率估计函数进行估计，并取M=p+1，aM=1。4.6.2Pisarenko谐波分解

Pisarenko分解中，假定x(n)是由p个复指数加上白噪声组成，且复指数的个数p是已知的。又假设已知p+1个自相关序列值，或它们已被估计出来，因此有一个(p+1)×(p+1)的自相关阵，噪声子空间的维数是1，它由最小特征值

λmin=σ2w所对应的特征矢量张成。用vmin表示该噪声特征矢量，则vmin将与每个信号矢量ei正交，即

(4.6.17)因此，

在每个复指数频率wi(i=1，2，…，p)处都等于零。这样，噪声特征矢量的Z变换(称为特征滤波器)就在单位圆上有p个零点，即

(4.6.18)而且由特征滤波器的根就可以获得各复指数的频率。若不想求Vmin(z)的根，也可以构造如下的频率估计函数:它是式(4.6.16)在M=p+1，ap+1=1时的特例。由于PPHD(ejw)在复指数的频率处其值很大(理论上是无穷大)，对PPHD(ejw)进行峰值定位就可得到频率的估计。虽然PPHD(ejw)是功率谱的形式，但它称为伪谱(有时称为特征谱)，因为它并没有包含复指数的任何功率信息，也不包含噪声的功率信息。^^^一旦确定了复指数的频率以后，就可以由Rx的特征值求功率谱Pi了。假设信号子空间的特征矢量v1，…，vp已被归一化，即vivHi=1，由于

Rxvi=λivi，i=1，2…，p

将上式两边左乘vHi，则得

vHiRxvi=λivHivi=λi，

i=1，2…，p

(4.6.19)将式(4.6.12)的Rx表达式代入上式，有简化后，得(4.6.20)注意，求和项中的|eHkvi|2是频率wk处信号子空间特征矢量vi的DTFT的平方幅度，即有

|eHkvi|2=|Vi(ejwk)|2

其中，。因此式(4.6.20)可以写成

上式代表p个线性方程，它有p个未知量Pk，写成矩阵形式为(4.6.21)求解该方程组就得到功率Pk。

例4.6(用Pisarenko法在噪声中估计两个复指数分量)假设有两个复指数加白噪声的随机过程，其自相关序列的前三个值是

rx(0)=6

rx(1)=1.92705+j4.58522

rx(2)=－3.42705+j3.49541

试用Pisarenko法求复指数分量的频率和功率。

解:

由于p=2，必须对3×3的自相关阵Rx进行谐波分解。首先有

其特征矢量为

特征值为λ1=15.8951，λ2=1.1049，λ3=1.0000因此，最小特征值是λ3=1.0000，对应的特征矢量是于是得特征滤波器Vmin(z)，求得其根为z1=0.5+j0.8660=ejπ/3，z2=0.3099+j0.9511=ej2π/5因此两个复指数的频率是为求复指数的功率，必须计算信号特征矢量v1和v2的DTFT的平方幅值在w1和w2处的值，有|V1(ejw1)|2=2.9685，|V1(ejw2)|2=2.9861|V2(ejw1)|2=0.0315，|V2(ejw2)|2=0.0139因此式(4.6.21)就变成其中，σ2w=λmin=1。求解上式得P1=2，P2=3。例4.7(用Pisarenko法估计单个正弦信号)设x(n)是随机相位正弦波加白噪声形成的随机过程:

x(n)=Asin(nw0+φ)+w(n)

已知rx(0)=2.2，rx(1)=1.3，rx(2)=0.8，则其3×3的自相关矩阵是

Rx=Toep{2.2，1.3，0.8}

可计算出其特征值为

λ1=4.4815，λ2=1.4，λ3=0.7185

因此白噪声功率是σ2w=λmin=0.7185。Rx的特征矢量是

注意，各特征矢量都有对称性。为估计正弦波的频率，应求特征滤波器Vmin(z)的根。Vmin(z)是v3的Z变换,即

Vmin(z)=0.4437(1－1.755z－1+z－2)

可求出Vmin(z)的根为z=e±jw0，其中w0满足

2cosw0=1.755，即w0=0.159π

最后用式(4.6.21)估计正弦波的功率。由于x(n)只包含一个正弦波，实际上可以简化计算，不需求解方程。由于该随机过程的自相关序列是因此，，而已知rx(0)=2.2，σ2w=0.7185，因此有A2=2.963。虽然Pisarenko谐波分解在数学上很精致，但它并不实用，原因之一是必须已知复指数的个数。在自相关值精确已知的理想情况下，可以先估计复指数的个数，然后检验最小特征值的重复性，就可以确定复指数的个数。但自相关是估计值时，就不便采用重复性检验准则了，因为很少会出现两个最小特征值相等的情况。Pisarenko分解的另一个局限性是它假定加性噪声是白噪声，但实际中一般都不是这样，因此会导致频率估计的偏差。虽然可以修正Pisarenko方法使其适用于非白噪声的情况，但要求已知加性噪声的功率谱。从计算量的角度来看，Pisarenko谐波分解要求出信号自相关矩阵的最小特征值和特征矢量，对高阶问题，这需要很大的计算量，但可以用基于Levinson-Durbin递归的迭代算法来求最小特征值和特征矢量，以提高计算效率。4.6.3MUSIC算法

假设随机过程x(n)有p个复指数分量和方差为σ2w的白噪声，Rx是其M×M的自相关矩阵，且M