数字信号处理-第四章

第四章功率谱估计4.1引言4.2经典谱估计4.3现代谱估计中的参数建模4.4

AR模型谱估计的性质4.5

AR谱估计的方法4.6最大熵谱估计与最大似然谱估计4.7特征分解法谱估计NO.24.1引言对信号和系统进行分析研究、处理有两类方法：一类是在时域进行，我们前面学习的维纳滤波、卡尔曼滤波和自适应滤波都属于这种方法；本章则是在频率域进行研究的一类方法。这两类方法都是信号处理的重要方法。对确定性信号傅里叶变换是在频率域分析研究的理论基础，但对于随机信号，其傅里叶变换不存在，因此转向研究它的功率谱。谱，就是信号的某些特征在频域随频率的分布。功率谱反映了随机信号功率的分布特性，有着很广泛的应用：在雷达信号处理中，回波信号的频谱提供了运动目标的位置、强度、速度等信息；在声纳系统中，为了寻找舰艇或潜艇也要对混有噪声的信号进行谱分析；在语音处理中，谱分析用来探测语音语调共振峰；在电子战中，还利用频谱对目标进行分类。按照Weiner-Khintchine定理，信号的功率谱和其自相关函数服从傅里叶变换关系：实际所能得到的随机信号的长度总是有限的，用有限长度的信号计算得到的功率谱只是真实功率谱的一个估计。功率谱估计分为经典谱估计和现代谱估计两大类。上式称做功率谱的定义，对于平稳随机信号，服从各态历经定理，集合平均可以用时间平均代替，由上式还可以推出功率谱的另一个定义：将计算自相关函数中的集合平均用时间平均代替：代入功率谱定义式，得令l=n+m,则

上式中x(n)是观测数据，Pxx(ejω)是随机变量，必须对Pxx(ejω)取统计平均值，得到该式被认为是功率谱的另一定义，周期图法谱估计

Weiner-khintchine定理表明功率谱是无限多个自相关函数的函数，但观测数据只有有限个,只能得到有限个自相关函数。按照上面的定义公式求功率谱,也需要无限多个观测数据。因此根据有限个样本数据，计算随机信号的功率谱，是一个功率谱的估计问题现代谱估计是以信号模型为基础，下图表示x(n)的信号模型，输入白噪声w(n)的均值为0，方差为σ2w，x(n)的功率谱为：如果由观测数据能够估计出信号模型的参数，则信号的功率谱可以按照上式计算出来，这样,估计功率谱的问题就变成了估计信号模型参数的问题。信号模型有很多种，如AR模型、MA模型等等，针对不同的情况，需要选择不同的模型。现代谱估计的质量比经典谱估计的质量有很大的提高。但遗憾的是，尚无任何理论能指导选择一个合适的模型，只能根据功率谱的一些先验知识，或者说一些重要的谱特性，来选择模型。4.2经典谱估计BT法是先估计自相关函数,然后进行傅里叶变换得到功率谱。设对随机信号x(n)，只观测到一段样本数据，n=0,1,2,…,N-1。根据这一段样本数据估计自相关函数,有两种估计方法，即有偏自相关函数估计和无偏自相关函数估计。有偏自相关函数估计的误差相对较小，这种估计是一种渐近一致估计：4.2.1

BT法对上式进行傅里叶变换，得到BT法的功率估计:为了减少谱估计的方差，经常用窗函数w(m)对自相关函数进行加权，此时谱估计公式为式中-(M-1)≤m≤(M-1)其它为了采用FFT计算傅里叶变换，必须将求和域(-M+1,M-1)移到(0~L-1)，功率谱的计算公式为：上式也被称为加权协方差谱估计。它要求加窗后的功率谱仍是非负的，这样窗函数w(m)的选择必须满足一个原则，即它的傅里叶变换必须是非负的，例如巴特利特窗就满足这一条件。

k=0,1,2,…,L-10≤m≤M-1M≤m≤L-ML-M+1≤m≤L-1按照有偏自相关函数公式估计自相关函数，已经证明这是渐近一致估计，但经过傅里叶变换后得到功率谱的估计，功率谱估计却不一定仍是渐近一致估计，可以证明它是非一致估计，是一种不好的估计方法。BT法中用有偏自相关函数进行估计时，它和用周期图法估计功率谱是等价的，因此BT法的估计质量和周期图法的估计质量是一样的。将功率谱的另一定义式重写如下:如果忽略上式中求统计平均的运算，假设观测数据为:x(n)

0≤n≤N-1，便得到周期图法的定义:4.2.2周期图法用周期图法计算功率谱框图由周期图法功率谱估计公式推导它与BT法的等价关系令

m=k-n,即k=m+n1.周期图与BT法的等价关系方括号中的部分是有偏自相关函数的计算公式：利用有偏自相关函数的BT法和周期图法是等价的关系已知自相关函数的估计值，m=-(N-1),…-1,0,1,…,N-1，按照BT法求功率谱的统计平均值：有偏自相关函数的统计平均值在第一章中已确定，将结果代入上式,得到2.周期图法谱估计质量分析1)周期图的偏移式中两序列乘积的傅里叶变换，其频域服从卷积关系WB(ejω)称为三角窗的谱函数。上式表明，周期图的统计平均值等于它的真值与三角窗函数频谱的卷积，因此周期图是有偏估计，但当N→∞时，wB(m)→1，三角窗函数的频谱趋近于δ函数，周期图的统计平均值趋于它的真值，因此周期图属于渐近无偏估计。周期图的方差的精确表达式很繁冗，为分析简单：假设x(n)是实的零均值的正态白噪声信号，方差是σx2，即功率谱是常数σx2，其周期图用IN(ω)表示，N表示观测数据的长度。按照周期图的定义:推导得周期图的方差：2)周期图的方差当N趋于无限大时，周期图的方差并不趋于0，而趋于功率谱真值的平方，即所以无论怎样选择N，周期图的方差总是和σ4x同一个数量级。信号功率谱的真值是σ2x，这说明周期图的方差很大。用这种方法估计的功率谱在σ2x附近起伏很大，所以周期图是非一致估计，是一种很差的功率谱估计方法。图4.2.2白噪声的周期图可以看到，随着N的增大，功率谱曲线并没有趋向真值，而是起伏变的越来越剧烈。周期图法估计功率谱不是一致估计，均方误差很大，使估计出的功率谱不可靠。其频率分辨力低是根本缺点，原因在于观测数据只有一段。由于BT法和周期图法具有等效的功率谱估计质量，因此BT法也不是一致估计，分辨率低，估计误差大。4.2.3经典谱估计方法的改进基本思想：对随机变量进行观测，得到L组独立的记录数据，用每一组数据求其周期图，然后将L个周期图加起来求平均。这样得到的周期图，其方差将是用一组数据得到的周期图的方差的1/L。1.平均周期图法假设随机信号x(n)的观测区间为：0≤n≤M-1，共进行了L次独立观测，得到L组记录数据，某一组记录数据用xi(n),i=1,2,…,L表示，第i组的周期图为：将得到的L个周期图进行平均，作为信号x(n)的功率谱估计：为了分析偏移，对上式求统计平均:其中，窗函数的频谱为表明平均周期图仍然是有偏估计，偏移和每一段的数据个数M有关；由于M<N,平均周期图的偏移比周期图的偏移大，表现在三角谱窗主瓣的宽度比周期图主瓣的宽度宽。所以其分辨率更加降低，因此也可以说，偏移的大小反映分辨率的低与高。对平均周期图求方差，由于是L次独立观测，L个周期图相互独立，因此平均周期图的方差为即平均周期图的估计方差是周期图的方差的1/L。即：是以分辨率的降低为代价换取了估计方差的减少平均周期图法这种方法是用一适当的窗函数W(ejω)与周期图进行卷积，来达到使周期图平滑的目的。式中是有偏自相关函数-(M-1)≤n≤M-12.窗函数法周期图的窗函数法和前面提到的BT法的加权协方差谱估计是类似的。窗函数法中，周期图和窗函数的频谱卷积得到功率谱，等效于在频域对周期图进行修正，使周期图通过一个线性系统，滤除掉周期图中的快变成分，谱窗函数需具有低通特性。对上式求统计平均,得到有偏自相关函数的期望等效于真值与三角函数的乘积上式表明,周期图的窗函数法仍然是有偏估计,其偏移和wB(m)、w(m)两个窗函数有关，如果w(m)窗的宽度比较窄，当M比N小得多时，|M|<<N，则wB(m)~1,上式可以近似写成：由于窗函数w(n)比wB(n)宽度窄，所以其频谱的主瓣更宽，利用窗函数法可以平滑周期图，减少估计误差，但是偏移加大了，使分辨率降低。这种方法和平均周期图法一样，首先把数据长度为N的信号x(n)分成L段，每一段数据长度为M，N=LM。然后把窗函数w(n)加到每一个数据段上，求出每一段的周期图，形成修正的周期图，再对每一个修正的周期图进行平均。第i段的修正周期图为3.修正的周期图求平均法同样,将每一段的修正的周期图之间近似看成互不相关,最后功率谱估计为这种在计算周期图之前，先对各数据段加窗函数的方法，使平均周期图的估计方差减少，当然分辨率同样变低。但这种方法对窗函数没有限制；此外，分段时，相邻的两段可以有重叠，进一步使方差减少，可以重叠50%。总之，传统的功率谱估计方法无论采取哪一种改进方法，总是以减少分辨率为代价来换取估计方差的减少，提高分辨率的问题无法根本解决。对于由白噪声和正弦信号或者窄带信号组成的随机信号，在计算周期图之前，一般应该给数据加窗，如果不加数据窗，相近的低电平信号可能被高电平信号的旁瓣淹没掉。在下图a)中，ω/π=0.12处的正弦信号的旁瓣几乎掩盖了在ω/π=0的信号，经过给数据加窗处理后，大大压低了旁瓣，使低电平信号清晰可见，如图b)所示。但由于主瓣加宽，功率谱波峰变宽了，降低了信号的分辨率。一般来说，两个等幅的正弦信号的频率相隔很近，可以不加数据窗，频率间隔应该大于2π/N才能分辨。利用数据窗减少窄带过程周期图的旁瓣a)没有数据窗b)加哈明数据窗4.3现代谱估计中的参数建模在第一章中学习了用参数模型来描述随机信号的方法。如果能根据观测数据求出信号的模型(即H(z))，假设模型的系统函数为H(z)，输入白噪声方差为σw2，则信号的功率谱就表示为：功率谱估计可分成三个步骤：1)选择合适的信号模型；2)根据x(n)有限的观测数据，或者它的有限个自相关函数估计值，估计模型的参数;3)由模型的参数计算功率谱。选择模型主要考虑的是模型能够表示谱峰、谱谷和滚降的能力。对于具有尖峰的谱，应该选用具有极点的模型，如AR和ARMA模型；对于具有平坦的谱峰和深谷的信号，可以选用MA模型；既有极点也有零点的谱应选用ARMA模型，相对地说,ARMA模型适用范围较宽。对于滚降太快的谱，没有一种模型可以准确地表示功率谱，可以选用高阶的AR模型近似表示。如果选择不合适，例如，选用MA模型去估计具有尖峰的功率谱，估计效果会很差。4.3.1模型选择下图表示的是用MA模型估计二阶AR信号功率谱的例子，图中a)、b)、c)中的AR信号谱峰较平坦，用二阶MA模型功率谱拟和真实谱时，差别较大，随着阶数的提高，估计的谱愈来愈近似于真实的谱。但是对于d)、e)、f)中AR信号的谱峰很窄，在用MA信号模型拟和时，直到MA模型阶数提高到10阶，其效果仍很差对AR(2)信号的模型选择对AR(2)信号的模型选择实际信号中一般都含有和信号不相关的噪声，对带有噪声的信号，如果信号是AR模型，由于噪声的存在需要用ARMA模型；如果用AR模型，则需要阶数更高。选择模型的另一个考虑是尽量减少模型参数。当然这和选择模型是否合适有关系，虽然三种信号模型均有普遍应用价值，当模型选择合适时，估计的功率谱和实际的谱拟合得好，如果不合适，只有提高阶数才能得到较近似的谱，这样需要估计的参数增多，同样也会降低谱估计的质量。因此应该在选择模型合适的基础上，尽量减少模型的参数。go根据观测数据确定模型的参数对各种功率谱估计方法不尽相同。介绍模型参数和信号自相关函数之间的关系，这些关系在功率谱估计中起着很重要的作用假设模型的差分方程和系统函数分别如下式所示：4.3.2模型参数和自相关函数之间的关系在第一章已推导出系统输出功率谱与输入功率谱之间的关系:式中,Pxx(z)和Pww(z)分别表示模型输出信号x(n)和输入信号w(n)功率谱的z变换形式;H(z)表示模型的系统函数。将上式两边同乘以A(z)，得1.ARMA模型的系数和信号自相关函数之间的关系go将上式进行z反变换，公式左边的反变换为公式右边的反变换为式中由于模型H(z)是因果的，h(n)=0，n＜0,可以得到m=0,1,2,3,…,q

m≥q+1令ARMA模型输出自相关函数与模型参数之间的关系式ARMA模型的系数和自相关函数之间的关系如果能由信号的观测数据估计出信号的自相关函数，可以按照上式求出ARMA模型的参数，再由模型参数求出信号的功率谱。但是在公式中由于项的存在，模型输出自相关函数和模型参数之间的关系是非线性的，从而增加了估计功率谱的困难。但是当m>q时，上式是一个线性方程go上式共有p个方程。可以用该方程首先计算出AR部分的p个系数hA(i),

i=1,2,3,…,p;然后代入关系式，设法求出MA部分的系数。将上式以矩阵形式表示：AR模型的系统函数为：H(z)=1/A(z)，相当于ARMA模型中B(z)=1的情况，这样在公式中m≥1m=0

2.AR模型的系数和信号自相关函数之间的关系m=0,1,2,3,…,q

m≥q+1h(0)=?整个关系式用矩阵表示为可以用模型参数表示为：为自相关矩阵，它满足 ，且沿任一对角线的元素相等，它也是正定矩阵上面推导出的关系式或模型系数关系式确定了AR模型参数(包括模型输入噪声方差)和信号自相关函数之间的关系。这是一个线性方程，如果能够由信号的观测数据求出信号的自相关函数，可以按照该公式，通过解线性方程得到模型参数.MA模型的系统函数H(z)=B(z)，相当于ARMA模型中A(z)=1,hA(n)=δ(n)的情况，此时h(n)=hB(n)，由ARMA模型可得到MA模型的系数和自相关函数的关系m=0,1,…,q

m≥q+1MA模型的参数和自相关函数之间也是非线性关系3.MA模型的系数和信号自相关函数之间的关系m=0,1,2,3,…,q

m≥q+1三种信号模型的参数和信号自相关函数之间的关系。这些关系式提供了一种估计功率谱的方法，即首先根据信号观测数据估计信号自相关函数，然后按照所选择的信号模型，解上面相应的方程，求出模型参数，最后按照下式求出信号的功率谱:参数模型信号的功率谱传统的功率谱估计方法无论采取哪一种改进方法，总是以减少分辨率为代价来换取估计方差的减少，其最大缺点是：频谱的分辨率低MA,ARMA模型：自相关函数和模型参数之间呈非线性的关系，增加了估计功率谱的困难用n时刻之前的p个数据为：x(n-1)，x(n-2),…,x(n-p)，预测n时刻的数据x(n)：预测误差为：4.4.1AR模型的线性预测4.4

AR模型谱估计的性质其中，为预测值，为预测系数e(n)表示线性一步预测误差api(i=1,2,3,…，p)表示预测器的系数，它和线性预测器单位脉冲响应h(n)差一负号对差分方程进行z变换输入是观测数据，输出的是预测误差令He(z)=E(z)/X(z)，得到预测误差滤波器的系统函数称He(z)为一步线性预测误差滤波器，其作用是将信号x(n)转换成预测误差e(n)。一般认为e(n)具有白噪声的性质，因此He(z)也称为白化滤波器。预测误差滤波器/白化滤波器AR模型的系统函数为:当api=ai(i=1,2,3,…,p)时，He(z)和H(z)互为逆滤波器，

He(z)=1/H(z)=A(z)白化滤波器的系统函数:预测误差滤波器AR模型当AR模型的阶数与线性预测误差滤波器的阶数相同时，二者互为逆滤波器。AR模型和预测误差滤波器的级联相当于一个全通网络。w(n)是AR模型输入白噪声，故预测误差e(n)具有白噪声的性质，所以，预测误差滤波器也叫白化滤波器。并且，最小预测误差功率就是激励源白噪声的方差σw2。由于AR模型具有这种特性，因而AR模型法也称为线性预测AR模型法。AR模型必须是因果稳定的，即极点均在单位圆内，才能保证信号x(n)是平稳随机信号，由于AR模型H(z)和预测误差滤波器He(z)互为逆滤波器，所以He(z)应为最小相位系统。4.4.2预测误差滤波器的最小相位特性自相关函数和AR模型的模型参数之间的关系服从尤勒-沃克(Yule-Walker)方程：m≥1m=0上式中，对于m≥1的情况，公式本身就是一个递推方程，如果已由观测数据计算出p+1个自相关函数，用,m

=

0,1,2,,…,p表示，对于m＞p的情况,可以用该公式外推得到：4.4.3AR模型隐含自相关函数延拓特性0≤m≤p

m

>p

式中，系数hA(l)是用前p+1个自相关函数求出的参数因此AR模型隐含着自相关函数外推的特性。在经典谱估计BT法中，自相关函数只能限于由观测数据计算出的有限个自相关函数，其它的认为是0，造成了谱估计分辨率低、模糊。正是AR模型具有自相关函数外推的特性，使AR估计的谱具有高分辨率的优点。信号频域的分辨率与时域长度的关系：4.5

AR谱估计的方法自相关法的出发点是选择AR模型的参数使预测误差功率最小，预测误差功率为假设信号x(n)的数据区在0≤n≤N-1范围，有p个预测系数，N个数据经过冲激响应为api(i=0,1,2,…,p)的滤波器，输出预测误差e(n)的长度为N+P4.5.1自相关法——列文森(Levenson)递推法由信号的观测数据，估计信号的功率谱显然,e(n)的长度长于数据的长度，上式中数据x(n)的两端需补充零点，这相当于无穷长的信号经过加窗处理，得到长度为N的数据。用上式对系数api求微分的方法使预测误差功率最小，得到上面的矩阵就是Yule-Walker方程式中的自相关函数采用有偏自相关估计m=0,1,2,…,pm=-p,-p+2,…,-1因此自相关法是基于解Yule-Walker方程的一种方法。首先由信号的观测数据估计出自相关函数，再解该方程，得到模型参数，便可求出信号的功率谱。因此该方法也称为Yule-Walker法。但是直接解该方程，需要计算逆矩阵，很不方便。利用Yule-Walker方程中自相关矩阵的性质，可以导出Levenson-Durbin递推法，这是一种高效的解方程方法。简称为列文森递推法。i=1,2,3,…,k-1由k=1开始递推，递推到k=p，依次得到{a11,σ21}，{a21,a22,σ22}，…，{ap1,ap2,…,app,σ2p}。AR模型的系数以及模型输入白噪声方差求出后，信号功率谱为：上式表明： ，说明随着阶数增加，预测误差功率将减少或者不变，为此要求|akk|≤1,akk称为反射系数。递推公式提供了一种确定模型阶数的方法，如模型的阶数未知，由低阶开始递推，当递推到M阶时，预测误差满足要求，则停止递推，选AR模型的阶数为M。递推法效率高，当阶数变化时，无需重新计算如果已知N个观测数据（x(n),0≤n≤N-1），利用列文森递推法计算功率谱的计算流程图如下图所示利用列文森递推法计算功率谱的流程图k=1,2,3,…,p-1该方法和自相关法一样，仍利用预测误差功率最小方法求模型参数，但求预测误差功率的公式不同：对比自相关法求预测误差功率的公式，不同的是求和限不同。该公式使用的观测数据是已知的，不需要在数据两端补充零点，因此与自相关法相比，去掉了加窗处理的不合理假设。4.5.2协方差法与修正协方差法1.协方差法白噪声的方差为仍然使用梯度最小的方法求模型参数：协方差函数观测数据x(n)(n=0,1,2,…,N-1)，利用上面公式可以求出模型的参数：{api(i=1,2,3,…,p);σ2w}。式中的协方差函数cxx(j,k)，有两个变量，因此也适合于非平稳随机信号。式中的协方差矩阵是埃尔米特（Hermitian）矩阵，是半正定的。这种方法近似于自相关法。一些实验结果说明它的分辨率优于自相关法，另外对于纯正弦信号数据，可以有效地估计正弦信号的频率。修正协方差法使用前向和后向预测误差平均值最小的方法，估计AR模型的参数，从而得到信号的功率谱。信号的前向和后向预测分别如下：2.修正协方差法式中apk是AR模型的参数最小预测误差平均功率是模型输入白噪声的方差，即ρp=σ2w，前、后向预测误差平均功率为前向和后向预测误差功率ρpe、ρpb分别用下式表示和协方差法一样，上式仅对用到的观测数据的预测误差求和。为了使预测误差平均功率最小，求ρp对apk(k=1,2,3,…,p)的微分,或者用复梯度法求，得到化简并写成矩阵形式为：白噪声的方差估计值为例:已知信号的四个观察数据为x(n)={x(0),x(1),x(2),x(3)}={2,4,1,3}，分别用自相关法和协方差法估计AR(1)模型参数。1)自相关法：解2)协方差法:x(n)={2,4,1,3}4.5.3伯格(Burg)递推法设信号x(n)观测数据区间为：0≤n≤N-1，前向、后向预测误差功率分别用ρp,e和ρp,b表示，预测误差平均功率用ρp表示，公式分别为前向、后向预测误差递推公式如下：将上式带入平均误差公式中，得到求预测误差平均功率ρp最小时的反射系数kp，上式就是利用伯格递推法求第p个反射系数的公式将伯格递推法求AR模型参数的递推公式总结如下：利用伯格递推法求AR模型参数的流程图如下图所示伯格递推法流程图书159页n=p…N-1p=0

计算AR模型参数应用较多的为以下三种方法：1、自相关法---（Levinson法）2、协方差法/改进协方差法3、Burg法自相关法（Levinson法）计算简单，但由于假设数据外为零使得分辨率相对较差协方差法与自相关法的主要区别是：求和范围不同。数据段两端不需要添加任何零取样值，没有假设数据段N以外的数据等于零。谱估计性能较好，但潜在不稳定因素改进协方差法谱估计性能最好，但计算过于复杂Burg法不直接估计AR参数，而是先估计反射系数，然后利用Levinson算法由反射系数计算AR参数计算反射系数时，使用前向、后向预测误差功率的平均值最小准则。计算不太复杂，且给出了较好的谱估计质量.是较为通用的方法在AR模型谱估计中，模型阶数的选择是一个关键问题。一般模型的最好选择是先验未知的，实际中需预先选定模型阶次。如果是纯p阶AR信号，选择模型阶次k＜p时，将产生对谱的平滑作用，降低谱的分辨率，如下图所示。图中，AR信号p=4，选择模型的阶次k=2，产生的平滑作用使两个峰变成一个峰，分辨率明显降低。如果选择k≥p，且假定观测的数据没有误差（没有干扰），估计的参数应是：4.5.4关于AR模型阶次的选择AR模型阶次太小时的平滑作用因此,对于纯p阶AR信号，应选择阶次k≥P。如果是白噪声中的AR信号(观测数据有误差或者信号中含有白噪声)，此时选择ARMA模型合适，如选择了AR模型，其阶次应加大，较低的阶次会使谱估计产生偏移，降低分辨率。当然,这也和信噪比有关，信噪比愈低，平滑作用愈严重，愈需要高的阶次，因此信噪比低应选高的阶次。一般来说，阶次愈高，分辨率愈高；但阶次太高，会因参数过多而使估计误差加大，谱峰分裂，因此，对于白噪中的AR信号，其阶次的选择应折衷考虑4.6最大熵谱估计与最大似然谱估计4.6.1最大熵谱估计按照Shannon对熵的定义，当随机变量取离散值时:式中pi是i的概率。当X取连续值时，熵的定义为1.利用最大熵的原则外推自相关函数式中,p(x)是X的概率密度函数，对于离散随机序列,概率密度函数用联合概率密度函数代替。显然,熵代表一种不确定性，最大熵代表最大的不确定性，或者说最大的随机性。下面我们研究对于有限的自相关函数值不作任何改变，对于未知自相关函数用最大熵原则外推，即不作任何附加条件的外推方法。假设x(n)是零均值正态分布的平稳随机序列，它的N维高斯概率密度函数为式中

按照（4.6.2）式，x(n)信号的熵为

（4.6.3）

式中det(Rxx(N))表示矩阵Rxx(N)的行列式，由上式表明为使熵最大，要求det(Rxx(N)最大。若已知N+1个自相关函数值rxx(0),rxx(1),…,rxx(N)，下面用最大熵方法外推rxx(N+1)。设rxx(N+1)确实是信号自相关函数的第N+2个值，根据自相关函数的性质，由N+2个自相关函数组成的矩阵为（4.6.4）

它必须是非负定的矩阵，

即

（4.6.5）

将行列式展开，det(Rxx(N+1))是rxx(N+1)的二次函数，该二次函数系数的符号是：(-1)1+N+2(-1)1+N+1=-1，且det(Rxx(N+1))对rxx(N+1)的二次导数是-2det［Rxx(N-1)］，它是负值，负值表示det(Rxx(N+1))对rxx(N+1)的一次导数是减函数，det(Rxx·(N+1))作为rxx(N+1)的函数，凹口向下，那么只有一个最大值。为选择rxx(N+1)使det(Rxx(N+1)最大，解下列方程：(4.6.6)用数学归纳法，得到

（4.6.7）

上式是rxx(N+1)的一次函数，可以解出rxx(N+1)。继续再将rxx(N+1)代入Rxx(N+2)和det(Rxx(N+2))中，求det(Rxx(N+2))对rxx(N+2)的最大值，得到rxx(N+2);以此类推，可推出任意多个其它自相关函数值，而不必假设它们为零，这就是最大熵谱估计的基本思想。2.最大熵谱估计与AR模型谱估计的等价性我们已经知道AR模型信号自相关函数与模型参数服从Yule-Walker方程，即

m≥1m=0将m≥1的情况写成矩阵形式：

式中,ai=hA(i),i=1,2,3,…,N,ai是AR模型系数。解该方程，可以得到模型系数ai,即（4.6.8）

（4.6.9）

（4.6.10）

(4.6.11)在（4.3.6）式中，令m=N+1，得到

(4.6.12)将以上求出的系数a1,a2,…,aN代入上式，求出rxx(N+1)。而最大熵外推自相关函数的公式是(4.6.7)式，按照该公式的最后一行展开，得到（4.6.13）

上式即是最大熵外推自相关函数的公式，对比（4.6.12）式，两公式完全一样，证明了AR模型功率谱估计和最大熵谱估计的等价性。这里最大熵外推自相关函数等价于已知N+1个自相关函数，匹配一个N阶AR信号模型的系数。一旦通过解Yule-Walker方程，解出模型参数，最大熵谱估计用下式计算信号功率谱：

（4.6.14）

最大似然谱估计是用一个FIR滤波器实现，该滤波器对所关心频率的正弦信号，可以无失真地通过，而对于其它频率的信号，让其频响尽可能地小，亦即将它们尽可能地滤除。此时,滤波器输出的均方值，就作为信号的功率谱估计。设实信号用x(n)表示，FIR滤波器系统函数用A(z)表示：

输出y(n)为

(4.6.15)4.6.2最大似然谱估计——最小方差谱估计式中

输出信号的均方值为

(4.6.16)上式中T表示转置,H表示共轭转置,Rp=E［XXT］是Toeplith自相关矩阵,为求,必须先求FIR滤波器的系数。求这些系数的原则是：在所关心频率ωi处，信号x(n)无失真地通过，即在ωi处的传输函数为1：

式中

(4.6.17)另外一个原则是在ωi附近的频率分量尽量衰减掉，即ω≠ωi处,滤波器输出y(n)的均方差最小,即(4.6.16)式最小,此时作为信号x(n)的功率谱估计 。因此,最大似然谱估计称为最小方差谱估计更为合适，但由于习惯也可以仍称为最大似然谱估计。在以上原则下，使方差最小的滤波器系数和分别为［30］、［31］

应该指出,此时并不是真正意义上的信号功率谱,只是描述了信号功率谱的相对强度。

下面分析最小方差谱估计与AR模型谱估计之间的关系：

（4.6.20）

(4.6.21)4.7特征分解法谱估计无论是实正弦波还是复正弦波，都可以用一个退化AR模型表示，设P个实正弦波组成的信号用下式表示：

（4.7.1）

式中，初相位θi是在区间（-π，π）均匀分布的随机变量，首先分析下面的三角恒等式：

-π＜ω＜π

4.7.1正弦波用退化AR模型表示令x(n)=sin(ωn+θ),则上式变为

（4.7.2）

将上式进行Z变换，得到

（4.7.3）

这样(4.7.2)式的特征多项式为

（4.7.4）

上式的两个根分别是:z1=ejω,z2=e-jω,它们共轭成对，且模为1。由这两个根可以确定正弦波的频率。对比AR模型的系统函数，可以把正弦波信号用一个特殊的AR（2）模型表示，括弧中的2表示模型是二阶的。该AR模型的激励白噪声方差趋于0，极点趋于单位圆。通常称为退化的AR模型。这一模型系数有两个，即2cosω和1，（4.7.2）式是模型的差分方程。

对于P个实正弦波,特征多项式是

（4.7.5）

上式是z-1的2P阶多项式，可以表示为

（4.7.6）

注意上式中的系数ak(k=1,2,3,…,2P),必须保证它的根共轭成对。考虑到根共轭成对，也可表示为

（4.7.7）

这样由(4.7.6)式，P个正弦波组合的模型用下面2P阶差分方程描述

（4.7.8）

对于复正弦波情况，P个复正弦波组成的信号是

（4.7.9）

用一个退化的AR（p）模型表示的差分方程为

（4.7.10）

其特征多项式为

(4.7.11)其根为

1≤i≤P

注意这里的根不是共轭成对出现的。总结以上P个正弦波组合是一个退化的AR(2P)过程，独立参量个数为P个；P个复正弦波的组合是退化的AR（P）过程，独立参量个数仍为P个。实正弦过程相应的退化AR过程的阶数比复正弦情况的阶数高1倍。

白噪声中正弦波组合的信号为

（4.7.12）

式中，w(n)为白噪声,且

将(4.7.12)式中x(n)的用AR(2P)表示，即将(4.7.8)式带入(4.7.12)式中，

得到

（4.7.13）

4.7.2白噪声中正弦波组合用一特殊的ARMA模型表示将（4.7.12）式中的n用n-i代替,x(n-i)=y(n-i)-w(n-i)再将上式带入(4.7.13)式,得到

(4.7.14)上式可以看成一个特殊的ARMA(2P,2P)模型，括弧中的两个2P分别表示ARMA模型系统函数分子和分母的阶次。它与一般的ARMA模型比较，有三方面不同：

(1)它的AR部分和MA部分具有相同的参数，它们存在共同的因子；(2)由于特征多项式(4.7.6)式的根的模为1，故AR部分特征多项式不满足平稳性条件，MA部分特征多项式也不满足可逆性条件；（3）AR部分的y(n)=x(n)+w(n)，y(n)是含白噪声的观测值，而通常为信号的x(n)不含白噪声。这种特殊的ARMA模型结构，不能用一般的ARMA模型结构求解。下面介绍特征分解技术。将（4.7.14）式写成矩阵形式：

YTA=WTA

(4.7.15)式中

4.7.3特征分解法谱估计用向量Y左乘(4.7.15)式,并取数学期望,得到

E［YYT］A=E［YW

T］A

(4.7.16)式中

将上面关系式带入(4