【正文】 實際應(yīng)用中所遇到的隨機過程常?？梢钥闯捎砂自肼?[]n? 經(jīng)意線性系統(tǒng)形成的。 基于 MATLAB 的譜估計實現(xiàn) 21 4 AR 模型法 AR 模型法的基本理論 任何具有功率譜密度的隨機信號都可以看成由一白噪聲 []n? 激勵一物理網(wǎng)絡(luò)所形成。 同周期圖法一樣， BT法也存在譜分辨率的問題。 仿真圖如圖 38所示 0 50 100 150 200 25020253035404550550 50 100 150 200 25025303540455055 圖 38（采樣點數(shù)為 500） 圖 39(采樣點數(shù)為 5000) 從前述討論可知， BT法可通過增加采樣點數(shù) N 來使估計的方差減小，從而使性能得到改進。 figure(1)。 k=index*Fs/nfft。 %計算序列的自相關(guān)函數(shù) CXk=fft(cxn,nfft) Pxx=abs(CXk)。 nfft=512。%采樣頻率 n=0:1/Fs:1。?[]XRm作為自相關(guān)函數(shù)的估計，而較少使用 ?[]XRm。?v a r [ [ ] ] v a r [ [ ] ] v a r [ [ ] ]X X XNmR m R m R mN? ? ??????? 由于 39。?lim [ [ ] ] [ ] 0XXN E R m R m?? ??， 39。?[]XRm是 []XRm的有偏估計；若 m 取有限值，則當(dāng) N?? 時，39。 上 式 表明 ： 當(dāng) 0m? 時，由于39。01? ?[ ] [ ] [ ] [ ] , 1NmXXnNmR m X n X n m R m m NNN????? ? ? ? ?? 其均值 39。 基于 MATLAB 的譜估計實現(xiàn) 19 式 (36)估計自相關(guān)函數(shù)的方法，雖然當(dāng) N 較大且 Nm?? 時能得到一致估計，但當(dāng) m 接近 N 時， ?[]XRm的估計方差變大，因而不能得到有用的估計。依此類推，對不同的 l 的情況，出現(xiàn)的次數(shù)將 為 ()N m l??，于是上式可寫成 1 22 ( 1 )1?v a r [ [ ] ] { ) [ [ ] [ ] [ ] }()NmX X X Xl N mR m N m l R l R l m R l mNm??? ? ? ?? ? ? ? ? ?? ? 1 22 ( 1 )1 1 [ [ ] [ ] [ ] ]()NmX X Xl N m ml R l R l m R l mN m N??? ? ? ?? ? ?? ? ? ? ??????? 1 22 ( 1 ) [ [ ] [ ] [ ] ]()NmX X Xl N mN R l R l m R l mNm??? ? ? ?? ? ? ?? ? 當(dāng) Nm?? 時，上式以 1N 趨于零。 其次求 ?[]XRm的方差，根據(jù)方差的定義 22? ? ?v a r [ [ ] ] [ [ ] ] [ [ ] ]X X XR m E R m E R m?? (37) 由式（ 36）可知，上式等號右邊第一項可表示為 1122 001?[ [ ] ] { [ ] [ ] [ ] [ ] }()N m N mX nkE R m E X n X n m X k X k mNm? ? ? ???? ? ?? ?? (38) 當(dāng) []Xn是零均值正態(tài)序列時，它的各高階矩都可以用其一階和二階矩來表示。 首先求 ?[]XRm的均值，有 101?[ [ ] ] { [ ] [ ] }NmXnE R m E X n X x mNm?????? ? 101 ? [ ] [ ]Nm XXn R m R mNm?????? ? 即自相關(guān)函數(shù)的估計值 ?[]XRm的數(shù)學(xué)期望等于序列 []Xn的自相關(guān)函數(shù)的真實值。 2，求 ?[]XRm的離散傅里葉變換，得到 []Xn功率譜的估值： ( ) [ ] , 1M jmXXmMG R m e M N??? ???? ? ?? 由于式 (36)僅利用了 []Xn的 N 個有限值得到自相關(guān)函數(shù)的估計 ?[]XRm，因此它與信號本身的自相關(guān)函數(shù) []XRm會有一定程度的差別。 BT 法 BT 法 BT 法首先計算隨機信號 []Xn的自相關(guān)函數(shù) [ ] [ ( ) ( ) ]XR m E X n X n m??，若[]Xm Rm????? ??? ，可由 []XRm的傅里葉變換的得到它的功率譜，即： ( ) [ ] jmXXmG R m e ???? ????? ? BT 法的估計步驟如下： 基于 MATLAB 的譜估計實現(xiàn) 17 1，先估計隨機信號 []Xn的自相關(guān)函數(shù)。 在采用周期圖法進行功率譜估計時，可以將平均周期圖法和加窗周期圖法這兩種方法結(jié)合在一起，效果會更好。但是在實際應(yīng)用中，幅度相差 100倍甚至更高的正弦分量同時存在于一個信號中是可能的，這種情況下，由于矩形窗的旁瓣幅度較大，難以區(qū)分這樣的兩個正弦分量，故估計的效果較差，故還應(yīng)考慮旁瓣幅度較小的窗函數(shù)。對圖 36的主瓣進行放大后如圖 37所示 . 基于 MATLAB 的譜估計實現(xiàn) 16 圖 37 從圖 37可以看出，只有運用矩形窗和凱撒窗才能分辨出這兩個頻率相差較小的信號。 仿真圖如下： 圖 36 從圖 36可以看出，采用 blackman窗和 hamming窗 后，功率譜的主瓣寬度比較窄，頻率泄露的情況相對不太明顯（也就是除主瓣以外的譜 線的幅度較主瓣幅度相比較?。?。triang窗 39。 figure(8) plot(f7,10*log10(Pxx7))。kaiser窗 39。 figure(7) plot(f6,10*log10(Pxx6))。hann窗 39。 figure(6) plot(f5,10*log10(Pxx5))。hamming窗 39。 figure(5) plot(f4,10*log10(Pxx4))。blackman窗 39。 figure(4) plot(f3,10*log10(Pxx3))。bartlett窗 39。 figure(3) plot(f2,10*log10(Pxx2))。boxcar窗 39。 figure(2) plot(f1,10*log10(Pxx1))。)。title(39。 subplot(337)。kaiser窗 39。plot(f6,10*log10(Pxx6))。)。title(39。 基于 MATLAB 的譜估計實現(xiàn) 15 subplot(335)。hamming窗 39。plot(f4,10*log10(Pxx4))。)。title(39。 subplot(333)。bartlett窗 39。plot(f2,10*log10(Pxx2))。)。title(39。 figure(1) subplot(331)。 [Pxx6,f6]=periodogram(xn,window6,nfft,Fs)。 [Pxx4,f4]=periodogram(xn,window4,nfft,Fs)。%直接法 [Pxx2,f2]=periodogram(xn,window2,nfft,Fs)。 window7=triang(length(xn))。 window5=hann(length(xn))。 window3=blackman(length(xn))。 window1=boxcar(length(xn))。%產(chǎn)生含有噪聲的序列 xn=cos(2*pi*40*n)+cos(2*pi*41*n)+3*cos(2*pi*90*n)+*randn(size(n))。 加窗周期圖的仿真 采用加窗周期圖方法進行功率譜估計的仿真程序如下： Fs=600。 figure(1) plot(k,plot_Pxx) 數(shù)據(jù)重疊方法的仿真圖如圖 35 所示。 k=index*Fs/nfft。 [Pxx,f]=pwelch(xn,window,noverlap,nfft)。 window=boxcar(100)。 xn=cos(2*pi*40*n)+5*cos(2*pi*90*n)+*randn(size(n))。 基于 MATLAB 的譜估計實現(xiàn) 14 程序如下： Fs=600。 end。 sx=sx+p。 for k=1:K [p,f]=periodogram(xn(n1:(n1+L1)))。 f=0。 function [sx,f]=bart(xn,K) L=floor(length(xn)/K)。%k相當(dāng)于頻率 ,256是 PERIODOGRAM（在 bart函數(shù)中調(diào)用的）函數(shù)中默認的 nfft的值 plot_Pxx=10*log10(Pxx(index+1))。 index=0:round(256/2 1)。 K=10。 %采樣頻率 n=0:1/Fs:1。如圖 41，真實的譜線應(yīng)該是在頻率為 40Hz 和 90Hz 處為兩個沖激，而周期圖法估計出的功率譜在這兩個頻率出的主瓣非常寬，失真較明顯。 0 50 100 150 6 0 5 0 4 0 3 0 2 0 1 001020 圖 33 由前述討論可知，在周期圖方法中，無論怎樣增大數(shù)據(jù)采樣點數(shù) N ，功率譜估計的方差一直 不為零。 [Pxx,f]=periodogram(xn,window,nfft,Fs)。 xn=cos(2*pi*40*n)+5*cos(2*pi*90*n)+*randn(size(n))。 仿真程序如下： Fs=600。 基于 MATLAB 的譜估計實現(xiàn) 12 采用加窗周期圖方法的數(shù)學(xué)表達式為 101 [ ] [ ]N jw nn X n n eN ?? ?????? 其中， []n? 為窗函數(shù)。 幾種常用窗函數(shù)的數(shù)學(xué)表達式如下： 矩形窗： ( ) 1, 0n n N? ? ? ? 三角窗（ bartlett） : 2 / ( 1 ) , 0()2 2 / ( 1 ) , ( 1 ) / 2n N n Nn n N N n N? ? ? ??? ? ? ? ? ? ?? 漢寧窗： 12( ) 1 c o s , 021 nn n NN ?? ????? ? ? ?????????? 哈明窗： 2( ) 0 . 5 4 0 . 4 6 c o s , 01nn n NN ?? ??? ? ? ?????? 布萊克曼窗： 24( ) 0 . 4 2 0 . 5 c o s 0 . 0 8 c o s , 011nnn n NNN??? ? ? ? ?? ? ? ? ?? ? ? ???? ? ? ? 在 MATLAB中還給出了幾種常用的窗函數(shù)，他們的時域和頻域特性分別如下： 基于 MATLAB 的譜估計實現(xiàn) 11 5 0 0 5000 . 51b o x c a r 窗的時域特性 5 0 0 5000 . 51b a r t l e t t 窗的時域特性 5 0 0 5000 . 51b l a c k m a n 窗的時域特性 5 0 0 5000 . 51h a m m i n g 窗的時域特性 5 0 0 5000 . 51hann 窗的時域特性 5 0 0 5000 . 51k a i s e r 窗的時域特性 5 0 0 5000 . 51t r i a n g 窗的時域特性 圖 31 5 0 0 0 5000204060b o x c a r 窗的頻域特性 5 0 0 0 5000102030b a r t l e t t 窗的頻域特性 5 0 0 0 5000102030b l a c k m a n 窗的頻域特性 5 0 0 0 5000102030h a m m i n g 窗的頻域特性 5 0 0 0 5000102030hann 窗的頻域特性 5 0 0 0 50002040k a i s e r 窗的頻域特性 5 0 0 0