1、12.5 关于线性预测的进一步讨论上一节使用的AR模型等效于一个 p 阶的线性预测器。即Yule-Walker方程等效于Wiener-Hopf 方程。但估计的功率谱的分辨率不理想，其原因是仅用了前向预测，即1()()pkkx nx nk()()()e nx nx n()x np(1)x n(1)x np()x n对同样一组数据，我们可以实现双向预测：()x np(1)x n(1)x np()x npkffknxkanx1)()()(Forward Prediction)()()(nxnxneff前向预测误差序列2|()|ffEen误差功率1()()()pbbkxnak x nk Backwar
2、d Prediction2|()|bbEe n后向预测误差功率1()()()pbbkxnpak x npk 对同一组数据的后向预测()()()()()bbbbe npx npxnpe npe n后向预测误差序列()()()bbe nx npe np令：0,1,2,()bbkpak可以得到使 最小的 及 。当然也可使用正交原理得：b(1),()bbaapminbmin11(0)()()()()(),1,2,.,pbbxxkpbxxkrak r krmak r mkmp 后向预测的Wiener-Hopef Eqfbminmin)()()()()(*kakakakakaffbfb可以证明：前、后向预
3、测对等关系上述结果表明，使用已知的 p 个数据，我们可以实现前向预测，也可以实现后向预测，两种情况下可各自得到对等的Wiener-Hopf方程。将它们单独使用，所得分辨率都不理想。可以设想，如将二者结合起来，即同时使前向、后向预测误差功率为最小，应能得到更好的分辨率。人们在线性预测方面进行了大量的研究。11*11()()(1)()(1)()ffbmmm mbbfmmm menenk enenenk en前、后向预测误差序列有如下的关系：00()()()fbene nx n1,2,mp初始条件11111111(1),()|(1)|*|(1)|cov(1),()var()var(1)bfmmmfb
4、mmbfmmfbmmenenkenenenenenen 反射系数上述关系引出了线性预测中的Lattice结构。这一结构在现代谱估计、语音信号处理中有着重要的应用。22|()|()|ffnbbnene n上述的关系还是集总平均。对实际的信号：单个样本有限长，求均值要简化，对(),0,1,1x n nN1()()()()()()ffpfkenx nxnx nak x nk22|()|()|ffbbEenEe n取代 的范围n1()()()()()()ffpfkenx nxnx nak x nk(0)0 x(0)(0)fex(1)(1)(0)fxax(1)(1)(1)fexx(1)x N(1)x N
5、p()x p(0)x(1)x p(1)()(1)fx Npap x N (1)(1)()(1)ffeNpx Npap x N (0)(0)(1)(0)(1)(1)(2)(0)(1)()(1)(1)(0)()(1)(2)()(1)(1)(1)(1)()()(1)(2)(1)fffffffxexxex px pxepx px pxxepx Nx Nx Npx NpeNx Nx Npx NpeNx Nx NeNpx 1(1)(2)()(1)fffaaapNN点数据，前向预测误差序列范围3X2X1X(0)(1)(0)(1)(2)(0)()(1)(1)(0)(1)(2)()(1)(1)(1)()(1)(
6、2)(1)xxxx px pxx px pxxx Nx Nx Npx Npx Nx Npx Npx Nx Nx N0:X上三角+中间块+下三角：上、下加窗；(0)()(1)(1)ffffeepeNeNP 0:X1:X()(1)ffepeN 中间块：上、下不加窗；(0)()(1)fffeepeN 2:X中间块+上三角：下不加窗、上加窗；()(1)(1)fffepeNeNp 3:X中间块+下三角：上不加窗、下加窗；12.6 AR模型系数求解算法AR模型系数求解算法很多，人们目前仍在探讨新的求解算法。目前，常用的算法是：1.自相关法 2.Burg算法 3.协方差(covariance)方法；4.改进
7、的协方差算法（modified),又称：Marple 算法 5.最大似然（Maximum Likelihood)估计 3.递推算法：由 求 ，由 递推，还是直接由 递推)(nx)(mrx)(mrx)(nx各算法之间的主要区别：1.的取值范围，即nnef),(10,XX23,XX选择那一个？2.仅用前向预测，还是前后向都预测？即 令 最小，还是 最小？ffbTfpfpfpfpfppNepeeee)1(,),(,),1(),0(一、自相关法(1),(2),()ffffTpaaaapfpfpaXe10fpHfppNnfpfeene|)(|102pffpHoaXXmin001令：使用0X使用前向预测使
8、最小，得注意：矩阵 的结果，即是对有限长数据求出的自相关函数，因此，上式等效于：pffppoaRmin100HXXN自相关法的特点：1.只用前向预测，且 等效前、后加窗,分辨率不好;)(nef2.用 ，得到的 是Toeplits阵，才可能用Levinson算法求解；00XXH1pR3.实际上是我们前面讨论过的Yule-Walker 方 程。方法最简单。12121|()|1|()|Nffppn pNbbppn penNpenNp11*1100()()(1)()(1)()1,2,()()()ffbmmm mbbfmmmmfbenenk enenenk enmpene nx n二、Burg算法使用前
9、、后向预测12fbfb前、后都不加窗Lattice 结构，递推算法1211211*11|)1(|)(|)1()(2NmnbmNmnfmNmnbmfmmnenenenek12*11)|1(1,2,1,)()()()(mmmmmmmmmkmkkmakmakkaka先求：()mmkam令：0fbmk得到 的求解公式：mk 再用Levinson 递推求其它递推步骤1.令：求出2.求 时的参数3.求出 ，再求4.用Levinson算法，求 时的5.重复上述过程，直到);()()(00nxnenebf1k1m1111(1),(1|)(0)xakkr11()()fbene n、2k2m22(1),apm B
10、urg算法：一个公认的较好的算法。Burg 算法的特点：112211|()|,|()|NNffbbppppn pn penenNpNp1.同时使用前向后后向预测，即使12fbfb最小2.的选择保证前、后不加 窗，即(),()fbppen en3.在每一级，仅对 最小，然后套用自相关法的Levinson递推算法，影响分辨率；fbmk4.直接用数据递推，方法简单。三、改进的协方差法Marple方法11221211|()|,|()|fbfbNNffbbppppn pn penenNpNp同Burg算法0()fbmai1,2,1,2,immp注意：这是Marple 算法和Burg算法的最大区别。Bur
11、g算法仅：/0,1,fbmkmp(1,1)(1,2)(1,)(1,0)(1)(2,1)(2,2)(2,)(2,0)(2)()(,1)(,2)(,)(,0)xxxxxxxxxxxxcccpcacccpcaa pcpcpcp pcp 上述最小化的结果是得到一个协方差方程：注意：该矩阵不是Toeplitz矩阵，因此不能用Levinson算法求解。Marple于1983年给出的求解上式的快速递归算法。所以，该算法称作“改进的协方差法，或Marple算法。该算法的估计性能最好，但计算复杂。（e）Burg算法 Burg算法10,()pf13p（g）Marple算法 Marple算法10,()ph13p12
12、.7 MA模型qkknukbnunx1)()()()(qkkzkbzH1)(1)(221()|1()|qjjkxkP eb k e()()H zB z)(nu)(nx(0)1b10()()()()()()()()()xqkqxukx n x nmb k u nmku nm x nb k rr mEEmk0,1,mq220()()()()()0,qq mk mkxb k b kmb k b kmr mmq再推导一步，有：非线性方程组MA模型的正则方程()()qjj mBTxmqPer m e222()()()jjMAqj mxmqPeB er m e从谱估计的角度，MA模型等效于经典法中的间接法
13、，所以分辨率低。因此，MA模型用于谱估计无优势。但，MA模型：1.常用于系统辨识；2.ARMA模型中包含了MA部分。令其等效为 模型求解算法：由于MA模型的正则方程是非线性方程，所以人们提出了很多的求解算法，如谱分解、基于迭代的方法、基于高阶AR模型近似的方法。后者最好用，基础是Wold分解定理。1()()1()qkqkHzHzb k z)MA(q1)(11)(1)(kkzkazAzH对 建立一个无穷阶的AR模型()x n()()1Az B z于是有：()()a kb k步骤：1.由 ,建立 得 ;2.对 建立 阶线性预测器，系数为 ,即建立两次AR模型。1,.1,0),(Nnnxqpp),A
14、R（pkkap1),()(kapqqkkb1),.,(1()()()()qka mb k a mkm1,m 1,mp 12()()()()|()|qppkMAmamb k amke me m近似12.8 ARMA(p,q)模型2101()()()()()()(),pq mxkkxpxka k r mkh k b mkr ma k r mkmq0,1,kqARMA模型的正则方程pm 对第二个式子，()(1)(1)(1)(1)(1)()(2)(2)(2)(1)(2)()()()xxxxxxxxxxxxr qr qr qpr qar qr qr qpr qar qpr qpr qr qpa p ,R
15、arrRae 可以先 求 ，然后再解第一个方程，求出 ；但这样做的效果不好，一是 的性能不好，二是第一个方程也不好求解。首先，建立一个超定方程（方程个数未知数）：paa 11qbb()(1)(1)(1)(1)(2)(1)(2)()()()(1)(2)()()xxxxxxxxxxxxr qr qr qpr qaar qpr qpr qr qpa pr Mr Mr Mpr M paa 11()()()(),1pxxkr ma k r mke mmq 21|()|MHn qe n e e1()HH aR RR r用求伪逆的方法可求出 ；注意，伪逆可用奇异值分解（SVD）的方法求解；求出 后，剩下的工
16、作是求b a apaaa,.,212.用 对 滤波；3.滤波输出 相当于一 MA（q)过程，按 上节MA模型的求解方法，可求 出ARMA(p,q)模型 的 参数。()A z)(nx)(nyARMA 模型系数求解的方法：1.先求出：，它们可构成 ；()A z()()B zA z()u n()y n)(nx()A z()u n()y n()B z（a）MA(10)（b）MA(16)（c）ARMA(10,10)（d）ARMA(10,13)12.10 基于矩阵特征分解的功率谱估计假定信号由 M 个复正弦加白噪声组成：1()exp()MkkkkX nAjnju n 21()exp()()Mxiiwir
17、kAjkk 21()2()MxiiwiPA 已知：(0),(1),()xxxrrr p不会奇异*(0)(1)()(1)(0)(1)()(1)(0)xxxxxxpxxxrrr prrr pRr pr pr(1)(1)pp可构成目标：1.由该矩阵估计 个正弦信号的频率和幅度；2.估计信号 的功率谱；()X nM定义：1,exp(),exp(),1,2,TiiijjpiMe为信号向量，它包含了 个复正弦，其频率和原信号的频率相同。M求解的关键是自相关矩阵的分解：信号相关阵的表示1()exp()()Mxiiwir kAjkk 因为：1MHpiiiwiARe eI所以：1MHpiiiwiARe eI1M
18、HpiiiiSA e epwWIpppRSW相关矩阵的分解：信号部分和噪声部分秩是秩为M秩为1p再定义11pHpiiiiSV V1,0,HijijV Vij特征分解121:,det()0,0,:pMMpNote ifpM thenSandso1MHpiiiiSV V11pHiiiIV V借用特征向量的特点 主特征向量1MVV构成的p+1维空间11pVV构成的M维信号空间1MVV构成的噪声空间11MpVV111()pMHHiwiiwiiii MVVVV信号空间特征值111pMHHpiiiwiiiiRVVVV基于噪声子空间的频率估计和功率谱估计：111MHHpiiiwMMiifMpthenRVVV
19、V噪声空间只有一个特征向量可以证明：1,0,1,2,MiiM Ve 和信号向量正交1M V110()exp()01,2,MHiMMike Vvkj kiM即：110()exp()01,2,MHiMMike Vvkj kiM求解上式，可得到 的 个根，它们都在单位圆上，因此可求出 12,M()V z实现了频率估计M10()()0Mjj kMkV evk e10()()0MkMkV zvk zM 阶多项式方法：1.由 估计 ，由2.构成 ，并假定 ；(),0,1x n nN()xr m(0),()xxrr ppRMp2.对 作特征分解，找最小的 ，及pR1p1pV3.代入上式，解出：实现了频率估计
20、。12,M 4.由下式，求12,MA AA12112212(1)exp()exp()exp()(2)exp(2)exp(2)exp(2)()exp()exp()exp()xMxMxMMrjjjArjjjAr MjMjMjMA求出5.由1(0)MxiwirAw按上述步骤，可求出正弦信号的参数 Pisarenko 谐波分解若噪声空间向量不止一个，估计信号的频率，可应用谱估计的方法。1211()()xpHkkk MPeV 111()()MUSICpHHkkk MPeV Ve1.若1,11kkMpMUSIC(Multiple Signal Classification)方法 111()1()EVpHH
21、kkk MKPeV Ve2.若1/,11kkkMpEV(Eigenvector)方法用特征分解求出的功率谱曲线与本章内容有关的MATLAB文件：1.pyulear.m 用AR模型的自相关法估 计信号的功率谱，其基本调用格式是：Px,F=pyulear(x,order,Nfft,Fs)2.pburg.m 用AR模型的Burg算法估计信 号的功率谱，其基本调用格式是：Px,F=pburg(x,order,Nfft,Fs)（一）有关功率谱估计的MATLAB文件3.pcov.m 用AR模型方差方法估计信号的 功率谱，其基本调用格式是：Px,F=pcov(x,order,Nfft,Fs)4.pmcov.
22、m 用AR模型的改进的方差方法估 计信号的功率谱，其基本调用格式是：Px,F=pmcov(x,order,Nfft,Fs)5.pmem.m 最大熵功率谱估计，其估计 性能类似pyulear,其基本调用格式是：Px,F=pmem(x,order,Nfft,Fs)6.pmusic.m 用自相关矩阵分解的MUSIC 算法估计信号的功率谱，其基本调用格 式是：Px,F=pmusic(x,order,Nfft,Fs)7.peig.m 用自相关矩阵分解的特征向量 法估计信号的功率谱，其基本调用格式是：Px,F=peig(x,order,Nfft,Fs),Px,F，V,E=peig(x,order,Nfft
23、,Fs),x:信号向量，order:模型的阶次，Fs:抽样频率，Nfft:对x作FFT时的长度。Px:估计出的功率谱，F是频率轴坐标。对peig,输出的E 是由自相关矩阵的特征值所组成的向量，V是由特征向量组成的矩阵。V的列向量张成了噪声子空间，V的行数减去列数即是信号子空间的维数。（二）有关（二）有关AR模型参数估计的文件：模型参数估计的文件：包括：aryule,arburg,arcov 及 armcov。8.aryule.m 用自相关法(即Yule-Walker法)估 计AR模型的参数，其基本调用格式是：a,E=aryule(x,order),a,E，k=aryule(x,order)9.
24、arburg.m 用Burg算法估计AR模型的参数，其基本调用格式是：a,E=arburg(x,order)a,E，k=arburg(x,order)10.arcov.m 用方差方法估计AR模型的参数，其基本调用格式是：a,E=arcov(x,order)11.armcov.m 用改进的方差方法估计AR模型 的参数，其基本调用格式是：a,E=armcov(x,order)x：信号向量；order：模型的阶次；a：AR模型系数向量；E：AR模型输入白噪声的功率，或order阶线 性预测器的最小预测误差。k：反射系数向量。（三）有关线性预测的MATLAB文件12.lpc：用来计算线性预测系数。13
25、 a=lpc(x,order);其作用等同于 aryule；13.ac2poly：由自相关函数求线性预测系数 14 a,E=ac2poly(R);14.Poly2ac：由线性预测系数求自相关函数。Rpoly2ac(a,E);15.ac2rc 由自相关函数得到反射系数及。16 k,R0=ac2rc(R);16.rc2ac 由反射系数及得到自相关函数。Rrc2ac(k,R0);17.poly2rc 由线性预测系数得到反射系数 k=poly2rc(a),或 k,R0=poly2rc(a,E);17.rc2poly 由反射系数及得到线性预测系数。a=rc2poly(k)，或 a,E=rc2poly(k,R0);18.Levinson：用Levinson-Durbin 算法求解 Toeplitz 矩阵，该文件是一个C-MEX 内部 文件。以上多个文件都要调用它。a,E,k=levinson(R,order)。
