清華大學(xué)數(shù)字信號處理課件第12章

1、12.5 關(guān)于線性預(yù)測的進一步討論上一節(jié)使用的ar模型等效于一個 p 階的線性預(yù)測器。即yule-walker方程等效于wiener-hopf 方程。但估計的功率譜的分辨率不理想，其原因是僅用了前向預(yù)測，即1( )()pkkx nx nk ( )( )( )e nx nx n()x np(1)x n(1)x np( )x n對同樣一組數(shù)據(jù)，我們可以實現(xiàn)雙向預(yù)測：()x np(1)x n(1)x np( )x npkffknxkanx1)()()(forward prediction)()()(nxnxneff前向預(yù)測誤差序列2|( )| ffeen誤差功率1 ( )( ) ()pbbkxnak

2、 x nk backward prediction2|( )| bbee n后向預(yù)測誤差功率1 ()( ) ()pbbkxnpak x npk 對同一組數(shù)據(jù)的后向預(yù)測()()()()( )bbbbe npx npxnpe npe n后向預(yù)測誤差序列( )()()bbe nx npe np令：0,1,2,( )bbkpak可以得到使 最小的 及 。當(dāng)然也可使用正交原理得：b(1),( )bbaapminbmin11(0)( ) ( )( )( ) (),1,2,.,pbbxxkpbxxkrak r krmak r mkmp 后向預(yù)測的wiener-hopef eqfbminmin)()()()(

3、)(*kakakakakaffbfb可以證明：前、后向預(yù)測對等關(guān)系上述結(jié)果表明，使用已知的 p 個數(shù)據(jù)，我們可以實現(xiàn)前向預(yù)測，也可以實現(xiàn)后向預(yù)測，兩種情況下可各自得到對等的wiener-hopf方程。將它們單獨使用，所得分辨率都不理想?？梢栽O(shè)想，如將二者結(jié)合起來，即同時使前向、后向預(yù)測誤差功率為最小，應(yīng)能得到更好的分辨率。人們在線性預(yù)測方面進行了大量的研究。11*11( )( )(1)( )(1)( )ffbmmm mbbfmmm menenk enenenk en前、后向預(yù)測誤差序列有如下的關(guān)系：00( )( )( )fbene nx n1,2,mp初始條件11111111(1),( )|(

4、1)|*|(1)|cov(1),( )var( ) var(1)bfmmmfbmmbfmmfbmmenenkenenenenenen 反射系數(shù)上述關(guān)系引出了線性預(yù)測中的lattice結(jié)構(gòu)。這一結(jié)構(gòu)在現(xiàn)代譜估計、語音信號處理中有著重要的應(yīng)用。22|( )|( )|ffnbbnene n上述的關(guān)系還是集總平均。對實際的信號：單個樣本有限長，求均值要簡化，對( ),0,1,1x n nn1( )( )( )( )( ) ()ffpfkenx nxnx nak x nk22|( )| |( )| ffbbeenee n取代 的范圍n1( )( )( )( )( ) ()ffpfkenx nxnx na

5、k x nk (0)0 x(0)(0)fex(1)(1) (0)fxax (1)(1)(1)fexx(1)x n (1)x np ( )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)fffaaap

6、nn點數(shù)據(jù)，前向預(yù)測誤差序列范圍3x2x1x(0)(1)(0)(1)(2)(0)( )(1)(1)(0)(1)(2)()(1)(1)(1) ()(1)(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模型系數(shù)求解

7、算法ar模型系數(shù)求解算法很多，人們目前仍在探討新的求解算法。目前，常用的算法是： 1. 自相關(guān)法 2. burg算法 3. 協(xié)方差(covariance)方法； 4. 改進的協(xié)方差算法（modified ) , 又稱：marple 算法 5. 最大似然（maximum likelihood)估計 3. 遞推算法：由 求 ，由 遞推，還是直接由 遞推)(nx)( mrx)( mrx)(nx各算法之間的主要區(qū)別：1. 的取值范圍，即nnef),(10, xx23,xx選擇那一個？ 2. 僅用前向預(yù)測，還是前后向都預(yù)測？即 令 最小，還是 最??？ffbtfpfpfpfpfppnepeeee)1(,)

8、,(,),1 (),0(一、自相關(guān)法(1),(2),( )fffftpaaaapfpfpaxe10fphfppnnfpfeene| )(|102pffphoaxxmin001令：使用0x使用前向預(yù)測使最小，得注意：矩陣 的結(jié)果，即是對有限長數(shù)據(jù)求出的自相關(guān)函數(shù)，因此，上式等效于：pffppoarmin100hxxn自相關(guān)法的特點：1. 只用前向預(yù)測，且 等效前、后加窗, 分辨率不好;)(nef2. 用 ，得到的 是toeplits陣，才可能用levinson算法求解；00xxh1pr3. 實際上是我們前面討論過的yule-walker 方 程。方法最簡單。12121|( )|1|( )|nff

9、ppn pnbbppn pennpennp11*1100( )( )(1)( )(1)( )1,2,( )( )( )ffbmmm mbbfmmmmfbenenk enenenk enmpene nx n二、burg算法使用前、后向預(yù)測12fbfb前、后都不加窗l(fā)attice 結(jié)構(gòu)，遞推算法1211211*11| ) 1(| )(|) 1()(2nmnbmnmnfmnmnbmfmmnenenenek12*11)|1 (1, 2 , 1,)()()()(mmmmmmmmmkmkkmakmakkaka先求： ( )mmkam令：0fbmk得到 的求解公式：mk 再用levinson 遞推求其它遞推

10、步驟1、 令： 求出2、求 時的參數(shù)3、求出 ，再求4、用levinson算法，求 時的5、重復(fù)上述過程，直到);()()(00nxnenebf1k1m1111(1),(1 |) (0)xakkr11( )( )fbene n、2k2m22 (1),apm burg算法：一個公認(rèn)的較好的算法。burg 算法的特點：112211|( )| ,|( )|nnffbbppppn pn penennpnp1. 同時使用前向后后向預(yù)測，即使12fbfb最小2. 的選擇保證前、后不加 窗，即( ),( )fbppen en3. 在每一級， 僅對 最小，然后套用自相關(guān)法的levinson遞推算法，影響分辨率

11、；fbmk4. 直接用數(shù)據(jù)遞推，方法簡單。三、改進的協(xié)方差法marple方法11221211|( )| ,|( )|fbfbnnffbbppppn pn penennpnp同burg算法0 ( )fbmai1,2,1,2,immp注意：這是marple 算法和burg算法的最大區(qū)別。burg算法僅：/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 上述最小化的結(jié)果是得到一個協(xié)方差方程：注意：該矩陣不是toepl

12、itz矩陣，因此不能用levinson算法求解。marple于1983年給出的求解上式的快速遞歸算法。所以，該算法稱作“改進的協(xié)方差法，或marple算法。該算法的估計性能最好，但計算復(fù)雜。（e）burg算法 burg算法10, ( )pf13p（g）marple算法 marple算法10, ( )ph13p12.7 ma模型qkknukbnunx1)()()()(qkkzkbzh1)(1)(221()|1( )|qjjkxkp eb k e( )( )h zb z)(nu)(nx(0)1b10( )( ) ()( ) ()() ( )( )()xqkqxukx n x nmb k u nmk

13、u nm x nb k rr meemk0,1,mq220( ) ()( ) ()( )0,qq mk mkxb k b kmb k b kmr mmq再推導(dǎo)一步，有：非線性方程組ma模型的正則方程()( )qjj mbtxmqper m e222()()( )jjmaqj mxmqpeb er m e從譜估計的角度，ma模型等效于經(jīng)典法中的間接法，所以分辨率低。因此，ma模型用于譜估計無優(yōu)勢。但，ma模型： 1. 常用于系統(tǒng)辨識； 2. arma模型中包含了ma部分。令其等效為 模型求解算法：由于ma模型的正則方程是非線性方程，所以人們提出了很多的求解算法，如譜分解、基于迭代的方法、基于高階

14、ar模型近似的方法。后者最好用，基礎(chǔ)是wold分解定理。1( )( )1( )qkqkhzhzb k z )ma(q1)(11)(1)(kkzkazazh對 建立一個無窮階的ar模型( )x n( ) ( )1az b z于是有：( )( )a kb k步驟：1、由 ,建立 得 ;2、對 建立 階線性預(yù)測器，系數(shù)為 ,即建立兩次ar模型。1,.1 , 0),(nnnxqpp),ar（pkkap1),()(kapqqkkb1),.,(1( )( ) ()( )qka mb k a mkm1,m 1,mp 12( )( )()( )| ( )|qppkmamamb k amke me m近似12.

15、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 ,rarrrae 可以先 求 ，然后再解第一個方程，求出 ；但這樣做的效果不好，一是 的性能不好，二是第一個方程也不好求解。首先，建立一個超定方程（方程個數(shù)未知數(shù)）： paa 11q

16、bb( )(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）的方法求解；求出 后，剩下的工作是求b a apaaa,.,212、用 對 濾波；3、 濾波輸出 相當(dāng)于一 ma（q) 過程，按 上節(jié)ma模型的求解方法，可求 出arma(p,q)模型 的

17、 參數(shù)。( )a z)(nx)(nyarma 模型系數(shù)求解的方法：1 先求出： ，它們可構(gòu)成 ；( )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 個復(fù)正弦加白噪聲組成：1( )exp( )mkkkkx najnju n 21( )exp()( )mxiiwir kajkk 21( )2()mxiiwipa 已知：(0),(1),( )xxxrrr p不會奇異*(0)(1)( )

18、(1)(0)(1)( )(1)(0)xxxxxxpxxxrrr prrr prr pr pr(1) (1)pp可構(gòu)成目標(biāo)：1. 由該矩陣估計 個正弦信號的頻率和幅度； 2. 估計信號 的功率譜；( )x nm定義：1,exp(),exp() ,1,2,tiiijjpime為信號向量，它包含了 個復(fù)正弦，其頻率和原信號的頻率相同。m求解的關(guān)鍵是自相關(guān)矩陣的分解：信號相關(guān)陣的表示1( )exp()( )mxiiwir kajkk 因為：1mhpiiiwiare ei所以：1mhpiiiwiare ei1mhpiiiisa e epwwippprsw相關(guān)矩陣的分解：信號部分和噪聲部分秩是秩為m秩為1

19、p再定義11phpiiiisv v1,0,hijijv vij特征分解121:,det()0,0,:pmmpnote ifpm thensandso1mhpiiiisv v11phiiiiv v借用特征向量的特點 主特征向量1mvv構(gòu)成的p+1維空間11pvv構(gòu)成的m維信號空間1mvv構(gòu)成的噪聲空間11mpvv111()pmhhiwiiwiiii mvvvv信號空間特征值111pmhhpiiiwiiiirvvvv基于噪聲子空間的頻率估計和功率譜估計：111mhhpiiiwmmiifmpthenrvvvv噪聲空間只有一個特征向量可以證明：1,0,1,2,miim ve 和信號向量正交1m v11

20、0( )exp()01,2,mhimmike vvkj kim即：110( )exp()01,2,mhimmike vvkj kim求解上式，可得到 的 個根，它們都在單位圓上，因此可求出 12,m ( )v z實現(xiàn)了頻率估計m10()( )0mjj kmkv evk e10( )( )0mkmkv zvk zm 階多項式方法：1. 由 估計 ，由 構(gòu)成 ，并假定 ；( ),0,1x n nn ( )xr m(0),( )xxrr pprmp2. 對 作特征分解，找最小的 ，及pr1p1pv3. 代入上式，解出： 實現(xiàn)了頻率估計。 12,m 4. 由下式，求12,ma aa12112212(1

21、)exp()exp()exp()(2)exp( 2)exp( 2)exp( 2)()exp()exp()exp()xmxmxmmrjjjarjjjar mjmjmjma求出5. 由1(0)mxiwiraw按上述步驟，可求出正弦信號的參數(shù) pisarenko 諧波分解若噪聲空間向量不止一個，估計信號的頻率，可應(yīng)用譜估計的方法。1211( )( )xphkkk mpev 111( )( )musicphhkkk mpev ve1. 若1,11kkmpmusic(multiple signal classification)方法 111( )1( )evphhkkk mkpev ve2. 若1/,1

22、1kkkmpev(eigenvector)方法用特征分解求出的功率譜曲線與本章內(nèi)容有關(guān)的matlab文件： 1. pyulear.m 用ar模型的自相關(guān)法估 計信號的功率譜，其基本調(diào)用格式是： px, f = pyulear(x, order, nfft, fs) 2. pburg.m 用ar模型的burg算法估計信 號的功率譜，其基本調(diào)用格式是： px, f = pburg(x, order, nfft, fs)（一）、 有關(guān)功率譜估計的matlab文件3. pcov.m 用ar模型方差方法估計信號的 功率譜，其基本調(diào)用格式是： px, f = pcov(x, order, nfft, fs

23、)4. pmcov.m 用ar模型的改進的方差方法估 計信號的功率譜，其基本調(diào)用格式是：px, f = pmcov(x, order, nfft, fs) 5. pmem.m 最大熵功率譜估計，其估計 性能類似pyulear, 其基本調(diào)用格式是： px, f = pmem(x, order, nfft, fs)6. pmusic.m 用自相關(guān)矩陣分解的music 算法估計信號的功率譜，其基本調(diào)用格 式是： px, f = pmusic(x, order, nfft, fs)7. peig.m 用自相關(guān)矩陣分解的特征向量 法估計信號的功率譜，其基本調(diào)用格式是：px, f = peig(x, or

24、der, nfft, fs), px, f，v, e = peig(x, order, nfft, fs),x :信號向量，order:模型的階次，fs:抽樣頻率，nfft:對x作fft時的長度。px:估計出的功率譜，f是頻率軸坐標(biāo)。對peig, 輸出的e 是由自相關(guān)矩陣的特征值所組成的向量，v是由特征向量組成的矩陣。v的列向量張成了噪聲子空間，v的行數(shù)減去列數(shù)即是信號子空間的維數(shù)。 （二）有關(guān)（二）有關(guān)ar模型參數(shù)估計的文件：模型參數(shù)估計的文件：包括：aryule, arburg, arcov 及 armcov。8. aryule.m 用自相關(guān)法(即yule-walker法)估 計ar模型的參數(shù)，其基本調(diào)用格式是： a, e = ar