1、第4章 功率谱估计第4章 功率谱估计4.1 经典功率谱估计4.2 AR模型功率谱估计的方法和性质4.3 最大熵谱估计方法4.4 最大似然谱估计4.5 互协方差估计与互谱估计4.6 特征分解法谱估计第4章 功率谱估计 4.1 经典功率谱估计经典功率谱估计4.1.1 BT法法根据维纳-辛钦定理，1958年Blackman和Tukey给出了这一方法的具体实现，即先由N个观察值xN(n)估计出自相关函数rx(m)，求其傅立叶变换，以此变换结果作为对功率谱Px()的估计。如果我们得到的是x(n)的N个观察值x(0)，x(1)，x(N1)，令 xN(n)=a(n)x(n)(4.1.1)其中a(n)是数据窗
2、，对于矩形窗 第4章 功率谱估计 1,01()0,0nNa n计算rx(m)的估计值的一种方法是10101()()()1()(),|1NxNNnNmNNnrnxn xnmNxn xnmmNN(4.1.2)的均值为()xrm110011()()()()()()NmNmNNNNmmExn xnmE xn xnm a n a nmNN 若a(n)是矩形窗，则()()xxNmE r mr mN第4章 功率谱估计所以，偏差(4.1.4)bias()()()()xxxxmr mE r mr mr mN 由此可以看出:(1)这种自相关函数的估计是一个有偏估计，且估计的偏差是，当N+时,。因此是rx(m)的渐
3、近无偏估计。(2)对于一个固定的N，当|m|越接近于N时，估计的偏差越大。(3)由式(4.1.3)可看出，是真值rx(m)和三角窗函数()xmr mNbias()0 xr m()xr m()xE r m第4章 功率谱估计 101()0mmNq mN其他,的乘积，q(m)的长度是2N1，它是由矩形数据窗ar(n)的自相关所产生的。的方差是()xr m(4.1.5)22var()()()xxxr mE rmE r m而(4.1.6)1122001()()()()()NmNmxmkE rmEx n x nmx k x kmN 21()()()()nkE x n x k x n m x k mN第4章
4、 功率谱估计假定x(n)是零均值的高斯随机信号，有Ex(n)x(k)x(n+m)=r2x(nk)+rx(nkm)rx(nk+m)所以22221()()()()()xxxxxnkE rmrnkrmr nkm r nkmN2221()()()()xxxxnkNmr mrnkr nkm r nkmNN将上式和式(4.1.3)代入式(4.1.5)得1122001var()()()()Nm Nmxxxxnkrmrnkr nkm r nkmN (4.1.7)第4章 功率谱估计令nk=i，式(4.1.7)可写成1211var()1()()()NmxxxxiNmmirmrir im r imNN (4.1.8
5、)在大多数情况下，rx(m)是平方可求和的，所以当N+时，varrx(m)0，又因为lim biasrx(m)=0，所以对于固定的延迟|m|，rx(m)是rx(m)的一致估计。对由式(4.1.2)得到的自相关函数估计rx(m)进行傅立叶变换:(4.1.9)jBT()()()e,|1MmxmMPv m rmMN第4章 功率谱估计其中，v(m)是平滑窗，其宽度为2M+1，以此作为功率谱估计，即为BT谱估计。因为用这种方法求出的功率谱是通过自相关函数的估计间接得到的，所以此法也称为间接法。第4章 功率谱估计4.1.2 周期图法周期图法周期图法是把随机信号的N个观察值xN(n)直接进行傅立叶变换，得到
6、XN(ej)，然后取其幅值的平方，再除以N，作为对x(n)真实功率谱Px()的估计。以Pper()表示周期图法估计的功率谱，则(4.1.10)2per1()NPXNj（e）其中je)=11jj00()e()()eNNnnNNNnXxnx n a n(4.1.11)a(n)为所加的数据窗，若a(n)为矩形窗，则j(e)=1j0()eNnNNXx n第4章 功率谱估计因为这种功率谱估计的方法是直接通过观察数据的傅立叶变换求得的，所以习惯上又称之为直接法。周期图法功率谱估计的均值为d)()21)(per(PPEx(4.1.12)(4.1.13)222sin2sin1|)(|1)(NNAN其中第4章
7、功率谱估计4.1.3 周期图法与周期图法与BT法的关系法的关系式(4.1.9)中取M为其最大值N1，且平滑窗v(m)为矩形窗，则 1jBT(1)()()eNnxmNPr m11j(1)01()()()()eNNnmNna n a nm x n x nmN 11jj()0(1)1()()e()()eNNnn mnmNa n x na nm x nmN令n+m=l，上式可变成11jjBT001()()()e()()eNNnlnlPa n x na l x lN2j1(e)NXN第4章 功率谱估计所以BTper1()()MNPP(4.1.14)由此可见，周期图法功率谱估计是BT法功率谱估计的一个特例
8、，当间接法中使用的自相关函数延迟M=N1时，二者是相同的。第4章 功率谱估计1.Welch法法假定观察数据是x(n)，n=0，1，N1，现将其分段，每段长度为M，段与段之间的重叠为Mk。如图 4.1所示，第i个数据段经加窗后可表示为 xiM(n)=a(n)x(n+ik)，i=0,1,L1,n=0,1,M1其中，k为一整数，L为分段数，它们之间满足如下关系:(L1)k+MN该数据段的周期图为(4.1.15)per1()()iiMPXMU其中1j0()()eMiinMMnXXn(4.1.16)第4章 功率谱估计U为归一化因子，使用它是为了保证所得到的谱是真正功率谱的渐进无偏估计。由此得到平均周期图
9、图 4.1 数据分段方法(4.1.17)1perper01()()LiiPPL第4章 功率谱估计如果x(n)是一个平稳随机过程，每个独立的周期图的期望值是相等的，根据式(4.1.15)和式(4.1.16)有(4.1.18)perper1()()()()d2ixE PE PPQ其中(4.1.19)21()QAMUA()是对应M个点数据窗a(n)的傅立叶变换，若M值较大，则Q()主瓣宽度较窄，如果Px()是一慢变的谱，那么认为Px()在Q()的主瓣内为常数，这样式(4.1.17)可以写成(4.1.20)per1()().()d2xE PPQ第4章 功率谱估计为了保证Welch法估计的谱是渐进无偏的
10、，必须保证(4.1.21)1()d12Q或(4.1.22)211.()d12AMU根据Parseval定理，式(4.1.22)可写成1201.()1MnanMU所以归一化因子U应取成(4.1.23)1201.()MnUanM第4章 功率谱估计的方差表达式为 per()P(4.1.24)如果x(n)是一个平稳随机过程，上式的协方差仅仅取决于il=r，令(4.1.25)11perperper2001var()cov(),()ilLLilPPPLperper()cov(),()ilrPP式(4.1.24)可写成单求和表示式(4.1.26)1perper(1)0()1var()var()1()iLrr
11、LrPPLL 其中表示某一数据段的周期图方差，即pervar()Pperper0var()var()(),0,1,1iPPiL 第4章 功率谱估计而是与的归一化协方差，如果各个数据段的周期图之间的相关性很小，那么式(4.1.26)可近似写成0()()rper()iPper()i rP(4.1.27)perper1var()var()iPPL这也就是说,平均周期图的方差减小为单数据段图方差的1/L。第4章 功率谱估计2.Bartlett对应Welch法，如果段与段之间互不重叠，且数据窗选用的是矩形窗，此时得到的周期图求平均的方法即为Bartlett法。可以从上面讨论的Welch法得到Bartle
12、tt法有关计算公式，第i个数据段可表示为 xiM(M)=x(n+iM)，i=0，1,L1,n=0,1,M1 (4.1.28)其中，LMN。该数据段的周期图为(4.1.29)2per1()()iiMPXM其中(4.1.30)1j0()()eMiinMMnXxn第4章 功率谱估计平均周期图为(4.1.31)1perper01()()iLiPPL其数学期望为(4.1.32)perper1()()()()d2ixE PE PPQ其中2sin12()sin2MQM(4.1.33)第4章 功率谱估计将式(4.1.33)与式(4.1.13)相比，取平均情况下A()的主瓣宽度是不取平均情况下A()的主瓣宽度的
13、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章 功率谱估计4.1.5 经典功率谱估计性能比较经典功率谱估计性能比较为了对几种经典功率谱估计
14、方法的性能进行比较，我们采用参考文献27中的数据。信号表示为 4j21()()ekf nkkx ny nA其中包含有四个复正弦，其归一化频率分别是f1=0.15，f2=0.16，f3=0.252，f4=0.16,Ak对应不同的系数，可得到不同的信噪比，本数据在f1处的信噪比为64 dB，在f2处的信噪比为54 dB，在f3处的信噪比为2 dB，在f4处的信噪比为30 dB。y(n)是一个复值的噪声序列，其功率谱为242j1()21()ekykPb k其中，2=0.01，b(k)是模型系数。第4章 功率谱估计x(n)的真实功率谱曲线如图4.2(a)所示，注意其频率范围是从0.50.5，即。令归一
15、化频率f1和f2相差0.01，目的是检验算法的分辨能力；f3的信噪比很低，目的是检验算法对弱信号的检测能力。现取N128，图4.2(b)示出了该数据段直接求出的周期图，所用数据窗为矩形窗，由于主瓣过零点宽度B=2/128=0.015 6250.01，所以f1和f2不能完全分开，只是在波形的顶部能看出两个频率分量。图4.2(c)是利用Welch平均法求出的周期图，共分四段，每段32点，没有重叠，使用Hamming窗，这时谱变得较平滑，但分辨率降低。第4章 功率谱估计图 4.2 功率谱估计方法比较第4章 功率谱估计图4.2(d)也是用Welch平均法求出的周期图，共分七段，每段32点，重叠16点，
16、使用Hamming窗，谱变得更加平滑，分辨能力和图4.2(c)大体一致。图4.2(e)是用BT法求出的功率谱曲线，M=32，没有加窗；图4.2(f)也是用BT法求出的欧尼功率谱曲线，M=16，使用了Hamming窗。第4章 功率谱估计4.2 AR模型功率谱估计的方法和性质模型功率谱估计的方法和性质4.2.1 AR模型功率谱估计的引出模型功率谱估计的引出由第2章讨论可知，p阶AR模型满足如下差分方程:xAn+a1xAn1+apxAnp=n其中，a1，a2，ap为实常数，且ap0；n是均值为0、方差为2的白噪声序列，也就是说，随机信号xAn可以看成是白噪声n通过一个系统的输出，如图4.3所示。第4
17、章 功率谱估计图 4.3 AR模型信号第4章 功率谱估计图4.3中 1()()H zA z(4.2.1)而A(z)=1+a1z1+apzp(4.2.2)已经证明(4.2.3)121(),0()(),0pk AkApk Aka r mkmr ma r km其中，rA(m)是AR模型的自相关函数，尤其对于0mp，由式(4.2.3)可写出矩阵方程为第4章 功率谱估计 (4.2.4)211(0)(1)()(1)(0)(1)0()(1)(0)0AAAAAApAAArrrparrrparprpr这就是AR模型的正则方程，又称Yule-Walker方程。对于一个p1阶预测器，预测值为(4.2.5)11()(
18、)()()()ppkkx na k x nkh k x nk其中，h(k)=a(k)，预测误差为(4.2.6)0()()()()()pke nx nx na k x nk其中，a(0)=1。p阶预测误差滤波器Ap(z)如图4.4所示。第4章 功率谱估计图 4.4 p阶预测误差滤波器第4章 功率谱估计图中，Ap(z)=1+a(1)z1+a(p)zp当Ee2(n)达到其最小值Ee2(n)min时，必满足Yule-Walker方程(4.2.7)2min(0)(1)()1(1)(0)(1)(1)0()(1)(0)()0 xxxxxxxxxE enrrr prrr par pr pra p当x(n)就是
19、图4.3所产生的p阶AR过程xAn，也即xn=xAn或rx(m)=rA(m)时，m=0，1，p，必满足关系式 22min()1,2,kaa kkpE en(4.2.8)第4章 功率谱估计或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模型，那么以方差为Ee2(n)min的白噪声n通过此线性系统，其输出功率谱必定
20、与待估计的随机信号的功率谱完全相同，因此，模型H(z)可以完全表示出AR(p)过程的xn功率谱，它们的关系即是第4章 功率谱估计(4.2.10)2AR2j1()()1expkkkPPa当然，实际待估计的随机信号xn可能是一个阶数大于p的AR过程，也可能根本就不是一个AR过程，但仍可以采用上述方法建立一个p阶的AR模型，作为对随机信号xn的功率谱估计，此功率谱估计可作为xn真实功率谱的一个近似，其步骤是:(1)对此随机信号xn建立p阶的线性预测误差滤波器，求得系数a(1)，a(2)，a(p)和Ee2(n)min；(2)令A(z)=1+a1z1+apzp，其中a1=a(1)，ap=a(p)，并构成
21、一线性系统(4.2.11)1()()H zA z第4章 功率谱估计那么将一方差为2的白噪声n通过该系统，其输出的功率谱可作为待估计随机信号xn的功率谱估计(4.2.12)22j1()1expkkkPa其中，2=Ee2(n)min。有下述几点需要注意:(1)预测误差滤波器Ap(z)实际上就是一个白化滤波器，但一般情况下，它不能将随机信号xn进行完全白化，所以图4.4中预测误差滤波器的输出en并不是白噪声n,因此Ee2(n)min表示的也是非白噪声的方差；只有当被估计的随机信号xn本身就是一个AR(p)过程,且当预测误差滤波器的阶数大于或等于p时，其输出en才是一个白噪声n，Ee2(n)min 即
22、表示此白噪声的方差。第4章 功率谱估计(2)对比式(4.2.4)和式(4.2.7)，当ak=a(k)，k=1，2，p和Ee2(n)min时，用这样的方法所建立的p阶AR模型，其自相关函数rA(m)和待估计随机信号xn的自相关函数rx(m)必有如下关系:(4.2.13)1()0()()xpAk Akr mmpr ma r mkmp由式(4.2.13)可见，用AR(p)模型来对任意的随机信号建模，一定具有自相关函数的匹配性质，随着模型阶数p的增大，匹配的程度会越来越好；当p+时，AR模型的自相关函数和随机信号xn的自相关函数可达到完全的匹配，即(4.2.14)lim()()0Axpr mr mm
23、第4章 功率谱估计而当p+时，p阶FIR预测误差滤波器退化为一个IIR预测误差滤波器，它一定可将随机信号xn进行完全的白化。为了说明这个问题，我们又重画出预测误差滤波器(如图4.5所示)，图中G(z)是随机信号xn的建立模型。图 4.5 AR模型与预测误差滤波器第4章 功率谱估计当p+时1()1()HzzG z因此，预测误差滤波器(4.2.15)1211()1(1)(2)111()1(1)()()Azazazz Hzz zG zG z 这证明了上述预测误差滤波器可将随机信号完全白化的结论。第4章 功率谱估计4.2.2 AR模型谱估计的性质模型谱估计的性质 1.隐含的自相关函数延拓的特性隐含的自
24、相关函数延拓的特性在前面讨论的经典BT法功率谱估计中，假定由给定的数据xN(n)，n=0，1，N1，可估计出自相关函数rx(m)，m=(N1)(N1)，在这个区间以外，用补零的方法将其外推，对此求其傅立叶变换(4.2.16)1jBT(1)()()eNmxmNPrm就可得到BT法的功率谱估计PBT()，此PBT()的分辨率显然是随着信号长度N的增加而提高的。第4章 功率谱估计而在AR模型谱估计中，上述限制不再存在。虽然给定的数据xN(n)，n=0，1，N1,是有限长度，但现代谱估计的一些方法，包括AR模型谱估计法，隐含着数据和自相关函数的外推，使其可能的长度超过给定的长度。前面讨论的AR模型的建
25、立，用到了单步预测的概念，预测值为1()()pkkx na x nk 这样x(n)可能达到的长度是1(N1+p)，如果在递推的过程中，用x(n)代替x(n)，那么还可继续不断地外推。同样，从AR模型的建立过程看，AR过程的自相关函数必满足1(),0()(),xpAk Akr mmpr ma r mkmp第4章 功率谱估计由此式可见，AR模型的自相关函数在0p范围内与rx(m)完全匹配，而在这区间外，可用递推的方法求得。自相关函数rA(m)实际上就是被估计信号xn的自相关函数的估计jj()()e()emmxxAmmPrmrm将其进行傅立叶变换，就可得到随机信号xn的功率谱(PSD)估计，即()(
26、)xArmr m(4.2.17)比较式(4.2.16)和(4.2.17)可见，AR模型法避免了窗函数的影响，因此它可得到高的谱分辨率，同时它所得出的功率谱估计Px()与真实的功率谱Px()偏差较小。图4.6示出了AR模型谱估计和BT谱估计法的比较。第4章 功率谱估计图 4.6 AR模型谱估计和BT谱估计性能比较第4章 功率谱估计2.AR模型的稳定性模型的稳定性AR模型稳定的充分必要条件是H(z)=1/A(z)的极点(即A(z)的零点)必须在单位圆内，而且这一条件也是保证x(n)是一个广义平稳过程所必需的。这很容易证明，如果H(z)有一个极点在单位圆外，那么x(n)的方差将趋于无穷，则x(n)是
27、非平稳的。因为AR模型的系数是由Yule-Walker方程解得的，可以证明如果(p+1)(p1)的全相关矩阵(0)(1)()(1)(0)(1)()(1)(0)xxxxxxpxxxrrr prrr pRr pr pr第4章 功率谱估计是正定的，AR模型的系数具有非零解，此时预测误差滤波器A(z)一定具有最小相位的特性；换句话说，AR模型H(z)一定是一个稳定的全极点滤波器。使用AR模型对纯正弦信号建模是不合适的，因为此时Rp可能出现奇异或非正定的情况，但在信号处理中经常要用正弦信号作为试验信号以检验某个算法或系统的性能。为克服自相关矩阵奇异的情况，最常用的方法是加上白噪声n，这样det(Rp)就
28、不会等于零。第4章 功率谱估计3.谱的平坦度谱的平坦度前面的讨论已经指出，AR(p)的系数ak就是预测误差功率最小时的p阶线性预测误差滤波器的系数，由于预测误差滤波器是一个白化滤波器，它的作用是去掉随机信号x(n)的相关性，在自己的输出端得到白噪声n，因此在这一节中，把白化的概念加以推广，表明AR参数也可以使用预测误差滤波器Ap(z)的输出过程具有最大的谱平坦度的方法得到。利用谱平坦度的概念可以把AR谱估计得到的结果看成是最佳白化处理的结果。功率谱密度的谱平坦度可定义为(4.2.18)1ln()d21()d2exxPxP第4章 功率谱估计它是Px()的几何均值与算术均值之比，可以证明 0 x1
29、如果Px()有很多峰(也就是它的动态范围很大)，例如，在由p个复正弦所组成的随机信号j()1()ekkpknkkx nA的式子中，Ak、k为常量，k是在(，+)范围内均匀分布的随机变量，则j21()ekpmxkkr mA21()()pxkkkP mA 第4章 功率谱估计此种信号的功率谱具有最大的动态范围。将上式代入式(4.2.18)，分子显见为零，因此x=0。但如果Px()是一个常数(也就是它的动态范围为零)，也即相当于xn是一个白噪声，则由式(4.2.18)显见x=1，由此可见，谱平坦度x直接度量了谱的平坦程度。现设预测误差滤波器1()1()pkpkAza k z 为最小相位，输入时间序列x
30、(n)是任意的(不一定是AR过程)，按照使输出误差序列e(n)的谱平坦度最大的准则来确定预测系数，为此，引入下述结果:如果Ap(z)是最小相位，则第4章 功率谱估计(4.2.19)21ln|()|d02pA计算预测误差滤波器输出过程e(n)的平坦度。因为 211lndlnd221lnd2epxxPAPP对上式两端取指数并除以 ，得 1lnd2eP 1lnd21d0e2110dd22ePxxexxeeePrrPP第4章 功率谱估计由于对于随机信号x(n)来说，x和rx(0)均是固定的，因此要使e最大，必须使re(0)最小，因为re(0)=Ee2(n)，因此使预测误差e(n)的功率谱平坦度最大和使
31、p阶预测误差滤波器输出的误差功率最小是等效的，亦即条件max e和条件a 22minminaBE enEx nx n 完全等效。如果x(n)本身就是一个AR(p)过程，也即 22exPA第4章 功率谱估计其中，A()=1+a1ej+apejp。现使其通过一个p阶预测误差滤波器Ap()=1+a(1)ej+a(p)ejp 在满足maxe的条件下，一定有A()=Ap()，也即aak=a(k),k=1，2，p此时预测误差滤波器输出的误差序列e(n)一定是一个白噪声序列。反之，如果将AR(p)通过一个k阶预测误差滤波器kp，同样在满足maxe的条件下，误差序列e(n)不可能是一个白噪声序列，这一结果与前
32、面讨论的AR模型谱估计的引出中所得的结果是完全一致的。a第4章 功率谱估计在这里有一个重要的概念需强调，即预测误差滤波器的输入、输出功率谱总满足关系式 2expPPA在满足Ee2(n)min的条件下，根据求得的Ap()建立AR模型:22AxpPPA第4章 功率谱估计其中，2=Ee2(n)min，是一个常数。比较这两式可见，用PA()作为Px()的一个估计，其估计的好坏完全取决于Pe()与一个常量相逼近的程度，换句话说，在建立AR模型时，正是由于用2代替了Pe()，才使得建立的AR模型功率谱PA()中丧失了很多Pe()的重要细节，而只有当误差序列e(n)是一个白噪声序列，Pe()是一个常量，并且
33、就等于2时，才能得到Px()=PA()=Px()。第4章 功率谱估计4.2.3 AR模型参数提取方法模型参数提取方法由前面所讨论的AR模型建立过程可知，AR模型参数必满足Yule-Walker方程:211011010100 xxxxxxpxxxrrrparrrparprpr 第4章 功率谱估计对于某随机过程来说，如果已知rx(0)、rx(1)、rx(p)，就可由上式求解出(2，a1，ap)，但实际上，对于一个所要估计的随机过程xn来说，通常已知的仅仅是有限的观察数据，如x(0)，x(1)，x(N1)，而其自相关函数值通常是未知的，因此通常是首先求得自相关函数的估计值，然后再通过某种算法，得到A
34、R模型参数的估计值。下面讨论常用的几种AR模型参数的提取方法。第4章 功率谱估计1.自相关法自相关法假定观察到的数据是x(0)，x(1)，x(N1)，而对于无法观察到的区间(即nN1)，x(n)的样本假定为零，预测误差功率的表示式为(4.2.20)2121011()()()1()()()pnkNppnkx na k x nkNx na k x nkN 式(4.2.20)的求和限又可以写为0，N1+p，这是因为在此区间外误差表示式总为0的缘故。为了得到预测误差功率的极小值，将式(4.2.20)对a(l)求微分并令其等于0，得第4章 功率谱估计(4.2.21)1011()()()()0,1,2,N
35、ppnkx na k x nkx nllpN 令(4.2.22)101()()0,1,.,()()1,2,.,Npxnxx n x nkkpNrkrkkp 代入式(4.2.21)得(4.2.23)1()()(),1,2,pxxkrla k rlk lp 写成矩阵形式，这组方程为(4.2.24)011(1)(1)102(2)(2)120()()xxxxxxxxxxxxrrrparrrrparrprpra prp 第4章 功率谱估计求出白噪声方差2的估计值2，它是 2min111()()()()()()ppnklx na k x nkx na k x nlN将式(4.2.21)代入上式，得2min
36、111()()()()(0)()()pnkpxxkx na k x nkx nNra k rk将式(4.2.25)与式(4.2.23)合并，得(4.2.26)(4.2.25)min101(1)10100()10 xxxxxxxxxrrrparrrpa prprpr 第4章 功率谱估计或写成(4.2.27)minp pRa0例例 4.1 利用有限长加窗数据()()0,1,.,1()0n xnnNx n其他其中，a(n)为矩形窗，x(n)=rnu(n)，|r|1，采用自相关法对其用一阶的AR模型建模。解解:首先计算rx(0)和rx(1):2120111(0)1NNnnxnrrrrNNr第4章 功率
37、谱估计由此得222120111(1)1NNnnxnrrrrrNNr 222(1)1(1)(0)1NxNxrrarrr 而a(1)的真正解应该是r，这说明当观察数据为有限值N时，其解总是有偏的，直到N+时，a(1)才收敛于r。第4章 功率谱估计2.协方差法协方差法假定观察到的数据是x(0)，x(1)，x(N1)，写出预测误差功率的表示式为(4.2.28)2111()()()pNn pkx na k x nkNp比较式(4.2.28)和式(4.2.20)可见，协方差法与自相关法的区别主要在于预测误差功率求和式的上下限取的不同。由于协方差法对于观察区间0，N1外的x(n)样本并未假定为零，这就要求式
38、(4.2.28)中x(nk)总是落在观察区间0，N1中，为此预测误差功率的求和上下限必须取在p，N1之间。第4章 功率谱估计将式(4.2.28)对a(l)求微分，井令其等于0，以得到预测误差功率的极小值，得(4.2.29)111()()()(),1,2,pNn pkx na k x nkx nl lpNp令11()()(,),1,2,(,)Nn pxxx nl x nkNprl kl kprk l将其代入式(4.2.29)，得1()(,)(,0),1,2,pxxka k rl krllp 写成矩阵形式，这组方程为第4章 功率谱估计1,11,21,(1)(1,0)2,12,22,(2)(2,0)
39、,1,2,()(,0)xxxxxxxxxxxxrrrparrrrparrprprp pa prp (4.2.30)求出白噪声方程2的估计值2:12min111()()()()()()ppNn pklx na k x nkx na l x nlNp将式(4.2.29)代入上式，得12min11()()()()pNn pkx na k x nkx nNp1(0,0)()(0,)pxxkra k rk第4章 功率谱估计将上式与式(4.2.30)合并，得min10,00,10,(1)1,01,11,00(),0,1,xxxxxxxxxrrrparrrpa prprprp p 或写成(4.2.31)mi
40、np pRa0比较式(4.2.31)和式(4.2.27)，表面上看，协方差法和自相关法最终所得到的正则方程具有相同的形式，但实际上两者具有完全不同的内容，主要表现在以下两个方面。第4章 功率谱估计式(4.2.31)中自相关矩阵Rp是对称的半正定矩阵，且不具有Toeplitz性质，故式(4.2.31)的求解不能采用Levinson递归算法来进行，虽然现在已经提出一系列协方差法的求解算法，但总的来说，其算法求解的复杂性仍远远超过自相关法。另外由于Rp的非正定性，它也不能保证所得到的预测误差滤波器具有最小相位特性，换句话说，利用协方差法估计出的AR模型极点不能保证在单位圆内。举个例子，当p=1和N2
41、时，有(1,0)(0,1)(0)(1)xxrrxx而2(1,1)(0)xrx第4章 功率谱估计因此(1)(1)(0)xax 显然a(1)的幅度也可能大于或等于1。采用协方差法对信号进行建模，能够较好地反映出信号真正的模型。例例 4.2 采用与自相关法讨论中相同的例子，根据x(n)，n=0，1，N1的观察数据采用协方差法对其用一阶的AR模型建模。首先计算rx(1，0)和rx(1，1)，即22121111(1,0)(1)()111NNxnrrx nx nrNNr 22121111(1,1)(1)()111NNxnrrx nx nNNr第4章 功率谱估计由此得(1,0)(1)(1,1)xxrarr
42、由此可见，只要N2，所建立的模型就与信号的真正模型严格相同。同理可算出2211(0,0)11 1NxrrNr由此求得min(0,0)(1)(0,1)0 xxrar这说明一阶AR模型和数据的真正模型完全一致。第4章 功率谱估计3.修正的协方差法修正的协方差法 假定观察到的数据是x(0)，x(1)，x(N1)，分别写出前向预测与后向预测的表示式 1()()()pfkxna k x nk 1()()()pbkxna k x nk 和其中a(k)是AR模型系数，然后分别写出前向预测误差功率和后向预测误差功率的表示式2111()()()pNfn pkx na k x nkNp21011()()()Npp
43、bnkx na k x nkNp 第4章 功率谱估计则前向预测误差功率f 和后向预测误差功率b 的平均值 为(4.2.32)1()2fb修正的协方差法采用前向和后向预测误差功率的平均值为极小的方法估计AR参数。为使式(4.2.32)极小化，可以将 相对于a(l)求微分，并令其等于0，得111011()()()()2()()()()()()01,2,.,pNn pkNppnkx na k x nkx nlNpa lx na k x nkx nllp 第4章 功率谱估计经整理得1110110()()()()()()()()()pNpNkn pnNpNn pna kx nk x nlx nk x n
44、lx n x nlx n x nl 令1101(,)()()()()2()(,)NpNxn pnxrl kx nl x nkx nl x nkNprk l 式(4.2.34)可写成(4.2.35)第4章 功率谱估计(4.2.36)1()(,)(,0)pxxka k rl krl l=1，2，p将式(4.2.36)写成矩阵形式，即1,11,21,(1)(1,0)2,12,22,(2)(2,0),1,2,()(,0)xxxxxxxxxxxxrrrparrrrparrprprp pa prp 求出白噪声方差2的估计值2:12min11()()()()2()pNnpkx na k x nkx nNp1
45、01()()()()Nppnkx na k x nkx n 第4章 功率谱估计其中利用了式(4.2.36)的结果，最后得2min1(0,0)()(0,)pxxkra k rk将上式与式(4.2.36)合并得min10,00,10,(1)1,01,11,00(),0,1,xxxxxxxxxrrrparrrpa prprprp p (4.2.37)或写成minp pRa0第4章 功率谱估计修正的协方差法的一个优点是，它的误差功率的计算(见式(4.2.32)是在相对于协方差法多一倍的数据点上进行，这在观察数据长度很短的情况下，是非常有利的，但是这要求信号在正反两个方向上呈现相同的特性。例如，噪声中叠
46、加上多个正弦信号就具有这种特性，因为正弦信号在正反两个方向上看起来是相同的，而噪声的自相关函数也是与方向无关的。另一方面，若对于一个在两个方向呈现不同特性的信号，如一个指数衰减的信号，从相反的方向看，就是一个指数增长的信号，对这种信号如仍采用修正的协方差法来进行AR模型的估计，可能就会得到不好的结果。第4章 功率谱估计例例 4.3 采用与自相关法和协方差法讨论中相同的例子。根据x(n)，n=0，1，N1 的观察数据，采用修正的协方差法对其用一阶的AR模型建模。解解:首先计算 rx(1，0)和rx(1，1)，即22211(1,0)22(1)1NxrrrNr和222211(1,1)(1)2(1)1
47、NxrrrNr由此得2(1,0)2(1)1(1,1)xxrarrr 第4章 功率谱估计因为模型的真正解是a(1)r，这说明此模型参量的估计是有偏的，且与数据的长度N无关。因为r|a(1)|，因此一阶AR模型的极点比真正的极点更接近于单位圆。换句话说，用修正协方差法对一随机信号进行建模，可以得到一“锐化”的谱分析结果，因此用修正协方差法对随机信号建模，可以很好地用来进行高分辨率谱估计。式(4.2.37)中自相关矩阵Rp是一个非Toeplitz矩阵，另外，虽然Rp是正定的(纯正弦信号情况除外)，但它也不能保证所得到的预测误差滤波器具有最小相位特性，即利用修正协方差法所估计出的极点与协方差法一样，不
48、能保证在单位圆内。但是有一个例外，在一阶的情况下，第4章 功率谱估计p=1时用修正协方差法所建立的AR模型又总是稳定的，换句话说，所估计出的极点总是在单位圆内。这是因为在p=1时，由式(4.2.35)可知12101(1,0)(1)()(1)()2(1)NNxnnrx nx nx nx nN111(1)()2(1)Nnx nx nN和12101(1,1)(1)(1)(1)(1)2(1)NNxnnrx nx nx nx nN12211(1)()2(1)NnxnxnN第4章 功率谱估计由此得1112212(1)()(1,0)(1)(1,1)(1)()NxnNxnx nx nrarxnxn 因为在观察
49、数据x(n)不恒为0的情况下，上式的分母总是大于分子，因此|a(1)|p 时的rx(k)值。用re(k)表示外推的值，显然应对re(k)加一些约束。例如，若(4.3.1)jjj(e)()e()epwkwkwxxekpkpPr kr k则Px(ejw)应对应于一个合法的功率谱，即Px(ejw)对所有w都是实值的且非负的。但一般来说,只约束Px(ejw)是实值和非负的还不能保证获得唯一的外推结果，因此必须对可允许的外推再加上一些约束。Burg就提出一个这样的约束，就是使随机过程的熵达到最大化。第4章 功率谱估计由于熵是随机性或不确定性的一个测度，最大熵外推等价于要找出使x(n)尽可能“白”(随机)
50、的自相关序列re(k)。在某种意义上，这种外推对x(n)加了尽可能少的约束或最少量的结构要求。对功率谱而言，这对应于约束Px(ejw)“尽可能平坦”。如图4.8所示，不同的外推方法会得到不同的功率谱形状，其中最上面的外推方法有最“平坦”的谱。第4章 功率谱估计图 4.8 自相关序列的不同外推方法及其对应功率谱第4章 功率谱估计对一个功率谱为Px(ejw)的高斯随机过程，其熵是(4.3.2)j1()ln(e)d2wxH xPw因此对于高斯过程，若已知部分自相关序列rx(k)(|k|p)，其最大熵功率谱是使式(4.3.2)最大化,同时约束Px(ejw)的逆DTFT在|k|p时等于给定的自相关值，即
