Chapter+4中國科學(xué)院大學(xué)現(xiàn)代數(shù)字信號(hào)處理課程課件

1、4.1 引言 4.2 經(jīng)典譜估計(jì) 4.3 現(xiàn)代譜估計(jì)中的參數(shù)建模 4.4 AR模型譜估計(jì)方法 4.5 最大熵譜估計(jì)方法4.6 特征分析法譜估計(jì)功率譜定義功率譜估計(jì)的方法功率譜估計(jì)的應(yīng)用信號(hào)的功率譜和其自相關(guān)函數(shù)服從一對(duì)傅里葉變換關(guān)系mjmxxjxxmrPe )()e (dPmrnjjxxxxe )e (21)(經(jīng)典譜估計(jì)方法間接方法：BT法直接方法：周期圖法現(xiàn)代譜估計(jì)方法參數(shù)法：ARMA模型法（AR模型、MA模型、ARMA模型）非參數(shù)法：Pisarenko諧波分解法、MUSIC方法在信號(hào)處理的許多場所，要求預(yù)先知道信號(hào)的功率譜密度(或自相關(guān)函數(shù))；常常利用功率譜估計(jì)來得到線性系統(tǒng)的參數(shù)估計(jì)；從

2、寬帶噪聲中檢測窄帶信號(hào)。 22)()(HPyyBT法周期圖法改進(jìn)的周期圖法BT法是先估計(jì)自相關(guān)函數(shù), 然后進(jìn)行傅里葉變換得到功率譜。有偏自相關(guān)函數(shù)估計(jì)的誤差相對(duì)較小，是一種漸近一致估計(jì)：1|0*)()(1)(mNnxxmnxnxNmr-BT(e)( )ejjmxxmPrmn 周期圖法的定義如下: 212j-j011(e )( )eNnjxxnPx nX eNN觀 測 數(shù) 據(jù) x(n)FFT取 模的 平 方1/NPxx(ej )1. 1. 周期圖與周期圖與BTBT法的等價(jià)關(guān)系法的等價(jià)關(guān)系 2101111*()00001(e)( )e11( )e( )e( )( )eNjj nxxnNNNNj k

3、j njk nknknPx nNx kx nx k x nNN令 m=k-n, 即k=m+n，則 1*01j(1)1(1)BT(e )ee(e1( ) () ( )Nj mxxmNmnxNNj mmxNjPPx n x mnNrm 利用有偏自相關(guān)函數(shù)的BT法和周期圖法是等價(jià)的。 2. 周期圖法譜估計(jì)質(zhì)量分析周期圖法譜估計(jì)質(zhì)量分析 1) 周期圖的偏移 1j-j(1)1-j(1)-j ( )|(e)ee( )( )e)NmxxmNNmmNmxxBxmxxxE rmNmrE Pwm rmmN其它0|)(NmNmNmwB式中 上式在頻域表示為： j-1(e)(e)ed2jjxxBxxE PWP式中 )

4、()e (mrFTPxxjxx2)2/sin()2/sin(1)()e (NNmwFTWBjB 周期圖的統(tǒng)計(jì)平均值等于它的真值卷積三角譜窗函數(shù)，因此周期圖是有偏估計(jì)，但當(dāng)N時(shí)，wB(m)1，三角譜窗函數(shù)趨近于函數(shù)，周期圖的統(tǒng)計(jì)平均值趨于它的真值，因此周期圖屬于漸近無偏估計(jì)。 2) 周期圖的方差 用這種方法估計(jì)的功率譜在2x附近起伏很大，故周期圖是非一致估計(jì)。 白噪聲的周期圖 Bartlett平均周期圖法窗口處理法平均周期圖Welch法（修正的周期圖求平均法） 主要思想：主要思想：對(duì)序列x(n)進(jìn)行L次獨(dú)立觀測或?qū)⑵浞殖蒐段，計(jì)算每組觀測數(shù)據(jù)的周期圖，再將L個(gè)周期圖加和后求平均。 LiixxIL

5、P1j)(1)e ( 假設(shè)隨機(jī)信號(hào)x(n)的觀測數(shù)據(jù)區(qū)間為：0nN-1，共進(jìn)行了L次獨(dú)立觀測，得到L組記錄數(shù)據(jù)，每一組記錄數(shù)據(jù)用xi(n), i=1, 2, 3, ,L表示； 或?qū)﹂L為N的數(shù)據(jù)x(n)分成L段，每段有M個(gè)數(shù)據(jù)，N=LM，第i段數(shù)據(jù)表示為xi(n)= x(n+iM-M)。 第i組的周期圖用下式表示： 210)(1)(MnnjiienxMIn 估計(jì)方法：估計(jì)方法： 將得到的L個(gè)周期圖進(jìn)行平均，作為信號(hào)x(n)的功率譜估計(jì)， 公式如下： LiixxILP1j)(1)e (n 估計(jì)效果分析：估計(jì)效果分析： 平均周期圖的估計(jì)方差是周期圖的方差的1/L，L越大方差越小，功率譜越平滑；相應(yīng)的

6、，M越小，偏移越大，分辨率越低；估計(jì)的均方誤差也減少； 以分辨率的降低換取了估計(jì)方差的減少，估計(jì)量的方差和分辨率是一對(duì)矛盾。平均周期圖法 Pxx(ej)/032101230321012303210123Pxx(ej)/Pxx(ej)/(a)(b)(c)N256， L2N256， L4N256， L8主要思想：主要思想：用一適當(dāng)?shù)墓β首V窗函數(shù)W(ej)與周期圖進(jìn)行卷積，來達(dá)到使周期圖平滑的目的。 (e)( )(e)jjxxNPIW()0 jW e其中，()1(e)( )(e)( )(e)d2jjjxxNNPIWIW 式中 mxxNNmNmrIj -1)1(e )()(de)e(21)(jjnWn

7、w-(M-1)nM-1 n 估計(jì)方法：估計(jì)方法：那么 mxxMMmjxxemwmrPj1)1()()()e (mxxMMmxxmwmrEPEj -1)1(je )()( )e (又 1j-j(1)()( )( ) ( )eMmxxxxBmME Perm wm w m 偏移分析： n 估計(jì)效果分析：估計(jì)效果分析：| ( )( )( )( )xxxxxxBNmE rmrmrmwmN可得 周期圖的窗函數(shù)法仍然是有偏估計(jì), 其偏移和wB(m)、w(m)兩個(gè)窗函數(shù)有關(guān)。 如果w(m)窗的寬度比較窄，M比N小得多，這樣|m|p AR譜估計(jì)方法可歸結(jié)為求解AR模型系數(shù)或線性預(yù)測器系數(shù)的問題。AR模型參數(shù)估計(jì)

8、方法：信號(hào)預(yù)測誤差最小原則（或預(yù)測誤差功率最?。┳韵嚓P(guān)法（Levison遞推法）Burg法協(xié)方差法修正協(xié)方差法（前后向線性預(yù)測最小二乘法）最大熵原則最大熵譜估計(jì)方法1、 自相關(guān)法自相關(guān)法列文森（列文森（Levinson）遞推）遞推 估計(jì)估計(jì)方法：方法：自相關(guān)法的出發(fā)點(diǎn)是選擇AR模型的參數(shù)使預(yù)測誤差功率最?。徊捎肔evison-Durbin遞推方法求解Yule-Walker方程得到AR模型參數(shù)。 21j -22j2je11| )e (|)e (PiiiwwxxaHP 預(yù)測誤差功率為 22111| ( )|( )()ppinnie nx na x niNN假設(shè)信號(hào)x(n)的數(shù)據(jù)區(qū)間在0nN-1范圍

9、，有P個(gè)預(yù)測系數(shù)，N個(gè)數(shù)據(jù)經(jīng)過沖激響應(yīng)為api(i=0,1, 2, , P)的濾波器， 輸出預(yù)測誤差e(n)的長度為N+P， 因此應(yīng)用下式計(jì)算: 2110210| )()(|1| )(|1pipiPNnPNninxanxNneN0( )( )( )()pepiie nx nh na x ni預(yù)測誤差功率最小，得到 )()2() 1 ()0()2() 1()2()0() 1 () 1() 1()0(21prrraaarprprprrrprrrxxxxxxppppxxxxxxxxxxxxxxxxxx0pkpkaa 采用Levinson-Durbin遞推法求解Yule-Walker方程： 由k=1開

10、始遞推，遞推到k=p，依次得到a11,21，a21,a22, 22，ap1,ap2,app,2p。 AR模型的各個(gè)系數(shù)以及模型輸入白噪聲方差求出后， 信號(hào)功率譜用下式計(jì)算： 21j -22j2je11| )e (|)e (PiiiwwxxaHP利用列文森遞推法計(jì)算功率譜的流程圖 輸入： x(n),n0,1,2,3,N1,2計(jì)算：)0(1),0(/ ) 1 (21 , 1211 , 1xxxxxxrarra)(mrxxp=22122*, 1, 12111, 111, 2 , 1 , 0)()(ppppkppppkppkppixxipxxppapkaaaaiprapra22ppiipiwxxaP1

11、j2je11)(e結(jié)束p=p+1NY性能分析：性能分析：該方法需要基于有限的觀測數(shù)據(jù)估計(jì)自相關(guān)序列，當(dāng)數(shù)據(jù)長度較短時(shí)，估計(jì)誤差會(huì)比較大，AR參數(shù)的計(jì)算就會(huì)引入很大的誤差。從而導(dǎo)致功率譜估計(jì)出現(xiàn)譜線分裂與譜峰頻率偏移等現(xiàn)象。2、 伯格（伯格（Burg）遞推法）遞推法 n 估計(jì)方法：估計(jì)方法：直接由時(shí)間序列計(jì)算AR模型參數(shù)的方法，求前、后向預(yù)測誤差平均功率最小時(shí)的反射系數(shù)kp，進(jìn)而求AR模型參數(shù)ak和2w。20,ppkwpkak 設(shè)信號(hào)x(n)觀測數(shù)據(jù)區(qū)間為：0nN-1，前向、后向預(yù)測誤差功率分別用p,e和p,b表示，預(yù)測誤差平均功率用p為 1122,11|( )|( )|1()2NNfbp ep

12、p bpn pn ppp ep benenNpNppkpkbppkpkfpkpnxapnxneknxanxne1*1)()()()()()(其中，前向、后向預(yù)測誤差公式分別為 求預(yù)測誤差平均功率p最小時(shí)的反射系數(shù)kp，令 0ppk1*11122112( )(1)(|( )|(1)| )Nfbppn ppNfbppn pen enkenen 基于反射系數(shù)kp，由Levinson-Durbin遞推關(guān)系求AR模型參數(shù)ak和2w，進(jìn)而求得功率譜Pxx 伯格遞推法流程圖 輸入： x(n),n0,1,2,3,N1階數(shù)IP1)(1)()()()(210000pnxNnxnenxneNnbfp = IP結(jié)束1

13、 pp12121*111)()() 1()(2NpnbpfpbpNpnfppnenenenekpwpppk2121pppipppippikaakaa*, 1, 11, 2, 1),() 1()(1, 2, 1),1()()(1*111NppnnekneneNppnneknenefppbpbpbpppfp2,;, 3 , 2 , 1:wippia輸出21j2je1)(epiipiwxxaPNY性能分析：性能分析：該方法避免了采用有限數(shù)據(jù)估計(jì)自相關(guān)函數(shù)的計(jì)算，適合短序列參數(shù)估計(jì)，克服了L-D遞推中的某些缺點(diǎn)，計(jì)算量小。但對(duì)正弦信號(hào)的譜估計(jì)，仍存在某些譜線分裂與頻率偏移現(xiàn)象。3 3、 協(xié)方差法與修正

14、協(xié)方差法協(xié)方差法與修正協(xié)方差法（1 1）. . 協(xié)方差法協(xié)方差法n 估計(jì)方法：估計(jì)方法：利用使預(yù)測誤差功率最小的方法求模型參數(shù)11212| )()(|1| )(|1NpnNpnpkpkknxanxpNnepN 該公式中使用的觀測數(shù)據(jù)均已得到，不需要在數(shù)據(jù)兩端補(bǔ)充零點(diǎn)， 因此比較自相關(guān)法去掉了加窗處理的不合理假設(shè)。0pkpkaan 性能分析：性能分析： 一些實(shí)驗(yàn)結(jié)果說明它的分辨率優(yōu)于自相關(guān)法，另外對(duì)于純正弦信號(hào)數(shù)據(jù)，可以有效地估計(jì)正弦信號(hào)的頻率。 （2 2）. . 修正協(xié)方差法（前后向線性預(yù)測最小二乘法）：修正協(xié)方差法（前后向線性預(yù)測最小二乘法）：n 估計(jì)方法：估計(jì)方法：修正協(xié)方差法使用前向和后

15、向預(yù)測誤差平均值最小的方法， 估計(jì)AR模型的參數(shù)，進(jìn)而估計(jì)信號(hào)的功率譜。,0pp kp kaa,2p fp bp 前向和后向預(yù)測誤差功率p,f、p,b分別用下式表示： 21,121*,011( )()1( )()pNp fpkn pkNppp bpknkx na x nkNpx na x nkNp 預(yù)測誤差平均功率最小1*,1,1*,011( )()()( )()()01,2,3,pNpp kn pkp lNppp knkx nax nkx nlaNpx nax nkx nllp n 性能分析：性能分析： 該方法去掉了Burg法所用的Levinson的約束條件，估計(jì)得到的譜在譜線分裂和頻率偏移

16、時(shí)較Burg法有較大改善。 自相關(guān)法可以用Levinson遞推算法，運(yùn)算量小，但分辨率受窗長度的限制；協(xié)方差法，去除了自相關(guān)法加窗處理的不合理假設(shè)，分辨率高，運(yùn)算量較大；修正協(xié)方差法，分辨率高，在譜線分裂和偏移上較Burg法有較大改善，運(yùn)算量大；Burg算法，可用改進(jìn)的Levinson遞推算法，分辨率高，但對(duì)正弦信號(hào)存在譜線分裂和偏移現(xiàn)象。 例例4.5.1 已知信號(hào)的四個(gè)觀察數(shù)據(jù)為x(n)=x(0), x(1), x(2), x(3)=2, 4, 1, 3， 分別用自相關(guān)法和協(xié)方差法估計(jì)AR（1）模型參數(shù)。 解解（1） 自相關(guān)法： 5 . 009)3()41 (4)24(20) 1()(21)

17、()(21) 1()()(| )(|4111111111114011401111402aaaaanxneaneneanxanxnenennn (2) 協(xié)方差法: 714. 00)3()41 (4)24(20) 1()(32)()(32) 1()()(| )(|31| )(|111111111311131111131212aaaanxneaneneanxanxnenenepNnnnNpn4 4、 關(guān)于關(guān)于ARAR模型階次的選擇模型階次的選擇n對(duì)于白噪中的AR信號(hào)，其階次的選擇應(yīng)折衷考慮。如選擇AR模型，其階次應(yīng)加大，較低的階次會(huì)使譜估計(jì)產(chǎn)生偏移， 降低分辨率。信噪比愈低，平滑作用愈嚴(yán)重，愈需要高

18、的階次， 因此信噪比低應(yīng)選高的階次。階次愈高，分辨率愈高；但階次太高，會(huì)使估計(jì)誤差加大，譜峰分裂。AR模型階次太小時(shí)的平滑作用50.0030.0010.0010.0030.000.0000.400.50真實(shí)AR(4)的PSDAR(2)的PSD模型PSD/dB/n 最終預(yù)測誤差（最終預(yù)測誤差（FPEFPE）準(zhǔn)則：）準(zhǔn)則： 用下式計(jì)算后再取最小，估計(jì)模型階次： kmin( )kNkFPE kNk式中, 表示k階AR模型的白噪聲方差（預(yù)測誤差功率）的估計(jì)值，N是觀測數(shù)據(jù)的長度，k是模型階次。上式中, 隨k增加而減少，但(N+k)/(N-k)卻隨k增加而加大。kkn 阿凱克信息論

19、準(zhǔn)則：阿凱克信息論準(zhǔn)則： kminAIC( )ln2kkNk 選擇使上式最小的k值作為模型的階次。但這種方法有過高估計(jì)模型階次的趨勢。為改進(jìn), 提出用k lnN代替公式中的2k的方法。 AIC和FPE估計(jì)階次的性能是相近的， 但對(duì)于短記錄數(shù)據(jù)，建議使用AIC準(zhǔn)則，對(duì)于長記錄（N）, 兩種方法得到相同的模型階次。 n 自回歸傳遞函數(shù)準(zhǔn)則（自回歸傳遞函數(shù)準(zhǔn)則（CATCAT）：）： k1111minCAT( )kiikkN式中 iiiNN同樣, 模型的階次取使上式最小的k值。最大熵譜估計(jì)與AR模型譜估計(jì)是一致的，采用最大熵原則外推自相關(guān)函數(shù)，得到自相關(guān)函數(shù)與模型參數(shù)服從Yule-Walker方程；通

20、過解Yule-Walker方程，解出模型參數(shù)，進(jìn)而計(jì)算信號(hào)功率譜。1 1 利用最大熵的原則外推自相關(guān)函數(shù)利用最大熵的原則外推自相關(guān)函數(shù) 按照Shannon對(duì)熵的定義， 當(dāng)隨機(jī)變量X取離散值時(shí)，熵的定義為 iiippHln式中pi是出現(xiàn)狀態(tài)i的概率。當(dāng)X取連續(xù)值時(shí)，熵的定義為 )(lnd)(ln)(xpExxpxpH式中, p(x)是X的概率密度函數(shù)， 假設(shè)x(n)是零均值正態(tài)分布的平穩(wěn)隨機(jī)序列，它的N維高斯概率密度函數(shù)為 XNRXNRxxxpxxHxxNN12/12/21)(21exp)(det)2(),(式中 H21,NxxxX)0() 1()() 1()0() 1 ()() 1()0()(

21、xxxxxxxxxxxxxxxxxxxxrNrNrNrrrNrrrNR先討論一維高斯分布的信號(hào)的熵，然后推廣到N維。 222222222( )ln( )( )ln22( )ln2( )21ln22ln2lnln2xHp xp xdxp xdxx p x dxp x dxee 同理可求得N維高斯分布信號(hào)的熵為/21/2ln(2e)(det()NxxHRN式中det(Rxx(N)表示矩陣Rxx(N)的行列式，由上式表明為使熵最大，要求det(Rxx(N)最大。 maxmaxdet()xxHRN由觀測序列x(n)估計(jì)自相關(guān)函數(shù) ；用 代替Levison-Durbin算法里的rxx(m)，求解AR模型

22、參數(shù) 和 ；計(jì)算功率譜 ；(e)jxxP xxrmka2w21(e)1ejwxxNj kkkPa xxrm該方法特別適用于對(duì)多正弦加白噪聲序列進(jìn)行譜分析，可以得到比AR模型法更高的分辨率和更準(zhǔn)確的頻率估計(jì)，尤其在信噪比低時(shí)更有效。估計(jì)思想：將白噪聲加正弦波作為一特殊的ARMA模型，用特征方程求它的參數(shù)，計(jì)算出正弦波的頻率和功率，及噪聲功率。1 1 正弦波用退化正弦波用退化ARAR模型模型表示表示令隨機(jī)相位實(shí)正弦信號(hào)x(n)=sin(n+), 由下面三角恒等式：) 1(sincos2)2(sin)sin(nnn- )2() 1(cos2)(nxnxnx可得：x(n)=-a1x(n-1)-a2x(

23、n-2)122cos,1aa 其中，正弦波信號(hào)用一個(gè)特殊的AR（2）模型表示： P個(gè)正弦波組合的模型用2P階差分方程描述 21( )()Pkkx na x nk Z變換 2111PkkkH za z 211PkkkA za z 對(duì)于多項(xiàng)式該方程有2P個(gè)根，這些根在z平面的單位圓上，它們是 iizjei=1, 2, 3, , 2P 式中的i即是正弦波的頻率。2 2 白噪聲白噪聲中正弦波組合用一特殊的中正弦波組合用一特殊的ARMAARMA模型模型表示表示)()sin()()()(1nwnqnwnxnyPiiii 21( )( )( )()Pkkx ny nw nx na x nk 將代入21221

24、1( )- ( )-()()( )()( )()PiiPPiiiiy n w nay niw niy na y niw na w ni即上述差分方程是一特殊的ARMA（2p,2p）模型，其AR和MA參數(shù)相同，這種對(duì)稱性是白噪聲中正弦信號(hào)的一個(gè)特點(diǎn)。白噪聲中正弦波組合的信號(hào)為 受到加性白噪聲干擾的AR(p)模型等價(jià)于一個(gè)ARMA(p,p)模型。3 3 特征特征分解法譜估計(jì)分解法譜估計(jì) YTA=WTA 式中 TT221T)2(,),1(),(, 1 )2(,),1(),(pnwnwnwWaaaApnynynyYpPiiPiiinwanwinyany2121)()()()(上式的矩陣形式 用向量Y左

25、乘上式, 并取數(shù)學(xué)期望, 得到 E EYY YY TA A=EYWYW TA 式中 yyyyyyyyyyyyyyyyyyyyRrprprprrrprrrYYE)0() 12()2() 12()0() 1 ()2() 1()0(TIWWEWWXEYWEw2TTT)( 特征向量A的求解：假定W是零均值白噪聲信號(hào)，X是確定性信號(hào)。AARwyy2式中,w2是自相關(guān)矩陣Ryy的特征值; A是對(duì)應(yīng)w2的特征向量。 對(duì)于白噪聲中含有p個(gè)正弦信號(hào)的情況，當(dāng)自相關(guān)矩陣Ryy的尺寸等于或大于(2p+1) (2p+1)的時(shí)候，方差w2是自相關(guān)矩陣Ryy的最小特征值。20yywRI A 求解各正弦波的頻率：特征方程為

26、 0122221120PPPkkkzazazaza該方程有2P個(gè)根，這些根在z平面的單位圓上，它們是 iizjei=1, 2, 3, , 2P 式中的i即是正弦波的頻率。4 Pisarenko4 Pisarenko諧波分解諧波分解mPmriPiiyycos)(1m0 式中， 22iiqP ，qi是頻率為i的正弦波的幅度，Pi是其功率。對(duì)于白噪聲中P個(gè)實(shí)正弦組合信號(hào), 它的自相關(guān)函數(shù)為 PiiwyyPr12)0( 求各正弦波的幅值和功率： 12112212coscoscos12cos2cos2cos2coscoscosPyyyyPyyPPrPrPFPrPPPPPFP=r 由此，可以得到其矩陣表征形式如下： 求解噪聲功率: piiyywPr12)0( 以上就是皮薩論科譜分解法的全過程。其中關(guān)鍵的一步是求自相關(guān)函數(shù)R Ryy的最小特征值及對(duì)應(yīng)的特征向量。 Pisarenko方法的計(jì)算過程， 具體計(jì)算步驟歸納如下： (1) 由y(n)的觀測數(shù)據(jù)估計(jì)出2p+1個(gè)相關(guān)函數(shù)ryy(m), m=0, 1,2, ,2p值。 (2) 求Ryy的最小特