功率譜估計淺談

1、功率譜估計淺談?wù)航榻B了幾種常用的經(jīng)典功率譜估計與現(xiàn)代功率譜估計的方法原理，并利用Matlab對隨機信號進行功率譜估計，對兩種方法做出比較，分別給出其優(yōu)缺點。關(guān)鍵詞：功率譜；功率譜估計；經(jīng)典功率譜估計；現(xiàn)代功率譜估計前言功率譜估計是從頻率分析隨機信號的一種方法 ,一般分成兩大類:一類是經(jīng)典譜估計;另一類是現(xiàn)代譜估計。由于經(jīng)典譜估計中將數(shù)據(jù)工作區(qū)以外的未知數(shù)據(jù)假設(shè)為零，這相當于數(shù)據(jù)加窗，導(dǎo)致分辨率降低和譜估計不穩(wěn)定?，F(xiàn)代譜估計則不再簡單地將觀察區(qū)外的未知數(shù)據(jù)假設(shè)為零，而是先將信號的觀測數(shù)據(jù)估計模型參數(shù)，按照求模型輸出功率的方法估計信號功率譜，回避了數(shù)據(jù)觀測區(qū)以外的數(shù)據(jù)假設(shè)問題。周期圖、自相關(guān)法

2、及其改進方法(Welch)為經(jīng)典(非參數(shù))譜估計方法, 其以相關(guān)和傅里葉變換為基礎(chǔ),對于長數(shù)據(jù)記錄較適用,但無法根本解決頻率分辨率低和譜估計穩(wěn)定性的問題,特別是在數(shù)據(jù)記錄很短的情況下,這一問題尤其突出。以隨機過程的參數(shù)模型為基礎(chǔ)的現(xiàn)代參數(shù)法功率譜估計具有更高的頻率分辨率和更好的適應(yīng)性,可實現(xiàn)信號檢測或信噪分離,對語音、聲納雷達、電磁波及地震波等信號處理具有重要意義,并廣泛應(yīng)用于通信、自動控制、地球物理等領(lǐng)域。在現(xiàn)代參數(shù)法功率譜估計方法中,比較有效且實用的是AR模型法，Burg譜估計法,現(xiàn)代譜估計避免了計算相關(guān),對短數(shù)據(jù)具有更強的適應(yīng)性,從而彌補了經(jīng)典譜估計法的不足,但其也有一些自身的缺陷。下面

3、就給出這兩類譜估計的簡單原理介紹與方法實現(xiàn)。經(jīng)典譜估計法經(jīng)典法是基于傳統(tǒng)的傅里葉變換。本文主要介紹一種方法：周期圖法。周期圖法由于對信號做功率譜估計，需要用計算機實現(xiàn)，如果是連續(xù)信號，則需要變換為離散信號。下面討論離散隨機信號序列的功率譜問題。連續(xù)時間隨機信號的功率譜密度與自相關(guān)函數(shù)是一對傅里葉變換對，即：若是的抽樣序列，由序列的傅里葉變化的關(guān)系，可得即與也是一對傅里葉變換對。顯然，由序列傅里葉的頻譜特性可知是以為周期的。而實際計算只能從離散隨機信號序列x(n)的有限長(長度為N)的數(shù)據(jù)來對與進行估計。設(shè)有限長離散序列為x(n),則：由DFT的下列卷積特性：若，則：從而：即綜上所述，先用FFT

4、求出隨機離散信號N點的DFT，再計算幅頻特性的平方，然后除以N，即得出該隨機信號的功率譜估計。由于這種估計方法在把離散化的同時，使其功率譜周期化，故稱之為“周期圖法”，也稱為經(jīng)典譜估計法。周期圖法進行譜估計，是有偏估計，由于卷積的計算過程會導(dǎo)致功率譜真實值的尖峰附近產(chǎn)生泄漏，相對地平滑了尖峰值，因此造成譜估計的失真。另外，當N時，功率譜估計的方差不為零，所以不是一致性估計。并且功率譜估計在等于整數(shù)倍的各數(shù)字頻率點互不相關(guān)。其譜估計的波動比較顯著，特別是當N越大、越小時，波動越明顯。但如果N取得太小，又會造成分辨率的下降。圖1. 原始信號1圖2 原始信號1的功率譜估計圖3. 原始信號2圖4. 原

5、始信號2的功率譜估計圖5. 平均周期圖法(4*256)圖6. 平均周期圖法(重疊一半)圖1所示的信號為，其中，randn是正態(tài)分布隨機數(shù)組，N為256，t是從0到1，dt為1/256。圖2為該信號的功率譜估計。圖2所示的信號為，其中，randn是正態(tài)分布隨機數(shù)組，N為1024，t是從0到1，dt為1/1024。圖4為該信號的功率譜估計。圖5是將圖2所示的信號分為四段，每段的范圍分別為(1,256),(257,512),(513,768),(768,1024).每一道都沒有重疊。然后對分段分別作傅里葉變換，再把功率譜加起來做平均，得到圖5。圖6是將圖2所示的信號分為六段，分別為(1:256)，(

6、129:384)，(257:512)，(385:640)，(513:768)，(641:896)，(769:1024)。每兩段之間都重合一半。圖1和圖3相比，圖1較為平滑，相應(yīng)的，圖1的功率譜也比較平滑。圖5和圖6比，圖6較為平滑，這是因為圖6的譜是六段的平均。對信號加入窗函數(shù)的話，功率譜的變化也是很明顯的。圖7. 加入矩形窗原始信號和512點、1024點功率譜圖8.Bartlett平均周期圖法現(xiàn)代譜估計法現(xiàn)代參數(shù)法功率譜估計方法中,比較有效且實用的是AR模型法，Burg譜估計法，在本文中介紹的是AR模型法。 AR模型法經(jīng)典譜的主要缺點是頻率分辨率低。這是由于周期圖法在計算中把觀測到的有限長的

7、N個數(shù)據(jù)以外的數(shù)據(jù)認為是零，這顯然與事實不符。如果把已觀測到的數(shù)據(jù)估計出一白噪聲激勵，就不必認為N個以外的數(shù)據(jù)全為零，就有可能克服經(jīng)典譜估計的缺點。一個實際中的隨機過程總是可以用以下模型很好的表示：當除外的所有均為零時的形式稱為p階自回歸模型即AR模型，又稱為全極點模型。當方差為的白噪聲通過AR模型時，輸出的功率譜密度為：若已知參數(shù)及,就可以得到信號的功率譜估計。它們之間是Yule-Walker方程。解這個方程是一個復(fù)雜的數(shù)學(xué)問題，這里不做討論。圖9. 原始信號3圖10. 自相關(guān)函數(shù)的無偏估計圖11. AR模型求出的功率譜圖7所示的信號長度為155個點，長度較周期圖法里的信號短，信號為，其中，

8、n為0155/100，間隔為1/100。圖8為自相關(guān)函數(shù)的無偏估計。圖9為功率譜。從圖9可以看出，這種方法具有極高的分辨能力。只是在20和30處有兩個峰值，在其它地方的值為零。將信號改變成以下信號：則圖11變成下圖：輸出誤差功率 0.0356 輸出階數(shù)P 2 （F1=20，F(xiàn)2=22）輸出誤差功率 0.0485 輸出階數(shù)P 11 （F1=20，F(xiàn)2=25）輸出誤差功率 0.0403 輸出階數(shù)P 9 （F1=20，F(xiàn)2=23）試驗一下發(fā)現(xiàn)兩個峰值頻率為20和23的時候開始可以分辨出來。與經(jīng)典譜估計方法進行一個對比（基于第一個信號）：相比與前面幾種方法得到的結(jié)果來看，相差非常大，不僅僅是分辨率的提

9、高，其余無效噪音或者說擾動都是非常小的。結(jié)論通過實驗可以直接看出以下特性：1) 經(jīng)典功率譜估計的方差大，譜分辨率差。采用經(jīng)典的傅里葉變換及窗口截斷，對長序列有良好估計。2) 現(xiàn)代譜的分辨率較高。這是由于在時域的開窗，使得在頻域發(fā)生“頻譜泄漏”，即功率譜的主瓣能量泄漏到旁瓣中，導(dǎo)致弱信號的主瓣被強信號的旁瓣所湮沒，造成譜的模糊。 3) 平均周期圖法的收斂性較好，曲線平滑，但是功率譜主瓣較寬，分辨率低。這是由于對隨機序列的分段處理引起了長度有限所帶來的 Gibbs 現(xiàn)象而造成的。 參考文獻1. 劉志剛,李錄明,趙冬梅. 現(xiàn)代譜估計法及應(yīng)用效果J. 石油地球物理勘探,2009,S1:5-9+167+

10、9.2. 樊劍,劉鐵,胡亮. 基于現(xiàn)代時頻分析技術(shù)的地震波時變譜估計J. 振動與沖擊,2007,11:79-82+98+184.3. 蔡希玲,劉學(xué)偉,呂英梅,曹錫娜. 統(tǒng)計F-X譜估計方法及應(yīng)用J. 勘探地球物理進展,2008,03:181-186+163.4. 李剛. 寬帶信號空間譜估計算法研究D.哈爾濱工程大學(xué),2007.5. 姚武川,姚天任. 經(jīng)典譜估計方法的MATLAB分析J. 華中理工大學(xué)學(xué)報,2000,04:45-47.6. 王鳳瑛,張麗麗. 功率譜估計及其MATLAB仿真J. 微計算機信息,2006,31:287-289.7. 王福杰,潘宏俠. MATLAB中幾種功率譜估計函數(shù)的

11、比較分析與選擇J. 電子產(chǎn)品可靠性與環(huán)境試驗,2009,06:28-31.附錄：譜估計的Matlab實現(xiàn)程序1：經(jīng)典法譜估計clf;Fs=1000;N=256;Nfft=256;%數(shù)據(jù)的長度和FFT所用的數(shù)據(jù)長度n=0:N-1;t=n/Fs;%采用的時間序列xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);Pxx=10*log10(abs(fft(xn,Nfft).2)/N);%Fourier振幅譜平方的平均值，并轉(zhuǎn)化為dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%給出頻率序列subplot(2,1,1),plot(f,

12、Pxx);%繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜/dB');title('周期圖 N=256');grid on;程序2：經(jīng)典譜加窗分析 Fs=1000;%采樣頻率n=0:1/Fs:1;%產(chǎn)生含有噪音的信號xn=sin(2*pi*50*n)+2*sin(2*pi*120*n)+randn(size(n);plot(n,xn);window=boxcar(length(xn);%矩形窗nfft=512;Pxx,f=periodogram(xn,window,nfft,Fs);%直接法figure(1)plot(f,1

13、0*log10(Pxx);window=boxcar(length(xn);%矩形窗nfft=1024;Pxx,f=periodogram(xn,window,nfft,Fs);%直接法figure(2)plot(f,10*log10(Pxx); 程序3：Bartlett平均周期圖法fs=1000;n=0:1/fs:1;xn=sin(2*pi*50*n)+2*sin(2*pi*120*n)+randn(size(n);nfft=1024;window=hamming(nfft);noverlap=0;p=0.9;Pxx,Pxxc=psd(xn,nfft,fs,window,noverlap,p

14、);index=0:round(nfft/2-1);k=index*fs/nfft;plot_Pxx=10*log10(Pxx(index+1);plot_Pxxc=10*log10(Pxx(index+1);figure(1)plot(k,plot_Pxx)figure(2)plot(k,plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc);程序4：現(xiàn)代譜估計法clear all% 產(chǎn)生信號 FS=100; %設(shè)采樣頻率為100，則根據(jù)F1/FS=0.2，F(xiàn)2/FS=0.3 OR 0.25 ，可以得到F1，F(xiàn)2便于計算；N=155; %數(shù)據(jù)長度 改變

15、數(shù)據(jù)長度會導(dǎo)致分辨率的變化；n=0:1/FS:N/FS; F1=20; %第一個SIN信號的頻率；F2=30; %第二個SIN信號的頻率，取25或者30；X=sin(2*pi*F1*n)+sin(2*pi*F2*n)+(1/20)*randn(size(n) ; %產(chǎn)生信號，由信噪比為10dB推出噪聲功率；信號長度從X（1）到X（N+1）；% 產(chǎn)生自相關(guān)函數(shù)數(shù)組 for m=1:N+1 %初始化R（m），R（m）用來存放自相關(guān)函數(shù)；由于R（0）在MATLAB里無效，所以從1開始到N+1； R(m)=0;endfor m=1:N+1 %做自相關(guān)函數(shù)的無偏估計； S=0; for n=1:N+2-

16、m H=X(n)*X(n+m-1); S=S+H; end R(m)=S/N;end %估計完畢 subplot(3,1,1);plot(X);subplot(3,1,2);plot(R);% Levinson算法主體部分 for k=1:N+1 %初始化后面算法中要用到的數(shù)組，其中a(k,k)用來存放AR模型參數(shù)，F(xiàn)PE(k)是最終預(yù)測誤差準則函數(shù)， a(k,k)=0; % U2(k)是AR模型中的另一參數(shù)即誤差功率； FPE(k)=0; U2(k+1)=0;endU2(1)=R(1); a(1,1)=-R(2)/R(1); %由Levinson算法，計算一階模型參數(shù)a11；U2(2)=(1

17、-(abs(a(1,1)2)*U2(1); %由Levinson算法，計算一階模型參數(shù) 誤差功率；FPE(1)=U2(2)*(N+2)/N; %計算一階模型的最終預(yù)測誤差準則函數(shù)；S=0;for k=2:N %由Levinson梯歸公式計算K階模型參數(shù)和FPE函數(shù)； for l=1:k-1 M=a(k-1,l)*R(k-l+1); S=S+M; end a(k,k)=-(R(k+1)+S)/U2(k); for i=1:k-1 a(k,i)=a(k-1,i)+a(k,k)*a(k-1,k-i); end U2(k+1)=(1-(abs(a(k,k)2)*U2(k); FPE(k)=U2(k+1)*(N+k+1)/(N-k+1); S=0;end %梯歸計算完畢；% 確定階數(shù)P min=FPE(1); %求出使得FPE函數(shù)取最小值的階數(shù)P；for k=2:N if FPE(k)<min min=FPE(k); p=k; endend %p=2 %為了調(diào)整效果可以在這里自行指定階數(shù)；disp('輸出模型參數(shù)a');for k=1:p disp(a(p,k);endd