第四章 功率谱估计

1、1 对于确定性信号，傅里叶变换是在频率分析研究的理论基础，但对于随机信号，其傅里叶变换不存在，因此研究它的功率谱。第四章功率谱估计4.1 引言 在实际应用中，通常只能采集或观测到平稳随机过程的一个抽样序列的一段(有限个)数据，如果根据这有限个已知数据来估计随机过程的功率谱的问题，简称谱估计(谱分析)问题。对信号和系统进行分析研究、处理有两类方法：一类是时域进行，维纳-卡尔曼滤波和自适应滤波都属于时域处理方法；另一类是频域方法。2按照weiner-khintchine定理，信号的功率谱和其自相关函数服从一对傅里叶变换关系，称为功率谱的定义。()( )jj mxxxxmp erm e1( )()2

2、jj mxxxxrmp eed*( )( ) ()xxrme x n x nm对于平稳随机信号，服从各态历经性，集合平均可以用时间平均代替，可推出功率谱的另一个定义。(4.1.1)(4.1.2)(4.1.3)将(4.1.3)式中的集合平均用时间平均代替，得*1( )lim( ) ()21nxxnnnrmx n x nmn(4.1.4)3将(4.1.4)式代入(4.1.1)式，得到()( )jj mxxxxmp erm e*1lim( ) ()21nj mnmnnx n x nm en*()1lim( )()21nj njn mnnnmx n ex nm en令 ，则lnm*1()lim( )

3、( )21njj nj lxxnnnlp ex n ex l en21lim( )21nj nnnnx n en 是随机变量，必须对 取统计平均值，得到()jxxp e()jxxp e21()lim( )21njj nxxnnnp eex n en将(4.1.4)式代入(4.1.1)式，得到上式被认为是功率谱的另一定义。(4.1.5)(4.1.6)4谱估计方法经典谱估计，也称为线性谱估计现代谱估计，也称为非线性谱估计经典谱估计bt法：1958年，r.blackmant和j.tukey提出， 先估计自相关函数，再计算功率谱。周期图法：1898年，schuster利用傅里叶级数去拟合待分析的信号，

4、提出周期图的术语，但直到fft出现，周期图法才受到人们的重视。这种方法直接对观测数据进行fft，取模平方，除以n得到功率谱。经典谱估计致命的缺点是频率分辨率低，原因是傅里叶变换域是无限大，而观测数据是有限长，观测不到的数据被认为是0。这相当于将信号在时域加了矩形窗，在频域使真正的功率谱卷积一个sinc函数。5现代谱估计( )h z( )w n( )x n22()()jjxxwpeh e如果由观测数据能够估计出信号模型的参数，信号的功率谱可以计算出来。谱估计问题变成了由观测数据估计信号模型参数的问题。模型种类很多，如ar模型、ma模型等。合适地选择模型，功率谱估计质量比经典谱估计的估计质量有很大

5、提高。参数模型法以信号模型为基础。参数模型法：ar、ma、arma模型非参数模型法：pisarenko、music、esprit法6图a bt法图b 最大熵谱估计法图c pisarenko谐波分解法(自回归psd法)71. 无偏自相关函数的估计4.2 经典谱估计4.2.1bt法bt法是先估计自相关函数，再进行傅里叶变换得到功率谱。下面介绍自相关函数的估计。利用观测到的实随机序列 ，估计自相关函数的两种方法是：无偏自相关函数估计和有偏自相关函数估计。( )(0,1,1)x n nn01mn估计器为( )()xxxxrmrm101 ( )( ) ()nmxxnrmx n x nmnm10nm8也可

6、写成一个表达式11nmn101 ( )( ) ()nmxxnrmx n x nmnm估计性能分析：估计量的偏差： ( )( )xxxxbe rmrm101 ( ) ( ) ()( )nmxxxxne rme x n x nmrmnm 0b，这是一种无偏估计。估计量的方差：222var( )( )( ) ( )( )xxxxxxxxxxrme rme rme r mrm9可以证明，12211var( )( )()()()nmxxxxxxxxkmnrmrkrkm rkmnmknm nmk一般观测数据量n很大，(1)mknmknnn1221var( )( )()()()nmxxxxxxxxkmnnr

7、mrkrkm rkmnm 此式表明，只有当 时，估计量的方差才趋于0，此估计是渐近一致估计。但是当 时，方差将很大，因此，这种方法在一般情况下不是一种好的估计方法。虽然是无偏的，但不能算是一致的。,nm n mn102. 有偏自相关函数的估计有偏自相关函数用 表示，估计器为 ( )xxrm11nmn101 ( )( ) ()nmxxnrmx n x nmn估计性能分析：估计量的偏差：( )( )( )( )xxxxxxxxnmnme rme rmrmrmnn因为 是无偏估计，两边取均值，得( )( )xxxxnmrmrmn ( )xxrm有偏估计11由 可见，只有当 时， 才是无偏的，其它 都

8、是有偏的，但当 时， ，因此 是渐近无偏的。偏移量为 ( )( )xxxxnme rmrmn( )( )( )xxxxxxmbrme rmrmn0m ( )xxrmmn ( )xxrm0b ( )xxrm估计量的方差：由两种估计的关系式( )( )xxxxnmrmrmn知估计量 的方差为122var( )() var( )xxxxnmrmrmn当 时， ，并且2var( )() var( )xxxxnmrmrmn21221var( )( )()()()nmxxxxxxxxkmnnmnrmrkrkm rkmnnm var( )var( )xxxxrmrmn var( )0 xxrm 虽然是有偏估

9、计，但是渐近一致估计，估计量的方差小于 的方差。实际应用中多用这种有偏自相关估计。也用符号 表示。 ( )xxrm ( )xxrm ( )xxrm13bt法功率谱估计采用有偏自相关函数估计法，101 ()( ) ()nmxxnrmx n x nmn对上式进行傅里叶变换，得到bt法的功率谱估计值为()()jj mbtxxmperm e为了减少谱估计的方差，经常用窗函数 对自相关函数进行加权，此时谱估计公式为( )w n1(1)()() ()mjj mbtxxmmperm w m e式中() (1)(1)( )mn0 w mmmmw n，其它加权协方差谱估计要求加窗后的功率谱仍是非负的，这样窗函数

10、的选择必须满足一个原则，即它的傅里叶变换必须是非负的，如巴特利特窗14为了采用fft计算 ，设fft的变换域为 ，必须将求和域 移到 ，功率谱的计算公式如下：1(1)()() ()mjj mbtxxmmperm w m e10()()ljj mbtxxmpesm e(1,1)mm(0 1)l (0 1)l ( )()btxxpkfft sm0,1,1kl () () 01()0 1 () () 1xxxxxxrm w mmmsmmmlmrml w mllmml这里自相关函数是渐近一致估计，但经过傅里叶变换得到功率谱的估计却不一定是渐近一致估计，可以证明它是非一致估计，是一种不好的估计方法。15

11、4.2.2周期图法功率谱的另一种定义式重写如下：21()lim( )21njj nxxnnnp eex n en忽略上式中求统计平均的运算，观测数据为 ，便得到周期图法的定义：( )(01)x nnn2101()( )njj nxxnp ex n en上式的绝对值符号内的部分可以用fft计算，这样周期图法的计算框图如下图所示。16下面证明周期图和bt法的等价关系，再分析周期图的估计质量。1.周期图和bt法的等价关系2101()( )njj nxxnp ex n en11*001( )( )nnj kj nknx k ex n en11*()001( )( )nnjk nknx k x n en

12、 令 ，即 ，则mknkmn11*(1)01()( ) ()nnmjj mxxmnnp ex n x mnen 上式中的方括号部分正是有偏自相关函数的计算公式，则1(1)()( )njj mxxxxmnp erm e()()jjxxbtp epe等价17已知自相关函数的估计值 ， ，求功率谱的统计平均值，(1)1nmn2.周期图法谱估计质量分析1)周期图的偏移 ( )xxrm1(1)()( )njj mxxxxmnp erm e1(1)()( )njj mxxxxmne p ee rm e1(1)( )nj mxxmnnmrm en( )( )j mbxxmw m rm e18 称为三角谱窗函

13、数。 式表明，周期图的统计平均值等于它的真值卷积三角谱窗函数，因此周期图是有偏估计，但当 时， ，三角谱窗函数趋近于 函数，周期图的统计平均值趋于它的真值，因此周期图属于渐近无偏估计。式 中乘积的傅里叶变换，在频域服从卷积关系，得到 ( )0 bnmmnw mn其它式中()-1()()()2jjjxxbxxe p ew ep ed式中()( )jxxxxp eft rm21sin(/2)()( )sin(/2)jbbnw eft w mn()jbw e( )1bw m n 192)周期图的方差 由于周期图方差的精确表示式很复杂，为分析简单，通常假设 是实的零均值的正态白噪声信号，方差是 ，即功

14、率谱是常数 ，其周期图用 表示， 表示观测数据的长度，则( )x n2x2x( )nin2()11( )() ( ) ( )( )( )jjn knnnnke ix ee x k x n rk rn enn 22var( )( )( )nnnie ie i下面先求周期图的均值，再求其均方值：1120011( )()( ) ( )nnjj kj nnnkix ex k x n eenn 20式中2 ( ) ( )()()xxxe x k x nrnknk 2221var( )( )nnxxnirnn1 ()0 nknknk12221221()()()() jjnnnne iiexexen上式说明

15、周期图是无偏估计，这是由于假设信号是实白噪声。求均方值时，先求两个频率 和 处的均方值，最后令121221 ( ) ( ) ( ) ( )nkpqe x k x n x p x qn 12()()( )( )( )( )jn kjp qnnnnrk rn rp rq ee21利用正态白噪声、多元正态随机变量多阶矩公式，有 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )e x k x n x p x qe x k x n e x p x q ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )e x n x p e x k x qe x n x q e x k x p4

16、 ,;,;, ( ) ( ) ( ) ( )0 xkn pq pn qk qn pke x k x n x p x q其它将上式代入周期图的均方值公式中，得12124()()()()2122()()=j n kj n kxnnnknke iineen 22412121212sin()/2)sin()/2)=1sin()/2sin()/2xnnnn22将 代入上式，得12224sin()( )=2sin( )nxne in2224sin()var( )= ( )-( )=1sin( )nnnxnie ie in显然，当n趋于无限大时，周期图的方差并不趋于0，而是趋于功率谱真值的平方，即4var(

17、 )nxni无论n取多长，周期图的方差总是和 同一数量级。信号功率谱真值是 ，说明周期图的方差很大，周期图的均方误差也非常大。用这种方法估计的功率谱在 附近起伏很大，故周期图是非一致估计，是一种很差的功率谱估计方法。显然，当n趋于无限大时，周期图的方差并不趋于0，而是趋于功率谱真值的平方，即4var( )nxni4x2x2x23图 白噪声的周期图24则第 组的周期图为： 假设随机信号 的观测数据区间为： ，共进行l次独立观测，得到l组记录数据，表示为 ，4.2.2经典谱估计方法改进一般有三种改进方法：平均周期图法、窗函数法和修正的周期图平均法。1.平均周期图法平均周期图的思想平均周期图的思想：

18、对一个随机变量进行观测，得到l组独立记录数据，用每一组数据求其均值，然后将l个均值加起来求平均。这样得到的均值，其方差将是用一组数据得到的均值的方差的1/l。i2101( )( )mj niinix n em01nm( )ix n1,2,il( )x n25可见，平均周期图仍然是有偏估计，偏移和每一段数据的个数m有关；由于 ，平均周期图的偏移比周期图的偏移大。由于平均周期图三角谱窗主瓣的宽度变宽，分辨率更加降低，因此，偏移的大小反映分辨率的低与高。偏移越大分辨率越低。将得到的l个周期图进行平均，作为信号 的功率谱估计，有mn( )x n11()( )ljxxiip eil为分析偏移，对上式求统

19、计平均，得11() ( ) ( )ljxxiiie p ee ie il前边已求出周期图的统计平均值为()-1()()()2jjjxxbxxe p ew ep ed21sin(/2)()( )sin(/2)jbbnw eft w mn26但实际中数据是连续的，数据段之间不能认为是完全不相关的，估计方差的减少小于 ，或认为近似为原来的 。平均周期图的方差是周期图的方差的 。显然是以分辨率的降低换取了估计方差的减少，当然，估计的均方误差也减少。 求平均周期图的方差，由于是l次独立观测，l个周期图相互独立，因此平均周期图的方差为1/ l1var()var ( )jxxip eil实际中，很难得到独立

20、的数据组，一般常用的方法是将长度为n的数据分成l段，每段有m个数据，n=lm，第i段数据表示为：1il 。( )(),ix nx nimm01,nm如果各段数据互不相关，利用平均周期图法估计功率谱的公式、偏移和方差的计算公式和上边结果一样，方差是周期图法的1/ l1/ l1/ l27按照 式，段数l愈多，方差愈小，即功率谱愈平滑，如果 ，则 。但如果数据量n一定，l加大，每一段的数据量m将减少，因此估计量方差减少了，使偏移加大，分辨率降低。估计量的方差和分辨率是一对矛盾。1var()var ( )jxxip eill var()0jxxp e可根据实际情况适当选择l和m，如果对分辨率要求不高，

21、可以取l大些。最好将n取大些，分辨率和估计误差都能适当满足要求。下图中，信号是均值为0，方差为1的白噪声，观察数据长度为n=256，l=2,4,8段，按平均周期图法估计功率谱曲线如图所示。28表明，随着分段数的增加，功率谱估计值在1附近摆动的幅度越来越小，显示出分段平均对周期图方差减少有明显效果。图 平均周期图法292.窗函数法 这种方法是用一适当的功率谱窗函数 与周期图进行卷积，来达到使周期图平滑的目的。即1( )()2jj nw nw eed()jw e()-1()( )()2jjxxnp eiw ed 1(1)( )( )nj mnxxmnirm e式中(1)1mnm 是有偏自相关函数。

22、 ( )xxrm那么1(1)()( ) ( )mjj mxxxxmmp erm w m e此式和前面讲的bt法的加权协方差谱估计方法是相同的。30 两个窗函数有关。将 和 代入上式，得在 式中，周期图和谱窗函数卷积得到功率谱，等效于在频域对周期图进行修正，使周期图通过一个线性非频变系统，滤除掉周期图中的快变成分，谱窗函数需具有低通特性。()-1()( )()2jjxxnp eiw ed 求 的统计平均，得()jxxp e()mn1(1)()( ) ( )mjj mxxxxmme p ee rm w m e ( )( )xxxxnme rmrmn( )bnmw mn1(1)()( )( ) (

23、)mjj mxxxxbmme p erm w m w m e上式表明，周期图的窗函数法仍然是有偏估计，其偏移和 、 ( )bw m( )w m31比较上式和周期图法的 式，等于将 式中的窗函数 改变成 。 ，则如果 窗的宽度比较窄，m比n小得多， ，则( ) 1bw m( )w mmn1(1)()( ) ( )mjj mxxxxmme p erm w m e()-1()()2jjxxp ew ed ()-1()()()2jjjxxbxxe p ew ep ed()jw e()jbw e如果 比 的主瓣宽度宽，利用窗函数法可进一步平滑周期图，减少估计方差，但是偏移加大，分辨率降低。()jw e(

24、)jbw e由于功率谱是频率的非负函数，要求 是非负函数，三角窗函数是适用的，哈明窗、汉宁窗并不适用。()jw e323.修正的周期图求平均法 这种方法和平均周期图法一样，首先把数据长度为n的信号 分成l段，每一段数据长度为m，n=lm，然后把窗函数 加到每一个数据段上，求出每一段的周期图，形成修正的周期图，再对每一个修正的周期图进行平均。( )x n2101( )( ) ( )mj niinix n w n eu( )w n第i段的修正周期图为1201( )mnuw nm式中u称为归一化因子。将每一段的修正周期图之间近似看成互不相关，最后功率谱估计为11()( )ljxxiip eil33对

25、上式求统计平均，得2101()( )mjj nnw ew n emu()-1()()()2jjjxxxxe p ep ew ed 式中 这种在计算周期图之前，先对各数据段加窗函数的方法，使平均周期图方法的估计方差减少，分辨率降低。但这种方法对窗函数没有限制；分段时，相邻的两段可以重叠，进一步使方差减少，可以重叠50%。为使功率谱估计渐近无偏，归一化因子u是不可缺少的。34由图可见，加哈明窗，大大压低了旁瓣，使低电平信号清晰可见，但由于主瓣加宽，功率谱波峰变宽了，降低了信号的分辨率。一般来说，两个等幅的正弦信号的频率相隔很近，可以不加数据窗，频率间隔应该大于 才能分辨。图 利用数据窗减少窄带过程

26、周期图的旁瓣(a)没有数据窗 (b)加哈明窗2 / n35 (1)选择合适的信号模型： (2)根据 有限的观测数据，或者它的有限个自相关函数估计值，估计模型的参数； (3)计算模型输出功率谱。 如果能确定信号 的信号模型，根据信号观测数据求出模型参数，系统函数用 表示，模型输入白噪声方差为 信号的功率谱可用下式求出22()()jjxxwp eh e4.3 现代谱估计中的参数建模( )x n( )h z2w功率谱估计可分成三个步骤： 这种方法隐藏着数据的外推，不是假设观测区以外的数据为0，它的估计分辨率都比经典谱估计高，故现代谱估计也统称为高分辨率谱估计。( )x n36 主要考虑原则：模型能够

27、表示谱峰、谱谷能力。对于具有尖峰的谱，应该选用具有极点的模型，如ar和arma模型；对于具有平坦的谱峰和深谷的信号，可以选用ma模型；既有极点也有零点的谱应该选用arma模型。对于没有模型准确反映其谱特性的信号，可以选用高阶的ar模型近似表示。 在选择模型合适基础上，应尽量减少模型的参数。4.3.1模型选择下图表示的是用ma模型估计二阶ar信号功率谱的例子。37图 对ar(2)信号的模型选择38假设模型的差分方程和系统函数分别用下式表示：11( )( )( )1qiiipiiibzb zh za za z4.3.2模型参数和自相关函数之间的关系11( )()()pqiiiix na x nib

28、w ni 式中，模型输入白噪声 均值为0，方差为 ，b(z)称为ma部分，是一q阶多项式，单位脉冲响应为 ， ， ( )w n( )bh n2w( )bnh nb0,nqa(z)称为ar部分，是一个p阶多项式， ， ，( )anh na0,np0(0)1aha39将上式两边同乘以a(z),得1.arma模型的系数和信号自相关函数之间的关系*1( )( )( )()xxwwpzpz h z hz*2*11( ) ( )( ) ( )()() ( )xxwwwpz a zpz b z hhb zzz将上式进行z反变换，左边的反变换为10( ) ( )( )( )( )()pxxwwaawwlz t

29、 pz a zrmhmh l rml公式右边的z反变换为12*2*01() ( )( ) ()qwwblz thb zh l q mlz式中1*1( )()()q mz t hhmz40此式即是arma模型输出自相关函数与模型参数之间的关系式。如果能由信号的观测数据估计出信号的自相关函数，即可求出arma模型的参数，再由 求出信号功率谱。由于模型h(z)是因果的， ， ，可以得到代入上式得12*2*01() ( )( )()qwwblz thb zh l hmlz( )0h n 这样得到2*00( )()( )()pqawwwbllh l rmlh l hml0n 2*101( )()( )(

30、) 0,1,( )( )() 1pq mawwwbllxxpawwlh l rmlh l h lmmqrmh l rmlmq2()()()jjjxxwwpepeh e41 ， ，然后再设法求出ma部分的系数。当 时，上式是一个线性方程，用矩阵表示如下( )(1)(1)(1)(1)(1)( )(2)(2)(2)(1)(2)( )()()xxxxxxxxaxxxxxxxxaxxxxxxxxarqrqrqprqhrqrqrqprqhrqprqprqrqphqp ( )ah imq上式共有p个方程。可以用该方程首先计算出ar部分的p个系数1,2,ip42ar模型的系统函数 ，相当于arma模型中 的情

31、况，这样公式中2.ar模型的系数和信号自相关函数之间的关系( )1/( )h za z( )1b z ( )( )bh nn因为 是因果的，因此 时 ，得1 ()()0 blmh lmlmlm *()0hm0m ( )h n121( )() 1( )( )() 0pawwlxxpawwwlh l rmlmrmh l rmlm43 式的矩阵形式如下也可以将上式中 的情况写成矩阵形式：1m (0)( 1)(1)(1)(1)(1)(0)(2)(2)(2)(1)(2)(0)( )( )xxxxxxxxaxxxxxxxxaxxxxxxxxarrrprhrrrprhrprprrphp 2(0)( 1)()

32、1(1)(0)(1)(1)0( )(1)(0)( )0 xxxxxxwxxxxxxaxxxxxxarrrprrrphrprprhp44 称为自相关矩阵，满足 ，是一个hermitian矩阵，也是一具toeplitz矩阵，是正定矩阵。令或者用模型参数表示211(0)( 1)()(1)(0)(1)0( )(1)(0)0 xxxxxxwxxxxxxpxxxxxxrrrparrrparprpr(1)(0)( 1)()(1)(0)(1)( )(1)(0)xxxxxxxxxxxxpxxxxxxxxrrrprrrprprprr(1)pxxrhxxxxrr45 的情况，此时， ，这样可得ma模型系数和信号相关

33、函数的关系为3.ma模型的系数和信号自相关函数之间的关系ma模型的系统函数 ，相当于arma模型中 ，( )( )h zb z( )1a z ( )( )ahnn( )( )bh nhn2*0( )() 0,1,( )0 1q mwblxxh l h lmmqrmmq222211()()1qjijjixxwwpjiibep eh eae464.4 ar模型谱估计的性质4.4.1ar的线性预测xxxxxx(0) ( 1) r ()(1) (0) r (1) ( ) (1) r (0)xxxxpxxxxxxxxrrparrprprp12min1( )00ppe e na 第二章，已推导出维纳线性一

34、步预测器系数和信号自相关函数之间的关系式，这就是著名的yule-walker方程。式中， 表示线性一步预测误差，公式为1( )( )( )( )()ppiie nx nx nx na x ni( )e n 表示预测器的系数，它和线性预测器的 关系为pia( )h n,( ) 1,2,pnh nanp 47令 ，由上式，得对 式进行z变换，得1( )( )()ppiie nx na x ni( )( )/( )ehze zx z1( )( )( )pipiie zx zx z a z1( )1piepiihza z 称 为线性一步预测误差滤波器，其作用是将 转换成预测误差 ，如图所示。一般认为

35、具有白噪声的性质，因此 也称为白化滤波器。( )ehz( )x n( )e n( )e n( )ehz图 预测误差滤波器48对比ar模型和预测误差滤波器的系统函数1( )1piepiihza z ( )( )( )x nx nw n11( )1piiih za z当 时， 和 互为逆滤波器， 。基于以上分析，可得ar模型定义为,(1,2, )piiaaip( )1/( )ehzh z( )h z( )ehz式中， 是基于信号前p个样本的最佳一步线性预测，为( )x n1( )()ppiix na x ni 是模型输入白噪声。( )w n49对比21(0)( 1)()(1)(0)(1)0( )(

36、1)(0)0 xxxxxxwxxxxxxpxxxxxxrrrparrrparprprxxxxxx(0) ( 1) r ()(1) (0) r (1) ( ) (1) r (0)xxxxpxxxxxxxxrrparrprprp12min1( )00ppe e na当ar模型的阶数与线性预测器的阶数相同时，得到piiaa( )( )w ne n22min( )we e n信号自相关函数和它的ar模型参数之间的关系服从yule-walker方程。因此ar模型法也称为线性预测ar模型法。50如果 具有白噪声性质，可以利用ar模型与预测滤波器之间的关系对上式进行解卷积，得到 信号。方法是：先由 的观测数

37、据估计ar模型的参数，得到ar模型的系统函数 ，再让 通过 的逆滤波器 ，便得到信号 。ar模型与线性预测之间的关系，可以被用来解卷积，假设信号 通过一个ar系统，系统单位脉冲响应 ，响应 ，( )s n( )h n( )x n( )( )( )x ns nh n( )s n( )s n( )x n( )h z( )x n( )h z1( )hz( )s n51 ar模型必须因果稳定，即极点均在单位圆内，才能保证信号 是平稳随机信号，于是ar模型 和预测误差滤波器 互为逆滤波，那么 应为最小相位系统。 但是由解yule-walker方程得到ar模型的参数，其极点不一定在单位圆内。 下面证明当最

38、佳p阶线性预测系数与ar模型参数相同时，由此得到的极点保证在单位圆内，ar滤波器稳定，预测误差滤波器 或者 是最小相位系统。 4.4.2预测误差滤波器的最小相位系统( )x n( )h z( )ehz( )ehz( )ehz( )a z52设 的第i个零点 在单位圆外部，即 ，用 代替 ，这时 的幅度函数不变，按照上式计算出的观测误差滤波器输出功率 应不变，仍是最小。下面将用证明预测误差滤波器输出功率 不是最小的，这一矛盾的结论只能说明 的零点不可能在单位圆外部。式中 ， 是最佳预测系数。 解yule-walker方程得到的是最佳线性预测滤波器的系数，此时预测误差滤波器输出功率 达到最小，用

39、表示，即epminep1( )1piiia za z 2min()()jjexxpa eped 下面用反证法证明 的全部零点不在单位圆外部。iaiz( )a z( )a z1iz *1/iziz( )a zep( )a z53说明将零点 换成 时，预测误差滤波器的输出功率将更小,与假设矛盾。 零点必须在单位圆内或单位圆上。又因为自相关矩阵是正定矩阵， 的零点必须在单位圆内部。因为 ，则1iz 111( )1(1)ppiijija za zz z 1111(1)(1) (1)( )pijijj iz zz zz za z式中 ，代入 ，得11( )(1)pjjj ia zz z2min()()j

40、jexxpa eped22min(1)( )()jjeixxpz ea zped式中2222222*11111jjjjiiiiiiiz ezezezezzz22*111jjiiz eez*1/iziz( )a z( )a z54对于 情况，公式本身就是一个递推方程，如果已由观测数据计算出p+1个自相关函数，用 ， 表示，对于 的情况，可以用该公式外推得到，ar模型的自相关函数和模型系数服从yule-walker方程4.4.3 ar模型隐含自相关函数延拓特性1m 121( )() 1( )( )() 0paxxlxxpaxxwlh l rmlmrmh l rmlm ( )xxrm0,1,mpmp

41、1 ( ) 0 ( )( )() xxpxxaxxlrmmprmh l rmlmpar模型的自相关函数隐含着自相关函数外推的特性，它具有很高分辨率。bt法中自相关函数只由观测数据计算出有限个值，其它认为是0。554.5 ar谱估计的方法 已知信号的观测数据，估计功率谱的方法：自相关法、伯格法、协方差和修正协方差法。 它们共同的特点是由信号观测数据求模型系数时，均使用信号预测误差最小的原则。 对于长记录数据，这些方法的估计质量相似；但对于短记录数据，不同的方法之间存在差别。4.5.1 自相关法-列文森(levinson)递推法自相关法的出发点是选择ar模型的参数使预测误差功率最小。56假设信号

42、的数据区在 范围，有 个预测系数， 个数据经过冲激响应为 的滤波器，输出预测误差 的长度为 ，因此有预测误差功率为( )x n211200111( )( )()npnpppinnie nx na x ninn ( )e n01nnnppn(0,1, )piaip22111( )( )()ppinnie nx na x ninn 的长度长于数据的长度，上式中数据 在 以外补充零点，相当于对无穷长的信号加窗处理，会引入误差。( )e n( )x n01nn57上式对系数 的实部和虚部求微分使预测误差功率最小，得pia12(0)( 1)(1)(1)(1)(0)(2)(2)(1)(2)(0)( )px

43、xxxxxxxpxxxxxxxxppxxxxxxxxarrrprarrrprarprprrp 式中自相关函数采用有偏自相关估计，即1*0*1( ) () 0,1, ( ) () 1,2, 1nmnxxxxx n x nmmprmnrmmpp yule-walker方程58levinson-durbin算法的推导：最直观的方法是根据ar(1)、 ar(2)、 ar(3)各阶模型的yule-walker方程的求解结果归纳出一般的迭代计算公式。第二种方法是利用自相关矩阵对称性和toeplitz性质，通过构造增广矩阵导出。 用线性方程组的常用解法(如高斯消元法)求解yule-walker方程，运算量数

44、量级为 。利用自相关矩阵的对称性和toeplitz性质，可构成一些高效算法(levinson-durbin算法)，运算量数量级为 。3p2plevinson-durbin算法：是一种按阶次递推的算法。它以ar(0)和ar(1)模型参数作为初始条件，计算ar(2)模型参数；再用ar(2)模型参数计算ar(3)模型参数，k阶模型参数由k-1阶模型参数计算得到。一直计算出ar(p)模型参数为止。59一阶ar模型(p=1)的yule-walker方程为211,11(0)(1)(1)(0)0 xxxxxxxxrrarr由该方程解出1,1(1)(0)xxxxrar 2211,1(1)(0)xxarxx,1

45、xxxx(0) (1) r ( )(1) (0) r (1) ( ) (1) r (0)xxxxpxxxxxxxxrrparrprprp12,00pp pa实平稳信号60由上面方程解出增加一阶，即令p=2，得222,12,2(0)(1)(2)1(1)(0)(1)0(2)(1)(0)0 xxxxxxxxxxxxxxxxxxrrrrrrarrra222,11,12,21,1(0)(1)(1)(2)/(0)(1)xxxxxxxxxxxxarrrrrraaa 22222,21(1)a21,11(2)(1)/xxxxra r 2222,2(0)(2)(1)/(0)(1)xxxxxxxxxxarrrrr

46、61 称为反射系数， 。 ，随着阶数增加，预测误差功率将减少或不变。然后令p=3,4,.,以此类推，可以得到一般递推公式如下：,1,1, 1,2,1p kpkppp kaak akp2221(1)pppk11,121( )()pxxpk xxkpprparpkk ,pp pka220(0)( )xxre xnpk1pk 22212p62由k=1开始递推，递推到k=p，依次得到各阶模型的参数， 递推公式提供了一种确定模型阶数的实验方法，如模型的阶数不知道，由低阶开始递推，当递推到m阶时，预测误差满足允许的值，停止递推，选ar模型的阶数为m。 这种方法递推效率高，当阶数变化时，无需从头计算。但需要

47、预先估计出信号自相关函数，当观测数据长度较短时，估计误差较大，会出现谱峰频率偏移和谱线分裂；如数据很长，估计自相关函数较准确。计算流程图如下2221112122212,pppppaaaaaaar模型的各个系数及模型输入白噪声方差求出后，信号功率谱用下式计算：22221()()1/(1)pjjj ixxwwpiipeh ea e63图 利用levinson-durbin递推法计算功率谱的流程图641. 协方差法4.5.2协方差法与修正协方差法 与自相关法一样，仍然利用使预测误差功率最小的方法求模型参数。但由观测数据求预测误差功率的公式不同，为2112111( )( )()pnnpkn pn pk

48、e nx na x nknpnp 与自相关法相比，求预测误差功率的求和限不同。该公式中使用的观测数据均已知，不需要在数据两端补充零点，因此比较自相关法去掉了加窗处理的不合理假设。65为求模型参数仍然应用复梯度法使上式达到最小，公式如下11( , )() ()nxxn pcj kx nj x nknpxxxxxx(1,1) (1,2) c (1, )(2,1) (2,2) c (2, ) ( ,1) ( ,2) c ( , )xxxxxxxxxxxxccpccpcpcpp p12(1,0)(2,0)( ,0)pxxpxxppxxacacacp 式中白噪声的方差为12min(0,0)(0, )nw

49、xxpkxxkpca ck称为协方差函数66 ，但不具有topelitz性质( )上式由 组成的 维矩阵是对称的，即xxmin121(0,0)(0,1)(0,2)c (0, )(1,0)(1,1)(1,2)(1, )0(2,0)(2,1)(2,2)(2, )0( ,0)( ,1)( ,2)( , )0 xxxxxxpxxxxxxxxpxxxxxxxxppxxxxxxxxcccpaccccpaccccpacpcpcpcp p ( , )( , )xxxxci jcj i(1) (1)pp( , )xxci j(1,1)( , )xxxxcijci j所以协方差方程的解不能借助于自相关法，可用乔里

50、斯基(choleskey)分解法求解上述方程组。用协方差法估计出的极点不能保证在单位圆内。67 和前向预测误差 用下式表示： 修正协方差方法使用前向和后向观测误差平均值最小的方法，估计ar模型参数，进而估计信号的功率谱。 该方法最初由nuttall在1976年提出，称为前向后向法；同年，ulrych和ctayton也独立提出，称之为最小二乘法，因此该法也称为前后向线性预测误差功率最小的最小二乘法。2. 修正协方差法下面对前向和后向线性预测误差滤波器进行分析。前向线性预测误差滤波器 为分析简单，假设信号是实平稳随机信号。前向线性预测误差滤波器由 预测 。其估计值(1), (2), ()x nx

51、nx np( )fpen( )x n( )x n68由于假设信号是实的，式中预测误差和系数均是实的。1( )()ppkkx na x nk 1( )( )( )( )()pfppkkenx nx nx na x nk对上式进行z变换，得到1( )( )( )pfkppkkezx za x z z令( )( )( )fpfezhzx z则010( )1 1ppkkfpkpkpkkhza za za 称为前向预测误差滤波器的系统函数( )fhz69将 代入，得用均方误差最小准则求前向预测误差滤波器的最佳系数图 前向预测误差滤波器pka2( ) 0 1,2,fppke enkpa1( )( )()p

52、fppkkenx na x nk( ) ()0 1,2,fpe en x nkkp前向预测误差与用于预测的数据正交，正交原理。70式中 。已知自相关函数，可解出 和 。前向预测误差滤波器的最佳系数 和信号的自相关函数之间的关系为yule-walker方程，如下pka121( )()0 1,2,(0)( )pxxpi xxippxxpi xxirka rkikpra ri用矩阵表示为22min( ) fppe enxx,1xxxx(0) (1) r ( )(1) (0) r (1) ( ) (1) r (0)xxxxpxxxxxxxxrrparrprprp12,00pp papka2p71 进行

53、预测，为此，将上式改为后向线性预测误差滤波器利用 数据预测 ，则称为后向线性预测。其估计值为(1), (2), ()x nx nx np( )x n1 ( )()ppkkx na x nk 一般前向、后向预测用同一数据进行，即利用( ),x n(1), (2), ()x nx nx np1 ()()ppkkx npa x npk 这样，前向预测是由 预测 ，后向预测是由 预测(), (1), (2), (1)x np x npx nx n()x np( )x n(1), (2), ( )x npx npx n预测 ，这两种预测数据之间的关系如图。72后向预测误差 表示的是信号在 时刻的预测误差

54、，图 前后向预测数据之间的关系( )bpennp( )()()bpenx npx np同样，利用最小均方准则，可得关于后向预测时的正交原理及yule-walker方程73式中， 是后向预测误差的最小误差功率。利用toeplitz性质，可得以下重要关系：121( )()0 1,2,(0)( )pxxpi xxippxxpi xxirka rkikpra ri( ) ()0 1,2,bpe en x npkkp 1,2, ,pkpkaakp22pp上式表明，前、后向预测的最小误差功率相等，系数也相等(如果是复数，则是共轭关系)。74式中，当 时， ，上式可写成由 ， ，得1 ()()ppkkx n

55、pa x npk ,1,0pkp p ( )()()bpenx npx np0,1,2,kp00( )() 1pbppkpkena x npka,00( )() 1pbpp p kpkenax nka图 后向预测误差滤波器后向预测误差滤波器的系数排序是前向预测误差滤波器系数排序的逆转排列。75对 进行z变换，得0( )()pbppkkena x npk1( )( )( )pbppkppkkezx z za x z zz后向预测误差滤波器的系统函数为1( )( )(1)( )bpppkbpkkezhzza zx z前、后向预测误差滤波器的系统函数之间关系为1( )()pbfhzzhz76前面已推

56、导出信号的前向和后向预测，分别为1( )()ppkkx na x nk 1( )()ppkkx na x nk 式中 是ar模型的参数。在这两种情况下，最小预测误差功率都是白噪声方差 。pka2w修正协方差使用估计的前向和后向预测误差功率的平均值为最小的方法来估计ar参数。即1()2ppfpb77式中前向和后向预测误差功率 、 分别用下式表示2111( )()pnpfpkn pkx na x nknppfpb21011( )()npppbpknkx na x nknp 和协方差法一样，仅对那些用到观测数据样本的预测误差求和。111( ( )() ()pnppkn pkplx na x nkx

57、nlanp101( ( )() ()0npppknkx na x nkx nl 1,2,lp78经过简化，得到令将上式写成矩阵形式1110() ()() ()pnpnpkkn pnax nk x nlx nk x nl 110( ) ()( ) ()npnn pnx n x nlx n x nl 1101( , )() ()() ()2()npnxxn pncj kx nj x nkx nj x nknp 79白噪声的方差估计值为xxxxxx(1,1) (1,2) c (1, )(2,1) (2,2) c (2, ) ( ,1) ( ,2) c ( , )xxxxxxxxxxxxccpccpc

58、pcpp p12(1,0)(2,0)( ,0)pxxpxxppxxacacacp 12,min(0,0)(0, )nwpxxpkxxkpca ck协方差法和修正协方差法相同，除了 的定义不同。修正协方差法还可以得到高分辨率的统计稳定的谱估计值。对于ar谱估计的谱峰由于附加观测噪声而偏离真实频率位置的现象很少。这种方法不能保证 ，因而不能保证滤波器的稳定性。( , )xxcj k1(1)kkakp80例：已知信号的四个观测数据为( ) (0), (1), (2), (3)x nxxxx ，分别用自相关法和协方差法估计ar(1)模型参数。2,4,1,3解：(1)自相关法14220011( )( )

59、4npnne ne nn 111( )( )()( )(1)ppiie nx na x nix na x n440011111( )1( )( ) (1)022nne ne ne n x naa111111112(42)4(14)(3)90aaaa110.5a 81(2)协方差法1322111( )( )3nn pne ne nnp111( )( )()( )(1)ppiie nx na x nix na x n330011112( )2( )( ) (1)033nne ne ne n x naa1111112(42)4(14)(3)0aaa110.714a 82levinson-durbin

60、递推法需要先由观测数据估计自相关函数，这是它的缺点。而伯格递推法则由信号观测数据直接计算ar模型的参数。4.5.3伯格(burg)递推法levinson-durbin递推法重写如下,1,1, 1,2,1p kpkppp kaak akp2221(1)pppk11,121( )()pxxpk xxkpprparpkk ,pp pka220(0)( )xxre xn,pp pka83 伯格递推法利用levinson-durbin递推公式，导出前向预测误差和后向预测误差，并按照使它们最小的原则求出 ，从而实现不用估计自相关函数，直接用观测数据得出结果。pk首先由前向、后向预测误差滤波器推导预测误差格