功率譜密度機(jī)器實(shí)現(xiàn)

1、1. 基本方法周期圖法是直接將信號(hào)的采樣數(shù)據(jù)x(n)進(jìn)行Fourier變換求取功率譜密度估計(jì)的方法。假定有限長(zhǎng)隨機(jī)信號(hào)序列為x(n)。它的Fourier變換和功率譜密度估計(jì)存在下面的關(guān)系：式中，N為隨機(jī)信號(hào)序列x(n)的長(zhǎng)度。在離散的頻率點(diǎn)f=kf,有：其中，F(xiàn)FTx(n)為對(duì)序列x(n)的Fourier變換，由于FFTx(n)的周期為N，求得的功率譜估計(jì)以N為周期，因此這種方法稱為周期圖法。下面用例子說(shuō)明如何采用這種方法進(jìn)行功率譜用有限長(zhǎng)樣本序列的Fourier變換來(lái)表示隨機(jī)序列的功率譜，只是一種估計(jì)或近似，不可避免存在誤差。為了減少誤差，使功率譜估計(jì)更加平滑，可采用分段平均周期圖

2、法（Bartlett法）、加窗平均周期圖法（Welch法）等方法加以改進(jìn)。2. 分段平均周期圖法（Bartlett法）將信號(hào)序列x(n)，n=0,1,N-1，分成互不重疊的P個(gè)小段，每小段由m個(gè)采樣值，則P*m=N。對(duì)每個(gè)小段信號(hào)序列進(jìn)行功率譜估計(jì)，然后再取平均作為整個(gè)序列x(n)的功率譜估計(jì)。平均周期圖法還可以對(duì)信號(hào)x(n)進(jìn)行重疊分段，如按2:1重疊分段，即前一段信號(hào)和后一段信號(hào)有一半是重疊的。對(duì)每一小段信號(hào)序列進(jìn)行功率譜估計(jì)，然后再取平均值作為整個(gè)序列x(n)的功率譜估計(jì)。這兩種方法都稱為平均周期圖法，一般后者比前者好。程序運(yùn)行結(jié)果為圖9-5，上圖采用不重疊分段法的功率譜估計(jì)，

3、下圖為2:1重疊分段的功率譜估計(jì)，可見(jiàn)后者估計(jì)曲線較為平滑。與上例比較，平均周期圖法功率譜估計(jì)具有明顯效果（漲落曲線靠近0dB）。3.加窗平均周期圖法加窗平均周期圖法是對(duì)分段平均周期圖法的改進(jìn)。在信號(hào)序列x(n)分段后，用非矩形窗口對(duì)每一小段信號(hào)序列進(jìn)行預(yù)處理，再采用前述分段平均周期圖法進(jìn)行整個(gè)信號(hào)序列x(n)的功率譜估計(jì)。由窗函數(shù)的基本知識(shí)（第7章）可知，采用合適的非矩形窗口對(duì)信號(hào)進(jìn)行處理可減小“頻譜泄露”，同時(shí)可增加頻峰的寬度，從而提高頻譜分辨率。其中上圖采用無(wú)重疊數(shù)據(jù)分段的加窗平均周期圖法進(jìn)行功率譜估計(jì)，而下圖采用重疊數(shù)據(jù)分段的加窗平均周期圖法進(jìn)行功率譜估計(jì)，顯然后者是更佳的，信號(hào)譜峰加

4、寬，而噪聲譜均在0dB附近，更為平坦（注意采用無(wú)重疊數(shù)據(jù)分段噪聲的最大的下降分貝數(shù)大于5dB，而重疊數(shù)據(jù)分段周期圖法噪聲的最大下降分貝數(shù)小于5dB）。4.  Welch法估計(jì)及其MATLAB函數(shù)Welch功率譜密度就是用改進(jìn)的平均周期圖法來(lái)求取隨機(jī)信號(hào)的功率譜密度估計(jì)的。Welch法采用信號(hào)重疊分段、加窗函數(shù)和FFT算法等計(jì)算一個(gè)信號(hào)序列的自功率譜估計(jì)(PSD如上例中的下半部分的求法)和兩個(gè)信號(hào)序列的互功率譜估計(jì)（CSD）。MATLAB信號(hào)處理工具箱函數(shù)提供了專門(mén)的函數(shù)PSD和CSD自動(dòng)實(shí)現(xiàn)Welch法估計(jì)，而不需要自己編程。（1） 函數(shù)psd利用Welch法估計(jì)

5、一個(gè)信號(hào)自功率譜密度，函數(shù)調(diào)用格式為：Pxx,f=psd(x,Nfft,Fs,window,Noverlap,dflag)式中，x為信號(hào)序列；Nfft為采用的FFT長(zhǎng)度。這一值決定了功率譜估計(jì)速度，當(dāng)Nfft采用2的冪時(shí)，程序采用快速算法；Fs為采樣頻率；Window定義窗函數(shù)和x分段序列的長(zhǎng)度。窗函數(shù)長(zhǎng)度必須小于或等于Nfft，否則會(huì)給出錯(cuò)誤信息；Noverlap為分段序列重疊的采樣點(diǎn)數(shù)（長(zhǎng)度），它應(yīng)小于Nfft;dflag為去除信號(hào)趨勢(shì)分量的選擇項(xiàng)：linear,去除線性趨勢(shì)分量，mean去除均值分量，none不做去除趨勢(shì)處理。Pxx為信號(hào)x的自功率譜密度估計(jì)。f為返回的頻率向量，它和Px

6、x對(duì)應(yīng)，并且有相同長(zhǎng)度。在psd函數(shù)調(diào)用格式中，缺省值為：Nfft=min(256,length(x),Fs=2Hz, window=hanning(Nfft),noverlap=0. 若x是實(shí)序列，函數(shù)psd僅計(jì)算頻率為正的功率注意程序前半部分中頻率向量f的創(chuàng)建方法。它與函數(shù)psd的輸出Pxx長(zhǎng)度的關(guān)系如下：若x為實(shí)序列，當(dāng)Nfft為奇數(shù)時(shí)，f=(0:(Nfft+1)/2-1)/Nfft;當(dāng)Nfft為偶數(shù)時(shí)，f=(0:Nfft/2)/Nfft。函數(shù)還有一種缺省返回值的調(diào)用格式，用于直接繪制信號(hào)序列x的功率譜估計(jì)曲線。函數(shù)還可以計(jì)算帶有置信區(qū)間的功率譜估計(jì)，調(diào)用格式為：Px

7、x,Pxxc,f=psd(x,Nfft,Fs,window,Noverlap,p)式中，p為置信區(qū)間，0<=p<=1。 由此可知，濾波器輸入白噪聲序列的輸出信號(hào)的功率譜或自相關(guān)可以確定濾波器的頻率特性。（2）函數(shù)csd利用welch法估計(jì)兩個(gè)信號(hào)的互功率譜密度，函數(shù)調(diào)用格式為：    Pxy,f=csd(x,y,Nfft,Fs,window,Noverlap,dflag)    Pxy,Pxyc,f=csd(x,y,Nfft,Fs,window,Noverlap,p)這里，x,y為兩個(gè)信號(hào)序列

8、；Pxy為x,y的互功率譜估計(jì)；其他參數(shù)的意義同自功率譜函數(shù)psd?？梢钥吹剑瑑蓚€(gè)白噪聲信號(hào)的互功率譜（上圖）雜亂無(wú)章，看不出周期成分，大部分功率譜在-5dB以下。然而白噪聲與帶有噪聲的周期信號(hào)的功率譜在其周期（頻率為1000Hz）處有一峰值，清楚地表明了周期信號(hào)的周期或頻率。因此，利用未知信號(hào)與白噪聲信號(hào)的互功率譜也可以檢測(cè)未知信號(hào)中所含有的頻率成分。 5 多 窗 口 法多窗口法(Multitaper method,簡(jiǎn)稱MTM法)利用多個(gè)正交窗口（Tapers）獲得各自獨(dú)立的近似功率譜估計(jì)，然后綜合這些估計(jì)得到一個(gè)序列的功率譜估計(jì)。相對(duì)于普通

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

10、折中，這需要在程序中反復(fù)調(diào)試來(lái)獲得。MATLAB信號(hào)處理工具箱中函數(shù)PMTM就是采用MTM法估計(jì)功率譜密度。函數(shù)調(diào)用格式為：Pxx,f=pmtm(x,nw,Nfft,Fs)式中，x為信號(hào)序列；nw為時(shí)間帶寬積，缺省值為4。通常可取2，5/2，3，7/2；Nfft為FFT長(zhǎng)度；Fs為采樣頻率。上面的函數(shù)還可以通過(guò)無(wú)返回值而繪出置信區(qū)間，如pmtm(x,nw,Nfft,Fs,option,p)繪制帶置信區(qū)間的功率譜密度估計(jì)曲線，0<=p<=1。6 最 大 熵 法（Maxmum entropy method， MEM法）如上所述

11、，周期圖法功率譜估計(jì)需要對(duì)信號(hào)序列“截?cái)唷被蚣哟疤幚?，其結(jié)果是使估計(jì)的功率譜密度為信號(hào)序列真實(shí)譜和窗譜的卷積，導(dǎo)致誤差的產(chǎn)生。最大熵功率譜估計(jì)的目的是最大限度地保留截?cái)嗪髞G失的“窗口”以外信號(hào)的信息，使估計(jì)譜的熵最大。主要方法是以已知的自相關(guān)序列rxx(0),rxx(1),rxx(p)為基礎(chǔ)，外推自相關(guān)序列rxx(p+1),rxx(p+2),，保證信息熵最大。最大熵功率譜估計(jì)法假定隨機(jī)過(guò)程是平穩(wěn)高斯過(guò)程，可以證明，隨機(jī)信號(hào)的最大熵譜與AR自回歸（全極點(diǎn)濾波器）模型譜是等價(jià)的。MATLAB信號(hào)處理工具箱提供最大熵功率譜估計(jì)函數(shù)pmem，其調(diào)用格式為：Pxx,f,a=pmem(x,p,Nfft,F

12、s,xcorr)式中，x為輸入信號(hào)序列或輸入相關(guān)矩陣；p為全極點(diǎn)濾波器階次；a為全極點(diǎn)濾波器模型系數(shù)向量；xcorr是把x認(rèn)為是相關(guān)矩陣。 比較最大熵功率譜估計(jì)(MEM)和改進(jìn)的平均周期圖功率譜估計(jì)，可見(jiàn)，MEM法估計(jì)的功率譜曲線較光滑。在這一方法中，MEM法選定全極點(diǎn)濾波器的階數(shù)取得越大，能夠獲得的窗口外的信息越多，但計(jì)算量也越大，需要根據(jù)情況折中考慮。 7 多信號(hào)分類法MATLAB信號(hào)處理工具箱還提供另一種功率譜估計(jì)函數(shù)pmusic。該函數(shù)執(zhí)行多信號(hào)分類法（multiple signal classification， Music法)。將數(shù)

13、據(jù)自相關(guān)矩陣看成由信號(hào)自相關(guān)矩陣和噪聲自相關(guān)矩陣兩部分組成，即數(shù)據(jù)自相關(guān)矩陣R包含有兩個(gè)子空間信息：信號(hào)子空間和噪聲子空間。這樣，矩陣特征值向量（Eigen vector）也可分為兩個(gè)子空間：信號(hào)子空間和噪聲子空間。為了求得功率譜估計(jì)，函數(shù)pmusic計(jì)算信號(hào)子空間和噪聲子空間的特征值向量函數(shù)，使得在周期信號(hào)頻率處函數(shù)值最大，功率譜估計(jì)出現(xiàn)峰值，而在其他頻率處函數(shù)值最小。其調(diào)用格式為：Pxx,f,a=pmusic(x,p,thresh,Nfft,Fs,window,Noverlap)式中,x為輸入信號(hào)的向量或矩陣；p為信號(hào)子空間維數(shù)；thresh為閾值，其他參數(shù)的意義與函數(shù)psd相同

14、。 功率譜密度相關(guān)方法的MATLAB實(shí)現(xiàn)%分段平均周期圖法（Bartlett法）%運(yùn)用信號(hào)不重疊分段估計(jì)功率譜 Nsec=256;n=0:sigLength-1;t=n/Fs;  %數(shù)據(jù)點(diǎn)數(shù)，分段間隔，時(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

15、(y(769:1024),Nsec).2)/Nsec; %第四段功率譜Pxx=10*log10(pxx1+pxx2+pxx3+pxx4)/4);   %平均得到整個(gè)序列功率譜f=(0:length(Pxx)-1)*Fs/length(Pxx);     %給出功率譜對(duì)應(yīng)的頻率%plot(f(1:Nsec/2),Pxx(1:Nsec/2);       %繪制功率譜曲線figure,plot(f(1:Nsec),Pxx(1:Nsec);&

16、#160;      %繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜 /dB');title('平均周期圖(無(wú)重疊) N=4*256');grid on%運(yùn)用信號(hào)重疊分段估計(jì)功率譜 pxx1=abs(fft(y(1:256),Nsec).2)/Nsec;   %第一段功率譜pxx2=abs(fft(y(129:384),Nsec).2)/Nsec; %第二段功率譜pxx3=abs(ff

17、t(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);

18、60;%功率譜平均并轉(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ù)長(zhǎng)度，分段數(shù)據(jù)

19、長(zhǎng)度、時(shí)間序列w=hanning(256);   %采用的窗口數(shù)據(jù)%采用不重疊加窗方法的功率譜估計(jì)pxx1=abs(fft(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; 

20、%第四段加窗振幅譜平方Pxx=10*log10(pxx1+pxx2+pxx3+pxx4)/4);    %求得平均功率譜，轉(zhuǎn)換為dBf=(0:length(Pxx)-1)*Fs/length(Pxx);  %求得頻率序列figuresubplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2); %繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜/dB');title('加窗平均周期圖(無(wú)重疊) N=4*256');gr

21、id on%采用重疊加窗方法的功率譜估計(jì)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(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),Ns

22、ec).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+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),P

23、xx(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=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),N

24、sec).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('功率譜/dB');title('Averaged Modified Periodogram (none overlap) N=4*256'

25、;);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(385:640),Nsec).2)/Nsec;Pxx5=abs(fft(w.*y(513:768),Nsec).2)/Nsec;Pxx6=abs(fft(w.*

26、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('Power Spectrum (dB)');title('Averaged Modified Periodog

27、ram (half overlap) N=1024');grid%PSD_WELCH方法%采樣頻率Nfft=256;n=0:sigLength-1;t=n/Fs; %數(shù)據(jù)長(zhǎng)度、時(shí)間序列window=hanning(256);   %選用的窗口noverlap=128;   %分段序列重疊的采樣點(diǎn)數(shù)（長(zhǎng)度）dflag='none'  %不做趨勢(shì)處理Pxx,Pxxc,f=psd(y,Nfft,Fs,window,noverlap,0.95); &

28、#160; %功率譜估計(jì),并以0.95的置信度給出置信區(qū)間，無(wú)返回值是繪制出置信區(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ù)長(zhǎng)度、分段長(zhǎng)度和時(shí)間序列window=

29、hanning(256);   %采用窗口Pxx1,f=pmem(x,20,Nfft,Fs);   %采用最大熵法，采用濾波器階數(shù)14，估計(jì)功率譜figure,subplot(2,1,1),plot(f,10*log10(Pxx1);   %繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');title('最大熵法 Order=20原始信號(hào)功率譜');grid onPxx1,f=pmem(y0,20,Nfft,Fs

30、);   %采用最大熵法，采用濾波器階數(shù)14，估計(jì)功率譜subplot(2,1,2),plot(f,10*log10(Pxx1);   %繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');axis(0 4000 -20 0)title('最大熵法 Order=20濾波后的信號(hào)功率譜');grid on%功率譜密度%PSD_WELCH方法Nfft=512;n=0:sigLength-1;t=n/Fs; %

31、數(shù)據(jù)長(zhǎng)度、時(shí)間序列window=hanning(256);   %選用的窗口noverlap=128;   %分段序列重疊的采樣點(diǎn)數(shù)（長(zhǎng)度）dflag='none'  %不做趨勢(shì)處理Pxx,Pxxc,f=psd(x,Nfft,Fs,window,noverlap,0.95);   %功率譜估計(jì),并以0.95的置信度給出置信區(qū)間，無(wú)返回值是繪制出置信區(qū)間figure;subplot(211);plot(f,10*log10(Pxx);  %繪制功率譜x

32、label('頻率/Hz');ylabel('功率譜/dB');grid on;title('PSDWelch方法的原始信號(hào)功率譜')subplot(212)Pxx,Pxxc,f=psd(y0,Nfft,Fs,window,noverlap,0.95);   %功率譜估計(jì),并以0.95的置信度給出置信區(qū)間，無(wú)返回值是繪制出置信區(qū)間plot(f,10*log10(Pxx);  %繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');axis(0 4000 -30 0)grid on;title('PSDWelch方法的濾波后的信號(hào)功率譜')%用多窗口法(MTM)n=0:sigLength-1;t=n/Fs;  %數(shù)據(jù)長(zhǎng)度、分段數(shù)據(jù)長(zhǎng)度，時(shí)間序列Pxx1,f=pmtm(x,2,Nfft,F