周期圖法估計功率譜

1、隨機信號譜估計方法的Matlab實現(xiàn)摘要：功率譜估計是隨機信號分析中的一個重要內容。從介紹功率譜的估計原理入手分析經典譜估計和現(xiàn)代譜估計兩類估計方法的原理、各自特點及在Matlab中的實現(xiàn)方法。經典功率譜估計的方差大、譜分辨率差，分辨率反比于有效信號的長度，但現(xiàn)代譜估計的分辨率不受此限制。給出了功率譜估計的應用。關鍵詞：功率譜估計；周期圖法；AR參數(shù)法； 1 引言在一般工程實際中，隨機信號通常是無限長的，例如，傳感器的溫漂，不可能得到無限長時間的無限個觀察結果來獲得完全準確的溫漂情況，即隨機信號總體的情況，一般只能在有限的時間內得到有限個結果，即有限個樣本，根據(jù)經驗來近似地估計總體的分布。有時

2、，甚至不需要知道隨機信號總體地分布，而只需要知道其數(shù)字特征，如均值、方差、均方值、相關函數(shù)、功率譜的比較精確的情況即估計值。功率譜估計(PSD)是用有限長的數(shù)據(jù)估計信號的功率譜，它對于認識一個隨機信號或其他應用方面都是重要的，是數(shù)字信號處理的重要研究內容之一。功率譜估計可以分為經典譜估計(非參數(shù)估計)和現(xiàn)代譜估計(參數(shù)估計)。2 .平均周期圖法和平滑平均周期圖法對于周期圖的功率譜估計, 當數(shù)據(jù)長度N 太大時, 譜曲線起伏加劇, 若N 太小, 譜的分辨率又不好,因此需要改進。兩種改進的估計法是平均周期圖法和平滑平均周期圖法。(1) Bartlett 法: Bartlett 平均周期圖的方法是將N

3、 點的有限長序列x(n)分段求周期圖再平均。Matlab 代碼示例1:fs=600;n=0:1/fs:1;xn=cos(2*pi*20*n)+3*cos(2*pi*90*n)+randn(size(n);nfft=512;window=hamming(nfft); %矩形窗noverlap=0;%數(shù)據(jù)無重疊p=0.9;%置信概率Pxx,Pxxc=psd(xn,nfft,fs,window,noverlap,p);index=0:round(nfft/2- 1);k=index*fs/nfft;plot_Pxx=10*log10(Pxx(index+1);plot_Pxxc=10*log10(P

4、xxc(index+1);figure(1)plot(k,plot_Pxx);figure(2)plot(k,plot_Pxx plot_Pxx- plot_Pxxc plot_Pxx+plot_Pxxc);matlab調試圖下圖(2) Welch 法:Welch 法對Bartlett 法進行了兩方面的修正, 一是選擇適當?shù)拇昂瘮?shù)w(n), 并在周期圖計算前直接加進去, 加窗的優(yōu)點是無論什么樣的窗函數(shù)均可使譜估計非負。二是在分段時, 可使各段之間有重疊,這樣會使方差減小。Matlab 代碼示例4:Fs=600;n=0:1/Fs:1;xn=cos(2*pi*40*n)+3*cos(2*pi*90

5、*n)+randn(size(n);nfft=512;window=boxcar(100);%矩形窗window1=hamming(100);%海明窗window2=blackman(100);%blackman 窗noverlap=20; %數(shù)據(jù)無重疊range='half' %頻率間隔為0 Fs/2, 計算一半的頻率Pxx,f=pwelch(xn,window,noverlap,nfft,Fs,range);Pxx1,f=pwelch(xn,window1,noverlap,nfft,Fs,range);Pxx2,f=pwelch(xn,window2,noverlap,n

6、fft,Fs,range);plot_Pxx=10*log10(Pxx);plot_Pxx1=10*log10(Pxx1);plot_Pxx2=10*log10(Pxx2);figure(1)plot(f,plot_Pxx);title('加矩形窗');figure(2)plot(f,plot_Pxx1);title('加海明窗');figure(3)plot(f,plot_Pxx2);title('加blackman 窗');matlab調試如下(3) AR模型法經典譜的主要缺點是頻率分辨率低。這是由于周期圖法在計算中把觀測到的有限長的N個數(shù)據(jù)

7、以外的數(shù)據(jù)認為是零，這顯然與事實不符。如果把以觀察到的為N個以外的數(shù)據(jù)全為零，就有可能克服經典譜估計的缺點。一個實際中的隨機過程總是可以用以下模型很好的表示：H(z)=i=1pbiz-i1+k=0panz-k當除b0外的所有bi，均為零時的形式稱為p階自回歸模型即AR模型，又稱為全極點模型。當方差為2的白自噪聲通過AR模型時，輸出的功率譜密度為：P()=21+k=0pake-jk在Matlab 仿真中可調用Pburg函數(shù)直接畫出基于burg 算法的功率譜估計的曲線圖。Matlab 代碼示例:用周期圖法求出的功率譜曲線和burg 算法求出的AR 功率譜曲線( p=50)fs=200;n=0:1/

8、fs:1;xn =cos (2*pi*40*n)+cos (2*pi*41*n)+3*cos (2*pi*90*n)+0.1*randn(size(n);window=boxcar(length(xn);nfft=512;pxx,f=periodogram(xn,window,nfft,fs);subplot(121)plot(f,10*log10(pxx)xlabel(' frequency(hz)' );ylabel(' power spectral density(Db/Hz)' );title(' periodogram psd estimate

9、' );order1=50;range=' half' ;magunits=' db' ;subplot(122)pburg(xn,order1,nfft,fs);周期圖法求出的功率譜曲線和burg 算法求出的AR 功率譜曲線( p=50) 如圖所示(4) 常見譜估計法的比較通過實驗仿真可以直觀地看出以下特性：a. 平均周期圖法和平滑平均周期圖法的收斂性較好, 曲線平滑, 估計的結果方差較小, 但是功率譜主瓣較寬, 分辨率低。這是由于對隨機序列的分段處理引起了長度有限所帶來的Gibbs 現(xiàn)象而造成的。b. 平滑平均周期圖法與平均周期圖法相比, 譜估值比較平滑, 但是分辨率較差。其原因是給每一段序列用適當?shù)拇翱诤瘮?shù)加權后, 在得到平滑的估計結果的同時, 使功率譜的主瓣變寬, 因此分辨率有所下降。c. 經典功率譜估計的分辨率反比于有效信號的長度, 但現(xiàn)代譜估計的分辨