22_AR谱估计的方法

§4.5AR谱估计的方法引言（1）主要方法本节要解决的问题是:如果不知道的自相关函数,而只知道它的N个观测数据,,如何求模型的参数?通常有下列两类方法:a)利用求的估计,再用Levinson算法求模型参数的估计值;b)利用最小二乘方准则,直接由信号时间序列计算模型参数.本节介绍以下三种方法:—自相关法(Levinson法,或Yule-Walker法);—协方差法及修正协方差法;—Burg法.)(nx)(mRxx)(nx1,,1,0NnAR)(nx)(mRxx)(ˆmRxx)(nxAR第一类方法:由先估计自相关函数——第二类方法:先估计反射系数以上方法都是基于线性预测原理,根据估计准则建立求解AR系数的数学模型.(2)估计(即模型参数)的准则理论上,采用预测误差功率最小准则,表示为(4.5.1)实际上,采用由时间平均代替集合平均的最小平方准则,即(4.5.2)式中的求和范围暂时未指定,取决于下面要讨论的算法.)(mRxxmin})]({[2neEmin)(2nenn4.5.1自相关法——Levinson递推法(1)估计准则假设信号的数据区间为,前向预测误差滤波器的冲激响应为,则预测误差为,(4.5.3)上式表明,是与的卷积,其长度为.因此,式(4.5.2)可具体表示为(4.5.4)(2)预测误差的计算式(4.5.3)的计算原理如图4.5.1所示.由于的长度大于的长度,因此计算中需对的两端补零.这相当于长度为点的已知数据,是由无限长数据()经加窗后得到的.)(nx10Nnpipiah),,1,0(pipipipipiinxainxhnxne01)()()()(10pa)(ne)(nxpiapNmin)(1|)(|12100102pNnpipipNninxaNneN)(ne)(ne)(nx)(nxN)(nx)(nxn用时间平均代替集合平均;用最小平方代替最小均方..参见§3.4(最小二乘自适应滤波器)关于最小二乘准则四种方法的讨论.由图可见,当数据长度时(通常都是满足的),计算的求和范围是0.(3)AR系数的求解由式(4.5.4),得图4.5.1自相关法计算预测误差e(n)的原理pppppaaa1,11)1(Nx)1(nx)(pnx)(nx)1(pnx)0(x000000第一个数据移入预测误差滤波器,前向补零最后一个数据移出预测误差滤波器最后一个数据刚进入预测误差滤波器,以后将补零进入.pppppaaa1,11pppppaaa1,11pppppaaa1,11数据长度N点系数长度p+1点1pN)(neppka(4.5.5)式中,——模型参数的维列矢量;——由取样自相关函数(其中)构成的取样自相关矩阵.自相关函数估计为(4.5.6)可见,是有偏自相关估计.将式(4.5.5)对求导并令其等于零(复梯度法),可得aRaxxpjipjxxpipNnpjpjpipipNnNajiRajnxainxaneˆ)(ˆ)()()]([T0,1000102式中的系数N是在定义.时引入的,见式(4.5.6).)(ˆjiRxxT1],,,1[pppaaaAR1pxxRˆ)(ˆmRxxjim)1(,),2(),1(),(ˆ,,2,1,0),()(1)(ˆ10ppmmRpmmnxnxNmRxxmNnxx)(ˆmRxxa(4.5.7)或者写成(4.5.8a)最小预测误差功率(即白噪声方差)为(4.5.8b)式(4.5.8)即为Yule-Walker方程.由此可见,自相关法也是基于解Yule-Walker方程的一种方法.Levinson-durbin递推法是求解Yule-Walker方程的高效算法,具体方法见§3.3.)(ˆ)2(ˆ)1(ˆˆˆˆ)0(ˆ)2(ˆ)1(ˆ)2(ˆ)0(ˆ)1(ˆ)1(ˆ)1(ˆ)0(ˆ21pRRRaaaRpRpRpRRRpRRRxxxxxxppppxxxxxxxxxxxxxxxxxxpmkmRamRpkxxpkxx,,3,2,1,)(ˆˆ)(ˆ1pkxxpkxxwkRaR12min)(ˆˆ)0(ˆˆˆ利用矢量及矩阵求导公式:1minˆpxx0aR由式(4.5.5)求导直接可得:此式即为Yule-Walker方程.(参见清华胡广书p311及p345).T()XAYXAY如果已知信号的个观测数据(,),利用Levinson递推法计算功率谱的流程见西电教材图4.5.1.(4)说明a)由个自相关函数的估值,利用Levinson递推法求解Yule-Walker方程所得的AR模型参数等效于前向预测器的系数,AR模型激励白噪声的方差等效于前向预测的最小预测误差功率.b)自相关法计算相对简单,但需事先根据已知观测数据估计自相关函数(采用有偏估计).在计算预测误差时因作了加窗处理,结果使得分辨率降低.数据越短,分辨率越低,可能还会出现谱峰频率偏移与谱线分裂(即在信号谱峰附近产生虚假谱线).N)(nx10Nn1p2ˆwminˆ4.5.2协方差法与修正协方差法1.协方差法(1)估计准则(4.5.9)注意:本算法仍采用前向预测误差,其求和范围是.(2)预测误差功率计算计算原理如图4.5.2所示.特点:为回避“加窗处理”而引入的频谱卷积效应,即避免在数据段以外补零,本算法“规定”在第一个数据移至位置时开始计算;而当最后一个数据移至时结束计算.计算时仅使用了已获得的观测数据,未在数据两端补零,因此,实际长度为.121)()(1NpnpkpkknxanxpN)(ne1Nnp)(ne)(nx)10(Nn)0(xppa)1(Nx)1(0pa)(ne1N(3)AR系数的求解使用复梯度法使预测误差功率达到最小,可得(4.5.10)图4.5.2协方差法计算预测误差e(n)的原理(三个不同时刻的计算状态:每个参数都对应一个实际的数据值)pppppaaa1,11结束计算pppppaaa1,11pppppaaa1,11)1(px开始计算)2(Nx)(pnx)1(pNx)1(Nx)(nx)1(nx)(pnx)1(pnx)(px)1(x)0(xpka)0,()0,2()0,1(ˆˆˆ),()2,()1,(),2()2,2()1,2(),1()2,1()1,1(21pcccaaappcpcpcpcccpcccxxxxxxppppxxxxxxxxxxxxxxxxxx协方差矩阵式中(4.5.11)称为协方差函数.白噪声的方差为(4.5.12)由式(4.5.10)(4.5.12)即可解出AR模型参数估值和功率谱.(4)说明式(4.5.10)和(4.5.12)构成协方差方程.由于式中不能表示为的函数,所以协方差矩阵不是Toeplitz矩阵,不能采用Levinson递推法求解.b)协方差函数有两个变量,故适合于非平稳随机信号.c)这种方法类似于自相关法,但其分辨率优于自相关法,并可有效地估计正弦信号频率.2.修正协方差法(1)估计准则1)()(1),(NpnxxknxjnxpNkjc1min2),0(ˆ)0,0(ˆˆNpkxxpkxxwkcac212ˆˆˆˆ{,,,;}ppppwaaa注意:模型参数包括:),(kjcxx)(kj使用前向和后向预测误差功率平均值最小准则,即(4.5.13)式中,和是前向和后向预测误差功率,分别表示为(4.5.14)(4.5.15)最小预测误差平均功率是模型输入的白噪声的方差,即.min)(5.0pbpeppepb12112)()(1)(1NpnpkpkNpnfppeknxanxpNnepNpNnpkpkNpnbppbknxanxpNnepN102112)()(1)(12wppkpkknxanx1)()(ˆ(4.5.14)利用前向预测公式:pkpkknxanx1)()(ˆ(4.5.15)利用后向预测公式:关于式(4.5.14)和(4.5.15)的说明:求和范围与协方差法相同,即为后向预测是利用某一时刻以后的个值,即预测的,因此,后向预测值为实际中,总是利用同一组数据来同时实现前向和后向预测的(见图).即利用向后一步预测,因此上式可改写为考虑到前向和后向预测具有相同的预测系数,即.于是预测误差功率令,上式求和范围应改为0,这样就得到式(4.5.15).1pnNnp)}(,),2(),1({pnxnxnxn1n2n1pnpn1pn前向预测后向预测pkpkbknxbnx1)()(ˆ),({nx),1(nx)}1(,pnx)(pnxpkpkbkpnxbpnx1)()(ˆpkpkab121)()(1NpnpkpkpbkpnxapnxpNpnipN1(2)参数的求解求对()的复梯度,并令其等于零,解得(4.5.16)其中:(4.5.17)最小预测误差功率:(4.5.18)由式(4.5.16)(4.5.18)即可解出AR模型参数估值和功率谱.ARpkappkapk,,2,1),0()2,0()1,0(ˆˆˆ),()2,()1,(),2()2,2()1,2(),1()2,1()1,1(21pcccaaappcpcpcpcccpcccxxxxxxppppxxxxxxxxxxxxxxxxxx110)()()()()(21),(NpnpNnxxknxjnxknxjnxpNkjcpkpkxxkcac1min),0(ˆ)0,0(ˆ模型参数包括:212ˆˆˆˆ{,,,;}ppppwaaa(3)说明a)修正协方差法除计算协方差函数不同外,两种计算功率谱的方法是一样的.b)同样不能采用Levinson递推法求解修正协方差的协方差方程.c)适合于估计正弦波频率,其频率位置与正弦相位及噪声的依赖关系小,可得到高分辨率且统计稳定的谱估计.),(kjcxx4.5.3伯格(Burg)法（1）Burg法的基本思路Burg法不直接估计AR参数,而是首先估计反射系数,然后再利用Levinson算法由反射系数求得AR参数.这样就免去了估计自相关函数,可以直接利用信号观测数据求模型参数.（2）估计准则与修正协方差法一样,采用前向和后向预测误差功率平均值最小准则,即(4.5.19)其中,前向,后向预测误差功率分别为(4.5.20)(4.5.21)pkmin)(5.0pbpep12)(1NpnfppenepN12)(1NpnbppbnepN以上二式的求和范围与协方差法相同,即为.式中,前向和后向预测误差分别为(4.5.22)(4.5.23)（3）反射系数的估计将前向和后向预测误差递推公式分别代入式(4.5.20)和(4.5.21),再由(4.5.19)得到1Nnppkpkfpknxanxne1)()()(pkpkbpkpnxapnxne1)()()(pk)()1()()1()()(1111nekn