【正文】 2. AIC 準則（信息論準則。粗略的講，AR譜的方差反比于數(shù)據(jù) )(nxN 的長度 N和信噪比 SNR[28]。因為 ? ?? ? 121 ???? ?? ?? ?? deP eP jARjx ，整個積分值為 1，所以又可得到譜匹配的局部性質(zhì)，也就是說 ? ??jAReP 是 ? ??jxeP 的包絡的一個好的近似。 模型譜的匹配性質(zhì) 由維納 辛欽定理及 AR 譜對應的一個無限長的自相關(guān)函數(shù)可得 ? ? ? ?? ?2???? j jjx eA ePeP ? ? ? ? ?2m in?? ? jjAR eAeP ? 若用 AR譜去匹配信號的譜，則誤差系列的譜應由常數(shù)譜來匹配，體現(xiàn) ??zA 的 32 白化性質(zhì)， 從整體上看 ? ?jARPe?將均勻地和 ? ?jxPe?相跟隨。 譜的分辨率 AR 譜估計的頻率分辨率，要優(yōu)于經(jīng)典譜估計方法。 AR 模型譜估計的實現(xiàn)及性質(zhì) 譜估計的步驟 我們可以按如下步驟估計 ??nx 的功率譜： ① 由 ??nxN 估計 ??nx 的自相關(guān)函數(shù)，得 ?xr ， m=0,1,....p； ②用 ?xr 代替上述地推算法中的 ??mrx ，重新求解 YuleWalker 方程，這時求出的 AR 模型參數(shù)是真實參數(shù)的估計值，即 ???1a 、 ???2a … ???pa 、 ?p? ； 將這些參數(shù)代入（ 3）式，得到 ??nx 的功率譜 ? ??jxeP 的估計。同 Burg 算法一樣，改進協(xié)方差算法進行功率譜估計時令前后向預測誤差功率之和最小，即對 ??neb 、 ??nef 前后都不加窗，但得到的協(xié)方差矩陣不是 Toeplitz 矩陣，因此正則方程不能用 Levinson 遞推算法求解。 計算步驟如下：①由初始條件 ? ? ? ? ? ?nxene nbb ?? 00 ，再由式? ? ? ?? ?? ? ? ?pmnenenenek NmnNmnbmfmNmnbmfm,.. .,2,11121 1 21211111 ??????? ?????????????? (47) 求出 ?1k ； ②由 ? ? ? ? 21010 ???? ? Nnx nxNr 得 1?m 時的參數(shù) 11,1 ?? ?ka , ? ?01 211??? ???????? ??xrk? ③求出 ??nef1 和 ??neb1 ，再 ?mk 由估計 ?2k ； ④依照 Levinson 遞推關(guān)系，求 2?m 時的參數(shù) , 2,21,2 ?? aa 及 2?? ； ⑤重復上述過程，直到 pm? ，求出所有階次時的 AR 參數(shù)。用 Burg算法進行功率譜估計時令前后向預測誤差功率之和最小，即對 ??neb 、 ??nef 前后都不加窗，使用 Levinson— Durbin 遞推可快速的求解 AR系數(shù)。因此自相關(guān)法也是已知所有 AR 系數(shù)求解方法中簡單的一種，但譜分辨率相對較差 [26]。模型參數(shù)算法就是基于 上述最小均方誤差時由模型參數(shù)估計信號功率的方法，主要有以下幾種經(jīng)典算法： 30 ? ? ? ? ? ?knxkanx pk f ???? ??1 ? ? ? ? ? ?fe n x n x n??? ? ?? ?2neEp ff ? ? ? ? ? ? ?knxkanx pk b ??? ??? 1 ? ? ? ? ? ?pnpnxne xb ???? ? ? ?? ?2neEp fbb ? 自相關(guān)法 (BT 法 )。 根據(jù)線性預測理論知：一個 P 階 AR 模型的 1p? 個參數(shù)同樣可用來構(gòu)成 P階的最佳線性預測 器，其預測的最小均方誤差等于 AR 模型激勵白噪聲的能量，即 AR 模型是在最小方差意義上對數(shù)據(jù) 的擬合。 AR 模型參數(shù)求解的典型算法 用線性方程組的常用解法 (例如高斯消元法 )解 Yule— Walker 方程，需要的運算量數(shù)量級為 3p ，但若利用系數(shù)矩陣的對稱性和 Toeplitz 性質(zhì)，則可構(gòu)成一些高效算法， Levinson— Durbin 算法是其中最著名、應用最廣泛的一種，這種算法的運算量數(shù)量級為 2p 。 將方程（ 41）兩邊同乘以 ? ?mnx ? ，并求得均值，最后得到 ? ? ? ? ? ?mrkmramr xuxpk kx ???? ?? 1 (44) 又因為 ? ? ? ?? ? ? ???? ???? 00 00 2 mh mmnxnuE ? (45) 29 由 Z 變換的定義， ? ? ? ?lim 0z H z h?? ?，在（ 42）式中，當 ??z 時，有 ?? 10?h ，綜合（ 44）和（ 45）兩式，有 ? ?? ?? ?? ? ? ???????????????????pkxkxpkkxmkramkmramr12101? (46) 在上面的推導中，應用了自相關(guān)函數(shù)的偶對稱性。 AR 模型又稱為自回歸模型，它是一個全極點模型，其當前輸出是現(xiàn)在輸入和過去輸入的加權(quán)和，表示如下 (其中 ??un為白噪聲序列； p 為 AR 模型的階數(shù) ): ? ? ? ? ? ?nuknxanX pk k ???? ??1 (41) ? ? ? ?kpk kzazAzH??????11112 (42) 由隨機信號通過線性系統(tǒng)理論知輸出序列的功率譜 ? ? 2121 ????? pkkjkjxeaep?? ? (43) 其中 2? 為白噪聲序列的方差，因此進行功率譜估計，必需求得 AR 模型的參數(shù) ka (k=l， 2? p)及 2? 。本章將會詳細的討論AR 模型參數(shù)的計算、譜的性能及其他算法（如線性預測、最大熵譜估計等）的關(guān)系，最后簡要給出 MA 模型及 ARMA 模型譜估計算法 [23]。不難想象， AR 模型易于反映譜中的峰值，而 MA 模型易于反映譜中的谷值。 工程實際中所遇到的功率譜大體分為三種，一種是“平滑”，即白噪聲的譜，另一種是“線譜”，這是由一個或多個純正弦信號所組成的信號的功率譜，這兩種是極端的情況；介于二者之間的是既有峰值又有谷值，這種譜稱為 ARMA 譜。輸出序列 ??xn可以是平穩(wěn)的隨機序列，也可以是確定性的時間序列。 （ 3）由 ??zH 的參數(shù)來估計 )(nx 的功率譜。參數(shù)模型法是現(xiàn)代譜估計的主要內(nèi)容，也是本章討論的重點，參數(shù)模型法的思路如下： （ 1）假定所研究的過程 )(nx 是由一個輸入序列 ??nu 激勵一個線性系統(tǒng) ??zH的輸出。當然，這種假定是不符合實際的，正是由于這些不符合實際的假定產(chǎn)生了經(jīng)典譜估計較差的分辨率。方差性能差的原因是無法實現(xiàn)功率譜密度原始定義中的求均值和求極限的運算。 3）弱信號被強信號的旁瓣淹沒； 4）需采用某種平滑或平均措施以改善譜估計的統(tǒng)計特性； 5）某些加窗的相關(guān)函數(shù) 會使功率譜估計值出現(xiàn)負值 [22]。矩形窗的頻譜主瓣不是無限窄的，且有旁瓣存在，這將導致能量向旁瓣中“泄漏”，主瓣變得模糊不清。 26 它們的主要缺點有： 1）頻率分辨率（區(qū)分兩個鄰近頻率分量的能力）不高。 歸結(jié)起來，自相關(guān)法和周期圖法的主要優(yōu)點是； （ 1） 計算量小 ： （ 2） 功率譜估計值正比于正弦波信號的功率。因此必須采取其他措施改善它的估計性能，才能使得這種方法有實用價值。 本章小結(jié) 根據(jù)相關(guān)函數(shù)和功率譜密度函數(shù)之間為傅立葉變換對的維納 欣欽關(guān)系，可從相關(guān)函數(shù)估計值經(jīng)傅立葉變換后得到功率譜的估計 值。 plot(k,plot_Pxx)。 k=index*Fs/nfft。 Pxx=abs(CXk)。)。 cxn=xcorr(xn,39。 %產(chǎn)生含有噪聲的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n))。 Fs=1000。 以下是利用 matlab 對自相關(guān)法進行仿真，并得到仿真圖（ 37)。對其改進的主要方法有兩種，即平均和平滑，平均就是將截取的數(shù)據(jù)段 ??nxN 再分成 L 個小段，分別計算功率譜取功率譜的平均，這種方法使估計的方差減少，但偏差加大，分辨率下降。 0 50 100 150 200 250 300 350 400 450 500 3 5 3 0 2 5 2 0 1 5 1 050 圖 34 加矩形窗仿真圖 23 0 50 100 150 200 250 300 350 400 450 500 3 5 3 0 2 5 2 0 1 5 1 050 圖 35 加漢明窗仿真圖 0 50 100 150 200 250 300 350 400 450 500 3 5 3 0 2 5 2 0 1 5 1 05 圖 36 加 blackman 窗仿真圖 24 從圖中可以看出，采用周期圖法估計得到的功率譜很不平滑，相應的估計協(xié)方差比較大，而且采用增加采樣點的辦法也不能使周期圖變得更加平滑，這是周期圖法的缺點。 pause。 pause。 plot_Pxx2=10*log10(Pxx2)。 22 plot_Pxx=10*log10(Pxx)。 [Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range)。half39。 %blackman 窗 noverlap=20。%矩形窗 window1=hamming(100)。 nfft=1024。 n=0:1/Fs:1。程序如下： clear。 figure(2) 20 plot(k,[plot_Pxx plot_Pxxplot_Pxxc plot_Pxx+plot_Pxxc])。 figure(1) plot(k,plot_Pxx)。 plot_Pxx=10*log10(Pxx(index+1))。 index=0:round(nfft/21)。 %數(shù)據(jù)無重疊 p=。 window=boxcar(length(n))。 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n))。 Fs=1000。 19 0 50 100 150 200 250 300 350 400 450 500 6 0 5 0 4 0 3 0 2 0 1 0010 圖 31 周期圖法仿真圖 以下是利用 matlab 對平均周期圖法進行仿真，并得到仿真圖（ 32）和（ 33）。%加矩形窗 [Pxx,f]=periodogram(xn,window,nfft,Fs)。 nfft=1024。%采樣頻率 n=0:1/Fs:1。程序如下： clear。 直接法和間接法的關(guān)系 由以下式子可以得出： 18 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?nkjNk Nnj w nNnkjNkjjjN enxkxNenxekxNeXeXNeXNI ???? ??????? ? ??? ???? ?????? 10 1010102 11*11 m=kn，則 k=n+m，得 ? ? ? ? ? ? ? ? ? ? ? ? mjN Nm NmjmN nN NmN emRemnxnxNI ??? ?? ?????? ?? ??? ??? ??? 1 11 01 1 ]1[ 可見周期圖法就是 BT法中 M=N1時的功率譜估計 [19]。 BT 法功率譜估計采用有偏自相關(guān)函數(shù)估計法。下面介紹自相關(guān)函數(shù)的估計。因為這兩個目的不可能同時實現(xiàn)，所以必須調(diào)整 L 在方差和偏差之間進行折衷 [18]。 顯然，為了獲得最大分辨率， L將盡可能選得大一些或就選 1??NL ，此時，平均周期圖就成為標準周期圖估計。 因為 ?? ?? 2/1 2/1 )()()](?[ ??? dPfWfPE xxBA V P E R，所以方差將減少 K倍。假定在區(qū)間 10 ??? Ln 上有 k 組獨立記錄數(shù)據(jù)，并且都是同一隨機過程的現(xiàn)實。這樣得到的均值，其方差將是用一組數(shù)據(jù) 得 到的均值的方差的 1/L。 從以上的分析可知，數(shù)據(jù)加窗用于周期圖譜估計可以降低譜估計值的旁瓣，但要降低譜估計的分辯率，而用數(shù)據(jù)加窗的辦法不能減小估計方差，因而無法降低分 辨率。但旁瓣的降低必然使主瓣加寬，而且降低了分辨率。顯然它不等于功率譜的真實值，因而是有偏估計。因此提出了周期圖的改進方法 ： 周期圖改進的方法之一是將長度為 N 的序列 )(nx 乘以同一長度的數(shù)據(jù)窗)(nw 。此外，如果序列是由多個正弦波信號組成的，而各分量強度不等，則弱信號分量可能淹沒在強信號譜的旁瓣中而無