【正文】 函數(shù) )(krxx 的另一種估計方法可表示為 ? ??? ?? 10 )()(1)(? kN nxx nxknxNkr （ 31） )()]()([1)](?[ 10 krN kNnxknxENkrE xxkN nxx ? ?? ? ???? 可見，它是相關(guān)函數(shù) )(krxx 的有偏估計。 當(dāng) ??N 時， 0)}(?var{ ?krxx ，即（ 31）式的 )(? krxx 也是 )(krxx 的一致估計。 功率譜密度的另一個定義 : 可以證明，功率譜密度（ PSD）的一個近似等效的定義是 })2(e xp)(12 1{l i m)( 2fnjnxMEfP M MnMxx ???? ????? （ 32） 上式定義的 PSD 與維納一辛欽定理 ?????? ?? k xxxx fkjkrfP )2e xp ()()( ?（ 33） 是等效的。 15 周期圖法 周期圖法的定義 周期圖法，它是把隨機序列 ??xn的 N 個觀測數(shù)據(jù)視為一能量有限的序列，直接計算 ??xn的離散傅立葉變換，得 ??Xk， 然后再取其幅值的平方，并除以 N，作為序列 ??xn真實功率譜的估計。 由式（ 34）可知周期圖的均值是真實 PSD和 Bartlett 窗傅里葉變換的卷積，在平均意義上得到真實功率譜密度（ PSD）的平滑形式。 )()](?[lim fPfPE xxP E RN ??? 這是由于 )(fWB 收斂到狄拉克 ? 函數(shù)。對于白高斯過程，)()](?[ 2 fPfPE xxxP E R ?? ?，方差 )}(?var{ fPPER ])2s in2s in(1)[( 22 fN NffPxx ???? 16 對任何非白過程，只要記錄數(shù)據(jù) 足夠長 , ])2s i n2s i n(1)[()](?va r[ 22 fN NffPfP xxP E R ???? 對于不靠近 0 或 21? 的頻率，且 ??N 時，上式近似地退化為： )()](?v ar[ 2 fPfP xxPER ? 可以看出，方差與數(shù)據(jù)長度 N無關(guān)，即方差不隨 N的增大而減小至零。 上述論證表明，我們不能寄希望于直接用周期圖方法獲得良好的譜估計，必須采用適當(dāng) 的修正措施減小估計方差，才能使之成為一種實用的方法 [17]。由此可以想到，當(dāng)用周期圖方法作譜估計時，它的譜分辨率約與 N 成反比，且和信號本身的特征，例如信噪比等無關(guān)。這種所謂信號能量（向旁瓣）泄漏現(xiàn)象如果不設(shè) 法消除，也將妨礙周期圖譜估計法的應(yīng)用。 數(shù)據(jù)加窗后的周期圖譜估計值的數(shù)學(xué)期望值等于譜的真實值與窗譜函數(shù)的平方的卷積。 若序列為單頻信號，則 )(fPxx 為 ? 函數(shù)，這樣，數(shù)據(jù)加窗后的譜估計值的均值與窗譜函數(shù)的平方 形狀相同，因此選用低旁瓣的數(shù)據(jù)窗可使得雜散響應(yīng)減少。 數(shù)據(jù)加窗后，周期圖譜估計值的方差大于或近似等于譜估計值的平方，且不隨數(shù)據(jù)長度的增大 而減小到 0。 平均周期圖 平均周期圖的思想：對一個隨機變量進行觀測，得到 L 組獨立記錄數(shù)據(jù)， 17 用每一組數(shù)據(jù)求其均值，然后將 L 個均值加起來求平均。 為了改進周期圖的統(tǒng)計特性，可以近似地用對一組周期圖進行平均的方法完成期望運算。 平均周期圖估計器定義為： ???? 10 )( )(?1)(? Km mP E RA V P E R fPKfP 210)( )2e xp ()(1)(? ??? ??Ln mmP E R fnjnxLfP ? 其中 )(? )( fPmPER 是第 m 個數(shù)據(jù)組的周期圖。 又因為 )](?va r[1)](?va r[ )( fPKfP mP E RA V P E R ? ，所以譜估計器具有 L/1 樣本的分辨率。但有時為了減小方差，應(yīng)該把 LNK /? 選得大一些，或等效地把 L取得小一些。 自相關(guān)法 自相關(guān)法是先估計自相關(guān)函數(shù)，再進行傅里葉變換得到功率譜。 利用觀測到的實隨機序列 ??nX ，估計自相關(guān)函數(shù)的兩種方法是：無偏自相關(guān)函數(shù)估計和有偏自相關(guān)函數(shù)估計。自相關(guān)法的理論基礎(chǔ)是維納 辛欽定理，即先由 ??nxN 估計出自相關(guān)函數(shù)?)(mr ，然后對 ?)(mr 求傅里葉變換得到 ??nxN 的功率譜，因為由這種方法求出的功率譜是通過自相關(guān)函數(shù)間接得到的，所以稱為間接法。 譜估計仿真與比較 以下是利用 matlab 對周期圖法 進行仿真，并得到仿真圖（ 31）。 Fs=1000。%產(chǎn)生含有噪聲的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n))。 window=boxcar(length(n))。%直接法 plot(f,10*log10(Pxx))。程序如下： clear。n=0:1/Fs:1。 nfft=1024。 %加矩形窗 noverlap=0。 %置信概率[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p)。 k=index*Fs/nfft。 plot_Pxxc=10*log10(Pxxc(index+1))。 pause。 0 50 100 150 200 250 300 350 400 450 500 4 0 3 0 2 0 1 0010203040 圖 32 周期圖法仿真圖 21 0 50 100 150 200 250 300 350 400 450 500 6 0 4 0 2 00204060 圖 33 平均周期圖法仿真圖 以下是 利用 matlab 對加窗周期圖法進行仿真，并得到仿真圖（ 34），（ 35）和（ 36）。 Fs=1000。 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n))。 window=boxcar(100)。%漢明窗 window2=blackman(100)。 %數(shù)據(jù)無重疊 range=39。%頻率間隔為 [0 Fs/2]，只計算一半的頻率 [Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range)。 [Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range)。 plot_Pxx1=10*log10(Pxx1)。 figure(1) plot(f,plot_Pxx)。 figure(2) plot(f,plot_Pxx1)。 figure(3) plot(f,plot_Pxx2)。周期圖法得出的估計譜方差特性不好：當(dāng)數(shù)據(jù)長度 N 太大 時，譜線的起伏加劇 [20]； N 太小時，譜的分辨率又不好。平滑是用一個適當(dāng)?shù)拇昂瘮?shù)與算出的功率譜進行卷積，使譜線平滑，這種方法得出的譜線是無偏的，方差也小，但分辨率下降。程序如 : clear。 %采樣頻率 n=0:1/Fs:1。 nfft=1024。unbiased39。 %計算序列的自相關(guān)函數(shù) CXk=fft(cxn,nfft)。 index=0:round(nfft/21)。 plot_Pxx=10*log10(Pxx(index+1))。 25 0 50 100 150 200 250 300 350 400 450 500505101520253035 圖 37 自相關(guān)法仿真圖 當(dāng) 1MN??時，自相關(guān)法與周期圖法估計出的功率譜是一樣的：當(dāng)1MN??時，自相關(guān)法的偏差大于周期圖法，在窗函數(shù)滿足一定條件時是漸進無偏估計，方差小于周期圖的方差，分辨率比周期圖法低，是與窗函數(shù)的選擇有關(guān) [21]。從原理上說這種方法是直截了當(dāng)?shù)?，可是由此得出的譜估計值不僅有偏，而且方差很大，它還不隨數(shù)據(jù)長度的增大而減小到零。常用的兩種辦法是相關(guān)函數(shù)加窗和分段平均周期圖法，這兩種辦法都有實用價值，它們都采用了 FFT 算法，使得計算量大大減小。 （ 3）對于長記錄數(shù)據(jù)，是一種良好的實用模型。 這是因為它們的頻率分辨率約為數(shù)據(jù)長度的倒數(shù)，且與數(shù)據(jù)的特征或其信噪比無關(guān)，而實際應(yīng)用中一般不可能獲得很長的數(shù)據(jù)記錄； 2）經(jīng)典譜估計方法在工程中都是以離散傅立葉變換為基礎(chǔ)的，它隱含著對無限長數(shù)據(jù)序列進行加窗處理（加了一個有限寬的矩形窗）。嚴(yán)重時，會使主瓣產(chǎn)生很大失真，甚至主瓣中的弱分量被旁瓣中的強泄漏所淹蓋。 27 4 現(xiàn)代譜估計 平穩(wěn)隨機信號的參數(shù)模型 由上一章討論可知，經(jīng)典功率譜估計方法的方差性較差，分辨率較低。分辨率低的原因，對周期圖法是假定了數(shù)據(jù)窗以外的數(shù)據(jù)全為零，對自相關(guān)法是假定了在延遲窗以外的自相關(guān)函數(shù)全為零。 在第一章已經(jīng)簡潔的介紹了現(xiàn)代譜估計的基本方法，這些方法技 術(shù)的目標(biāo)在于努力改善譜估計的分辨率。 （ 2）由已知的 )(nx ，或其自相關(guān)函數(shù) ??mrx 來估計 ??zH 的參數(shù)。 ??Hz是一個因果的線性移不變離散時間系統(tǒng)，當(dāng)然，它應(yīng)該是穩(wěn)定的，其單位抽樣響應(yīng) ??hn是確定的。若 ??xn是確定性的，那么 ??un是一個沖激序列，若 ??xn 是隨機序列，那么 ??un應(yīng)是一個白噪聲序列。顯然，由于 ARMA 模型是一個零極點模型，它易于反映功率譜中的峰值和谷值。 AR,MA 和 ARMA 是功率譜估計中最主要的參數(shù)模型。 28 AR 模型的正則方程與參數(shù)計算 正則方程的求導(dǎo) 參數(shù)模型法功率譜估計的主要思想是：將廣義平穩(wěn)的過程 ??xn表示成一個輸入序列 ??un激勵線性系統(tǒng) ??Hz的輸出；由已知的 ??xn或其自相關(guān)函數(shù) ??mrx來估計 ??Hz的參數(shù)；由 ??Hz的參數(shù)估計 ??xn的功率譜。 假定 ??un、 ??un都是平穩(wěn)的隨機信號， ??un為白噪聲，方差為 2? ，現(xiàn)在，我們希望建立 AR 模型的參數(shù) ka 和 ??xn的自相關(guān)函數(shù)的關(guān)系，也即 AR 模型的正則方程 [24]。上式寫成矩陣形式，即 ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?2120 1 2 11 0 1 1 02 0 0 2 01 2 0 0 0x x x xx x x xx x x xx x x xr r r r pr r r r p ar r r r p ar p r p r p r??? ?????? ????????????????????????????? ?? ???? 上述兩式即是 AR模型的正則方程，又稱 YuleWalker 方程 [25]。這是一種按階次進行遞推的算法，即首先以 ? ?0AR 和 ??1AR模型參數(shù)作為初始條件，計算 ? ?2AR 模型參數(shù)；然后根據(jù)這些參數(shù)計算 ??3AR 模型參數(shù)等，一直到計算出 ? ?ARp 模型參數(shù)為止，當(dāng)整個迭代計算結(jié)束后，不僅求得了所需要的 P 階 AR 模型參數(shù)，而且還得到了所有各低階模型的參數(shù)。 “ 前向預(yù)測”是利用 n 之前的 P 個值對 ??xn性預(yù)測，如公式 (48)、 (49)、(410)所示；與之對應(yīng)的“后向預(yù)測”公式為 (411)、 (412)、 (413)，其中 ??en為為預(yù)測誤差， P 預(yù)測誤差功率， f 表示前向預(yù)測， b 表示后向預(yù)測。用自相關(guān)法進行功率譜估計，但估計時令前向預(yù)測誤差功率最小，即對 ??nef 前后都加窗構(gòu)成， Wiener— Hopf 方程系數(shù)為 Toepli tz 矩陣，使用 Levinson— Durbin 算法可方便快速的求解 AR 系數(shù)。 Burg算法。 Burg 算法是建立在數(shù)據(jù)基礎(chǔ)之上的，避免了先計算自相關(guān)函數(shù)從而提高計 算速度；是較知為通用的方法，計算不太復(fù)雜，且分辨率優(yōu)于自相關(guān)個法，但對于白噪聲加正弦信號有時會出現(xiàn)譜線分裂現(xiàn)象。 31 (1)改進協(xié)方差算法。Marple 于 1980 年提出實現(xiàn)協(xié)方差方程求解的快速算法，大大提高了譜估計的性能 [27]。 AR 模型譜估計的性質(zhì) 譜的平滑特性 AR 模型是一有理分式，估計出的譜平滑，不需要像周期圖那樣再做平滑或平均，因此，不需要為此去犧牲分辨率。其原因在于求解 AR 模型 參 數(shù) 的 過 程 ， 實 際 上 意 味 著 將 根 據(jù) ? ? ? ? ? ?1,...,1,0 ?Nxxx 估計的? ? ?? ? ?Mrrr xxx ,...,1,0 按一定準(zhǔn)則進行了外推。當(dāng)均值為 1 時， ? ?jARPe?將在 ? ?jxPe? 的上下波動，即在有的區(qū)域， ? ? ? ?jjAR xP e P e???，而在另外的區(qū)域 ，? ? ? ??? jxjAR ePeP ? 。 譜的統(tǒng)計性質(zhì) 嚴(yán)格的分析 AR 譜的方法比較困難，目前尚未有一個解析表達式。 AR 模型階次 p 的選擇 準(zhǔn)則（最終預(yù)測誤差準(zhǔn)則） ? ? ? ?? ?112 ?? ??? mN mNmFPE m? 隨著 m的增加，使 ? ?mFPE 達到最小值時的 p