功率譜估計(jì)淺談

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

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

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

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

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

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

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

8、n為0155/100，間隔為1/100。圖8為自相關(guān)函數(shù)的無(wú)偏估計(jì)。圖9為功率譜。從圖9可以看出，這種方法具有極高的分辨能力。只是在20和30處有兩個(gè)峰值，在其它地方的值為零。將信號(hào)改變成以下信號(hào)：則圖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）試驗(yàn)一下發(fā)現(xiàn)兩個(gè)峰值頻率為20和23的時(shí)候開(kāi)始可以分辨出來(lái)。與經(jīng)典譜估計(jì)方法進(jìn)行一個(gè)對(duì)比（基于第一個(gè)信號(hào)）：相比與前面幾種方法得到的結(jié)果來(lái)看，相差非常大，不僅僅是分辨率的提

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

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

11、比較分析與選擇J. 電子產(chǎn)品可靠性與環(huán)境試驗(yàn),2009,06:28-31.附錄：譜估計(jì)的Matlab實(shí)現(xiàn)程序1：經(jīng)典法譜估計(jì)clf;Fs=1000;N=256;Nfft=256;%數(shù)據(jù)的長(zhǎng)度和FFT所用的數(shù)據(jù)長(zhǎng)度n=0:N-1;t=n/Fs;%采用的時(shí)間序列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);%繪制功率譜曲線(xiàn)xlabel('頻率/Hz');ylabel('功率譜/dB');title('周期圖 N=256');grid on;程序2：經(jīng)典譜加窗分析 Fs=1000;%采樣頻率n=0:1/Fs:1;%產(chǎn)生含有噪音的信號(hào)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)代譜估計(jì)法clear all% 產(chǎn)生信號(hào) FS=100; %設(shè)采樣頻率為100，則根據(jù)F1/FS=0.2，F(xiàn)2/FS=0.3 OR 0.25 ，可以得到F1，F(xiàn)2便于計(jì)算；N=155; %數(shù)據(jù)長(zhǎng)度 改變

15、數(shù)據(jù)長(zhǎng)度會(huì)導(dǎo)致分辨率的變化；n=0:1/FS:N/FS; F1=20; %第一個(gè)SIN信號(hào)的頻率；F2=30; %第二個(gè)SIN信號(hào)的頻率，取25或者30；X=sin(2*pi*F1*n)+sin(2*pi*F2*n)+(1/20)*randn(size(n) ; %產(chǎn)生信號(hào)，由信噪比為10dB推出噪聲功率；信號(hào)長(zhǎng)度從X（1）到X（N+1）；% 產(chǎn)生自相關(guān)函數(shù)數(shù)組 for m=1:N+1 %初始化R（m），R（m）用來(lái)存放自相關(guān)函數(shù)；由于R（0）在MATLAB里無(wú)效，所以從1開(kāi)始到N+1； R(m)=0;endfor m=1:N+1 %做自相關(guān)函數(shù)的無(wú)偏估計(jì)； 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 %估計(jì)完畢 subplot(3,1,1);plot(X);subplot(3,1,2);plot(R);% Levinson算法主體部分 for k=1:N+1 %初始化后面算法中要用到的數(shù)組，其中a(k,k)用來(lái)存放AR模型參數(shù)，F(xiàn)PE(k)是最終預(yù)測(cè)誤差準(zhǔn)則函數(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算法，計(jì)算一階模型參數(shù)a11；U2(2)=(1

17、-(abs(a(1,1)2)*U2(1); %由Levinson算法，計(jì)算一階模型參數(shù) 誤差功率；FPE(1)=U2(2)*(N+2)/N; %計(jì)算一階模型的最終預(yù)測(cè)誤差準(zhǔn)則函數(shù)；S=0;for k=2:N %由Levinson梯歸公式計(jì)算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 %梯歸計(jì)算完畢；% 確定階數(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