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

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

2、寬帶噪聲中檢測窄帶信號。 22)()(HPyyBT法周期圖法改進的周期圖法BT法是先估計自相關(guān)函數(shù), 然后進行傅里葉變換得到功率譜。有偏自相關(guān)函數(shù)估計的誤差相對較小，是一種漸近一致估計：1|0*)()(1)(mNnxxmnxnxNmr-BT(e)( )ejjmxxmPrmn 周期圖法的定義如下: 212j-j011(e )( )eNnjxxnPx nX eNN觀 測 數(shù) 據(jù) x(n)FFT取 模的 平 方1/NPxx(ej )1. 1. 周期圖與周期圖與BTBT法的等價關(guān)系法的等價關(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法和周期圖法是等價的。 2. 周期圖法譜估計質(zhì)量分析周期圖法譜估計質(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)計平均值等于它的真值卷積三角譜窗函數(shù)，因此周期圖是有偏估計，但當(dāng)N時，wB(m)1，三角譜窗函數(shù)趨近于函數(shù)，周期圖的統(tǒng)計平均值趨于它的真值，因此周期圖屬于漸近無偏估計。 2) 周期圖的方差 用這種方法估計的功率譜在2x附近起伏很大，故周期圖是非一致估計。 白噪聲的周期圖 Bartlett平均周期圖法窗口處理法平均周期圖Welch法（修正的周期圖求平均法） 主要思想：主要思想：對序列x(n)進行L次獨立觀測或?qū)⑵浞殖蒐段，計算每組觀測數(shù)據(jù)的周期圖，再將L個周期圖加和后求平均。 LiixxIL

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

6、，M越小，偏移越大，分辨率越低；估計的均方誤差也減少； 以分辨率的降低換取了估計方差的減少，估計量的方差和分辨率是一對矛盾。平均周期圖法 Pxx(ej)/032101230321012303210123Pxx(ej)/Pxx(ej)/(a)(b)(c)N256， L2N256， L4N256， L8主要思想：主要思想：用一適當(dāng)?shù)墓β首V窗函數(shù)W(ej)與周期圖進行卷積，來達到使周期圖平滑的目的。 (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 估計方法：估計方法：那么 mxxMMmjxxemwmrPj1)1()()()e (mxxMMmxxmwmrEPEj -1)1(je )()( )e (又 1j-j(1)()( )( ) ( )eMmxxxxBmME Perm wm w m 偏移分析： n 估計效果分析：估計效果分析：| ( )( )( )( )xxxxxxBNmE rmrmrmwmN可得 周期圖的窗函數(shù)法仍然是有偏估計, 其偏移和wB(m)、w(m)兩個窗函數(shù)有關(guān)。 如果w(m)窗的寬度比較窄，M比N小得多，這樣|m|p AR譜估計方法可歸結(jié)為求解AR模型系數(shù)或線性預(yù)測器系數(shù)的問題。AR模型參數(shù)估計

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

9、，有P個預(yù)測系數(shù)，N個數(shù)據(jù)經(jīng)過沖激響應(yīng)為api(i=0,1, 2, , P)的濾波器， 輸出預(yù)測誤差e(n)的長度為N+P， 因此應(yīng)用下式計算: 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模型的各個系數(shù)以及模型輸入白噪聲方差求出后， 信號功率譜用下式計算： 21j -22j2je11| )e (|)e (PiiiwwxxaHP利用列文森遞推法計算功率譜的流程圖 輸入： x(n),n0,1,2,3,N1,2計算：)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ù)估計自相關(guān)序列，當(dāng)數(shù)據(jù)長度較短時，估計誤差會比較大，AR參數(shù)的計算就會引入很大的誤差。從而導(dǎo)致功率譜估計出現(xiàn)譜線分裂與譜峰頻率偏移等現(xiàn)象。2、 伯格（伯格（Burg）遞推法）遞推法 n 估計方法：估計方法：直接由時間序列計算AR模型參數(shù)的方法，求前、后向預(yù)測誤差平均功率最小時的反射系數(shù)kp，進而求AR模型參數(shù)ak和2w。20,ppkwpkak 設(shè)信號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ù)kp，令 0ppk1*11122112( )(1)(|( )|(1)| )Nfbppn ppNfbppn pen enkenen 基于反射系數(shù)kp，由Levinson-Durbin遞推關(guān)系求AR模型參數(shù)ak和2w，進而求得功率譜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ù)估計自相關(guān)函數(shù)的計算，適合短序列參數(shù)估計，克服了L-D遞推中的某些缺點，計算量小。但對正弦信號的譜估計，仍存在某些譜線分裂與頻率偏移現(xiàn)象。3 3、 協(xié)方差法與修正

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

15、向預(yù)測誤差平均值最小的方法， 估計AR模型的參數(shù)，進而估計信號的功率譜。,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的約束條件，估計得到的譜在譜線分裂和頻率偏移

16、時較Burg法有較大改善。 自相關(guān)法可以用Levinson遞推算法，運算量小，但分辨率受窗長度的限制；協(xié)方差法，去除了自相關(guān)法加窗處理的不合理假設(shè)，分辨率高，運算量較大；修正協(xié)方差法，分辨率高，在譜線分裂和偏移上較Burg法有較大改善，運算量大；Burg算法，可用改進的Levinson遞推算法，分辨率高，但對正弦信號存在譜線分裂和偏移現(xiàn)象。 例例4.5.1 已知信號的四個觀察數(shù)據(jù)為x(n)=x(0), x(1), x(2), x(3)=2, 4, 1, 3， 分別用自相關(guān)法和協(xié)方差法估計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對于白噪中的AR信號，其階次的選擇應(yīng)折衷考慮。如選擇AR模型，其階次應(yīng)加大，較低的階次會使譜估計產(chǎn)生偏移， 降低分辨率。信噪比愈低，平滑作用愈嚴重，愈需要高

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

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

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

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

22、參數(shù) 和 ；計算功率譜 ；(e)jxxP xxrmka2w21(e)1ejwxxNj kkkPa xxrm該方法特別適用于對多正弦加白噪聲序列進行譜分析，可以得到比AR模型法更高的分辨率和更準確的頻率估計，尤其在信噪比低時更有效。估計思想：將白噪聲加正弦波作為一特殊的ARMA模型，用特征方程求它的參數(shù)，計算出正弦波的頻率和功率，及噪聲功率。1 1 正弦波用退化正弦波用退化ARAR模型模型表示表示令隨機相位實正弦信號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 其中，正弦波信號用一個特殊的AR（2）模型表示： P個正弦波組合的模型用2P階差分方程描述 21( )()Pkkx na x nk Z變換 2111PkkkH za z 211PkkkA za z 對于多項式該方程有2P個根，這些根在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ù)相同，這種對稱性是白噪聲中正弦信號的一個特點。白噪聲中正弦波組合的信號為 受到加性白噪聲干擾的AR(p)模型等價于一個ARMA(p,p)模型。3 3 特征特征分解法譜估計分解法譜估計 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是零均值白噪聲信號，X是確定性信號。AARwyy2式中,w2是自相關(guān)矩陣Ryy的特征值; A是對應(yīng)w2的特征向量。 對于白噪聲中含有p個正弦信號的情況，當(dāng)自相關(guān)矩陣Ryy的尺寸等于或大于(2p+1) (2p+1)的時候，方差w2是自相關(guān)矩陣Ryy的最小特征值。20yywRI A 求解各正弦波的頻率：特征方程為

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