功率譜密度相關(guān)方法的MATLAB實現(xiàn)_第1頁
功率譜密度相關(guān)方法的MATLAB實現(xiàn)_第2頁
功率譜密度相關(guān)方法的MATLAB實現(xiàn)_第3頁
功率譜密度相關(guān)方法的MATLAB實現(xiàn)_第4頁
功率譜密度相關(guān)方法的MATLAB實現(xiàn)_第5頁
已閱讀5頁,還剩5頁未讀 繼續(xù)免費閱讀

下載本文檔

版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請進(jìn)行舉報或認(rèn)領(lǐng)

文檔簡介

1、1. 基本方法周期圖法是直接將信號的采樣數(shù)據(jù)x(n)進(jìn)行Fourier變換求取功率譜密度估計的方法。假定有限長隨機(jī)信號序列為x(n)。它的Fourier變換和功率譜密度估計存在下面的關(guān)系:                     式中,N為隨機(jī)信號序列x(n)的長度。在離散的頻率點f=kf,有:  其中,F(xiàn)FTx(n)為對序列x(n)的Fourier變換,由于FFTx(n)的周期為N,求得的功率譜估計以N為周期,因此這種方法稱為周期圖法。下面用例子說明如何采用這種方

2、法進(jìn)行功率譜用有限長樣本序列的Fourier變換來表示隨機(jī)序列的功率譜,只是一種估計或近似,不可避免存在誤差。為了減少誤差,使功率譜估計更加平滑,可采用分段平均周期圖法(Bartlett法)、加窗平均周期圖法(Welch法)等方法加以改進(jìn)。2. 分段平均周期圖法(Bartlett法)將信號序列x(n),n=0,1,N-1,分成互不重疊的P個小段,每小段由m個采樣值,則P*m=N。對每個小段信號序列進(jìn)行功率譜估計,然后再取平均作為整個序列x(n)的功率譜估計。平均周期圖法還可以對信號x(n)進(jìn)行重疊分段,如按2:1重疊分段,即前一段信號和后一段信號有一半是重疊的。對每一小段信號序列進(jìn)行

3、功率譜估計,然后再取平均值作為整個序列x(n)的功率譜估計。這兩種方法都稱為平均周期圖法,一般后者比前者好。程序運行結(jié)果為圖9-5,上圖采用不重疊分段法的功率譜估計,下圖為2:1重疊分段的功率譜估計,可見后者估計曲線較為平滑。與上例比較,平均周期圖法功率譜估計具有明顯效果(漲落曲線靠近0dB)。3.加窗平均周期圖法加窗平均周期圖法是對分段平均周期圖法的改進(jìn)。在信號序列x(n)分段后,用非矩形窗口對每一小段信號序列進(jìn)行預(yù)處理,再采用前述分段平均周期圖法進(jìn)行整個信號序列x(n)的功率譜估計。由窗函數(shù)的基本知識(第7章)可知,采用合適的非矩形窗口對信號進(jìn)行處理可減小“頻譜泄露”,同時可增加頻峰的寬度

4、,從而提高頻譜分辨率。其中上圖采用無重疊數(shù)據(jù)分段的加窗平均周期圖法進(jìn)行功率譜估計,而下圖采用重疊數(shù)據(jù)分段的加窗平均周期圖法進(jìn)行功率譜估計,顯然后者是更佳的,信號譜峰加寬,而噪聲譜均在0dB附近,更為平坦(注意采用無重疊數(shù)據(jù)分段噪聲的最大的下降分貝數(shù)大于5dB,而重疊數(shù)據(jù)分段周期圖法噪聲的最大下降分貝數(shù)小于5dB)。4.  Welch法估計及其MATLAB函數(shù)Welch功率譜密度就是用改進(jìn)的平均周期圖法來求取隨機(jī)信號的功率譜密度估計的。Welch法采用信號重疊分段、加窗函數(shù)和FFT算法等計算一個信號序列的自功率譜估計(PSD如上例中的下半部分的求法)和兩個信號序列的互功率譜估

5、計(CSD)。MATLAB信號處理工具箱函數(shù)提供了專門的函數(shù)PSD和CSD自動實現(xiàn)Welch法估計,而不需要自己編程。(1) 函數(shù)psd利用Welch法估計一個信號自功率譜密度,函數(shù)調(diào)用格式為:Pxx,f=psd(x,Nfft,Fs,window,Noverlap,dflag)式中,x為信號序列;Nfft為采用的FFT長度。這一值決定了功率譜估計速度,當(dāng)Nfft采用2的冪時,程序采用快速算法;Fs為采樣頻率;Window定義窗函數(shù)和x分段序列的長度。窗函數(shù)長度必須小于或等于Nfft,否則會給出錯誤信息;Noverlap為分段序列重疊的采樣點數(shù)(長度),它應(yīng)小于Nfft;dflag為去除信號趨勢

6、分量的選擇項:linear,去除線性趨勢分量,mean去除均值分量,none不做去除趨勢處理。Pxx為信號x的自功率譜密度估計。f為返回的頻率向量,它和Pxx對應(yīng),并且有相同長度。在psd函數(shù)調(diào)用格式中,缺省值為:Nfft=min(256,length(x),Fs=2Hz, window=hanning(Nfft),noverlap=0. 若x是實序列,函數(shù)psd僅計算頻率為正的功率注意程序前半部分中頻率向量f的創(chuàng)建方法。它與函數(shù)psd的輸出Pxx長度的關(guān)系如下:若x為實序列,當(dāng)Nfft為奇數(shù)時,f=(0:(Nfft+1)/2-1)/Nfft;當(dāng)Nfft為偶數(shù)時,f=(0:Nfft/2)/Nf

7、ft。函數(shù)還有一種缺省返回值的調(diào)用格式,用于直接繪制信號序列x的功率譜估計曲線。函數(shù)還可以計算帶有置信區(qū)間的功率譜估計,調(diào)用格式為:Pxx,Pxxc,f=psd(x,Nfft,Fs,window,Noverlap,p)式中,p為置信區(qū)間,0<=p<=1。 由此可知,濾波器輸入白噪聲序列的輸出信號的功率譜或自相關(guān)可以確定濾波器的頻率特性。(2)函數(shù)csd利用welch法估計兩個信號的互功率譜密度,函數(shù)調(diào)用格式為: Pxy,f=csd(x,y,Nfft,Fs,window,Noverlap,dflag) Pxy,Pxyc,f=csd(x,y,Nfft,Fs,window,No

8、verlap,p)這里,x,y為兩個信號序列;Pxy為x,y的互功率譜估計;其他參數(shù)的意義同自功率譜函數(shù)psd。可以看到,兩個白噪聲信號的互功率譜(上圖)雜亂無章,看不出周期成分,大部分功率譜在-5dB以下。然而白噪聲與帶有噪聲的周期信號的功率譜在其周期(頻率為1000Hz)處有一峰值,清楚地表明了周期信號的周期或頻率。因此,利用未知信號與白噪聲信號的互功率譜也可以檢測未知信號中所含有的頻率成分。 5 多 窗 口 法多窗口法(Multitaper method,簡稱MTM法)利用多個正交窗口(Tapers)獲得各自獨立的近似功率譜估計,然后綜合這些估計得到一個序列的功率譜估計。相對于

9、普通的周期圖法,這種功率譜估計具有更大的自由度,并在估計精度和估計波動方面均有較好的效果。普通的功率譜估計只利用單一窗口,因此在序列始端和末端均會丟失相關(guān)信息,而且無法找回。而MTM法估計增加窗口使得丟失的信息盡量減少。MTM法簡單地采用一個參數(shù):時間帶寬積(Time-bandwidth product)NW,這個參數(shù)用以定義計算功率譜所用窗的數(shù)目,為2*NW-1。NW越大,功率譜計算次數(shù)越多,時間域分辨率越高,而頻率域分辨率降低,使得功率譜估計的波動減小。隨著NW增大,每次估計中譜泄漏增多,總功率譜估計的偏差增大。對于每一個數(shù)據(jù)組,通常有一個最優(yōu)的NW使得在估計偏差和估計波動兩方面求得折中,

10、這需要在程序中反復(fù)調(diào)試來獲得。MATLAB信號處理工具箱中函數(shù)PMTM就是采用MTM法估計功率譜密度。函數(shù)調(diào)用格式為:Pxx,f=pmtm(x,nw,Nfft,Fs)式中,x為信號序列;nw為時間帶寬積,缺省值為4。通??扇?,5/2,3,7/2;Nfft為FFT長度;Fs為采樣頻率。上面的函數(shù)還可以通過無返回值而繪出置信區(qū)間,如pmtm(x,nw,Nfft,Fs,option,p)繪制帶置信區(qū)間的功率譜密度估計曲線,0<=p<=1。6 最 大 熵 法(Maxmum entropy method, MEM法)如上所述,周期圖法功率譜估計需要對信號序列“截斷”或加窗處理,其結(jié)果是使估

11、計的功率譜密度為信號序列真實譜和窗譜的卷積,導(dǎo)致誤差的產(chǎn)生。最大熵功率譜估計的目的是最大限度地保留截斷后丟失的“窗口”以外信號的信息,使估計譜的熵最大。主要方法是以已知的自相關(guān)序列rxx(0),rxx(1),rxx(p)為基礎(chǔ),外推自相關(guān)序列rxx(p+1),rxx(p+2),,保證信息熵最大。最大熵功率譜估計法假定隨機(jī)過程是平穩(wěn)高斯過程,可以證明,隨機(jī)信號的最大熵譜與AR自回歸(全極點濾波器)模型譜是等價的。MATLAB信號處理工具箱提供最大熵功率譜估計函數(shù)pmem,其調(diào)用格式為:Pxx,f,a=pmem(x,p,Nfft,Fs,xcorr)式中,x為輸入信號序列或輸入相關(guān)矩陣;p為全極點濾

12、波器階次;a為全極點濾波器模型系數(shù)向量;xcorr是把x認(rèn)為是相關(guān)矩陣。 比較最大熵功率譜估計(MEM)和改進(jìn)的平均周期圖功率譜估計,可見,MEM法估計的功率譜曲線較光滑。在這一方法中,MEM法選定全極點濾波器的階數(shù)取得越大,能夠獲得的窗口外的信息越多,但計算量也越大,需要根據(jù)情況折中考慮。 7 多信號分類法MATLAB信號處理工具箱還提供另一種功率譜估計函數(shù)pmusic。該函數(shù)執(zhí)行多信號分類法(multiple signal classification, Music法)。將數(shù)據(jù)自相關(guān)矩陣看成由信號自相關(guān)矩陣和噪聲自相關(guān)矩陣兩部分組成,即數(shù)據(jù)自相關(guān)矩陣R包含有兩個子空間信

13、息:信號子空間和噪聲子空間。這樣,矩陣特征值向量(Eigen vector)也可分為兩個子空間:信號子空間和噪聲子空間。為了求得功率譜估計,函數(shù)pmusic計算信號子空間和噪聲子空間的特征值向量函數(shù),使得在周期信號頻率處函數(shù)值最大,功率譜估計出現(xiàn)峰值,而在其他頻率處函數(shù)值最小。其調(diào)用格式為:Pxx,f,a=pmusic(x,p,thresh,Nfft,Fs,window,Noverlap)式中,x為輸入信號的向量或矩陣;p為信號子空間維數(shù);thresh為閾值,其他參數(shù)的意義與函數(shù)psd相同。 功率譜密度相關(guān)方法的MATLAB實現(xiàn)%分段平均周期圖法(Bartlett法)%運用信號不重疊

14、分段估計功率譜 Nsec=256;n=0:sigLength-1;t=n/Fs; %數(shù)據(jù)點數(shù),分段間隔,時間序列pxx1=abs(fft(y(1:256),Nsec).2)/Nsec; %第一段功率譜pxx2=abs(fft(y(257:512),Nsec).2)/Nsec; %第二段功率譜pxx3=abs(fft(y(515:768),Nsec).2)/Nsec; %第三段功率譜pxx4=abs(fft(y(769:1024),Nsec).2)/Nsec; %第四段功率譜Pxx=10*log10(pxx1+pxx2+pxx3+pxx4)/4); %平均得到整個序列功率譜f=(0:length

15、(Pxx)-1)*Fs/length(Pxx); %給出功率譜對應(yīng)的頻率%plot(f(1:Nsec/2),Pxx(1:Nsec/2); %繪制功率譜曲線figure,plot(f(1:Nsec),Pxx(1:Nsec); %繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜 /dB');title('平均周期圖(無重疊) N=4*256');grid on%運用信號重疊分段估計功率譜 pxx1=abs(fft(y(1:256),Nsec).2)/Nsec; %第一段功率譜pxx2=abs(fft(y(129:384),Nsec

16、).2)/Nsec; %第二段功率譜pxx3=abs(fft(y(257:512),Nsec).2)/Nsec; %第三段功率譜pxx4=abs(fft(y(385:640),Nsec).2)/Nsec; %第四段功率譜pxx5=abs(fft(y(513:768),Nsec).2)/Nsec; %第五段功率譜pxx6=abs(fft(y(641:896),Nsec).2)/Nsec; %第六段功率譜pxx7=abs(fft(y(769:1024),Nsec).2)/Nsec; %第七段功率譜Pxx=10*log10(pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7)

17、; %功率譜平均并轉(zhuǎn)化為dBf=(0:length(Pxx)-1)*Fs/length(Pxx); %頻率序列figure,plot(f(1:Nsec/2),Pxx(1:Nsec/2); %繪制功率譜曲線xlabel('頻率/Hz'); ylabel('功率譜/dB');title('平均周期圖(重疊一半) N=1024');grid on %Nsec=256;n=0:sigLength-1;t=n/Fs; %數(shù)據(jù)長度,分段數(shù)據(jù)長度、時間序列w=hanning(256); %采用的窗口數(shù)據(jù)%采用不重疊加窗方法的功率譜估計pxx1=abs(fft(

18、w.*y(1:256),Nsec).2)/norm(w)2; %第一段加窗振幅譜平方pxx2=abs(fft(w.*y(257:512),Nsec).2)/norm(w)2; %第二段加窗振幅譜平方pxx3=abs(fft(w.*y(513:768),Nsec).2)/norm(w)2; %第三段加窗振幅譜平方pxx4=abs(fft(w.*y(769:1024),Nsec).2)/norm(w)2; %第四段加窗振幅譜平方Pxx=10*log10(pxx1+pxx2+pxx3+pxx4)/4); %求得平均功率譜,轉(zhuǎn)換為dBf=(0:length(Pxx)-1)*Fs/length(Pxx)

19、; %求得頻率序列figuresubplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2); %繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜/dB');title('加窗平均周期圖(無重疊) N=4*256');grid on%采用重疊加窗方法的功率譜估計pxx1=abs(fft(w.*y(1:256),Nsec).2)/norm(w)2; %第一段加窗振幅譜平方pxx2=abs(fft(w.*y(129:384),Nsec).2)/norm(w)2; %第二段加窗振幅譜平方pxx3=abs(

20、fft(w.*y(257:512),Nsec).2)/norm(w)2; %第三段加窗振幅譜平方pxx4=abs(fft(w.*y(385:640),Nsec).2)/norm(w)2; %第四段加窗振幅譜平方pxx5=abs(fft(w.*y(513:768),Nsec).2)/norm(w)2; %第五段加窗振幅譜平方pxx6=abs(fft(w.*y(641:896),Nsec).2)/norm(w)2; %第六段加窗振幅譜平方pxx7=abs(fft(w.*y(769:1024),Nsec).2)/norm(w)2; %第七段加窗振幅譜平方Pxx=10*log10(pxx1+pxx2+

21、pxx3+pxx4+pxx5+pxx6+pxx7)/7);%平均功率譜轉(zhuǎn)換為dBf=(0:length(Pxx)-1)*Fs/length(Pxx); %頻率序列subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2); %繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜/dB');title('加窗平均周期圖(重疊一半)N=1024');grid on%4分段平均周期圖法(hanning窗)Nsec=256;n=0:sigLength-1;t=n/Fs;w=hanning(256);Pxx1

22、=abs(fft(w.*y(1:256),Nsec).2)/Nsec;Pxx2=abs(fft(w.*y(257:512),Nsec).2)/Nsec;Pxx3=abs(fft(w.*y(513:768),Nsec).2)/Nsec;Pxx4=abs(fft(w.*y(769:1024),Nsec).2)/Nsec;Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4)/4);f=(0:length(Pxx)-1)*Fs/length(Pxx);figuresubplot(2,1,1)plot(f,Pxx);xlabel('頻率/Hz');ylabel('功

23、率譜/dB');title('Averaged Modified Periodogram (none overlap) N=4*256');grid%4分段(2:1重疊)平均周期圖法(hanning窗)Nsec=256;n=0:sigLength-1;t=n/Fs;w=hanning(256);Pxx1=abs(fft(w.*y(1:256),Nsec).2)/Nsec;Pxx2=abs(fft(w.*y(129:384),Nsec).2)/Nsec;Pxx3=abs(fft(w.*y(257:512),Nsec).2)/Nsec;Pxx4=abs(fft(w.*y(3

24、85:640),Nsec).2)/Nsec;Pxx5=abs(fft(w.*y(513:768),Nsec).2)/Nsec;Pxx6=abs(fft(w.*y(641:896),Nsec).2)/Nsec;Pxx7=abs(fft(w.*y(769:1024),Nsec).2)/Nsec;Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);f=(0:length(Pxx)-1)*Fs/length(Pxx);subplot(2,1,2);plot(f,Pxx);xlabel('頻率/Hz');ylabel('Powe

25、r Spectrum (dB)');title('Averaged Modified Periodogram (half overlap) N=1024');grid%PSD_WELCH方法%采樣頻率Nfft=256;n=0:sigLength-1;t=n/Fs; %數(shù)據(jù)長度、時間序列window=hanning(256); %選用的窗口noverlap=128; %分段序列重疊的采樣點數(shù)(長度)dflag='none' %不做趨勢處理Pxx,Pxxc,f=psd(y,Nfft,Fs,window,noverlap,0.95); %功率譜估計,并以0.9

26、5的置信度給出置信區(qū)間,無返回值是繪制出置信區(qū)間figureplot(f,10*log10(Pxx); %繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');title('PSDWelch方法'); grid on  %最大熵法(MEM法)Nfft=256;n=0:sigLength-1;t=n/Fs; %數(shù)據(jù)長度、分段長度和時間序列window=hanning(256); %采用窗口Pxx1,f=pmem(x,20,Nfft,Fs); %采用最大熵法,采用濾波器階數(shù)14,估計功率譜figure,subplot(

27、2,1,1),plot(f,10*log10(Pxx1); %繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');title('最大熵法 Order=20原始信號功率譜');grid onPxx1,f=pmem(y0,20,Nfft,Fs); %采用最大熵法,采用濾波器階數(shù)14,估計功率譜subplot(2,1,2),plot(f,10*log10(Pxx1); %繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');axis(0 4000 -20 0)title(&#

28、39;最大熵法 Order=20濾波后的信號功率譜');grid on%功率譜密度%PSD_WELCH方法Nfft=512;n=0:sigLength-1;t=n/Fs; %數(shù)據(jù)長度、時間序列window=hanning(256); %選用的窗口noverlap=128; %分段序列重疊的采樣點數(shù)(長度)dflag='none' %不做趨勢處理Pxx,Pxxc,f=psd(x,Nfft,Fs,window,noverlap,0.95); %功率譜估計,并以0.95的置信度給出置信區(qū)間,無返回值是繪制出置信區(qū)間figure;subplot(211);plot(f,10*l

29、og10(Pxx); %繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');grid on;title('PSDWelch方法的原始信號功率譜')subplot(212)Pxx,Pxxc,f=psd(y0,Nfft,Fs,window,noverlap,0.95); %功率譜估計,并以0.95的置信度給出置信區(qū)間,無返回值是繪制出置信區(qū)間plot(f,10*log10(Pxx); %繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');axis(0 4000 -30 0)grid on;title('PSDWelch方法的濾波后的信號功率譜')%用多窗口法(MTM)n=0:sigLe

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會有圖紙預(yù)覽,若沒有圖紙預(yù)覽就沒有圖紙。
  • 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
  • 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護(hù)處理,對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對任何下載內(nèi)容負(fù)責(zé)。
  • 6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時也不承擔(dān)用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

評論

0/150

提交評論