1.?基本方法

周期圖法是直接將信號的采樣數據x(n)進行Fourier變換求取功率譜密度估計的方法。假定有限長隨機信號序列為x(n)。它的Fourier變換和功率譜密度估計存在下面的關系：

??????????????????????????

式中，N為隨機信號序列x(n)的長度。在離散的頻率點f=kΔf,有：

??????????

其中，FFT[x(n)]為對序列x(n)的Fourier變換，由于FFT[x(n)]的周期為N，求得的功率譜估計以N為周期，因此這種方法稱為周期圖法。下面用例子說明如何采用這種方法進行功率譜

用有限長樣本序列的Fourier變換來表示隨機序列的功率譜，只是一種估計或近似，不可避免存在誤差。為了減少誤差，使功率譜估計更加平滑，可采用分段平均周期圖法（Bartlett法）、加窗平均周期圖法（Welch法）等方法加以改進。

2.?分段平均周期圖法（Bartlett法）

將信號序列x(n)，n=0,1,…,N-1，分成互不重疊的P個小段，每小段由m個采樣值，則P*m=N。對每個小段信號序列進行功率譜估計，然后再取平均作為整個序列x(n)的功率譜估計。

平均周期圖法還可以對信號x(n)進行重疊分段，如按2:1重疊分段，即前一段信號和后一段信號有一半是重疊的。對每一小段信號序列進行功率譜估計，然后再取平均值作為整個序列x(n)的功率譜估計。這兩種方法都稱為平均周期圖法，一般后者比前者好。程序運行結果為圖9-5，上圖采用不重疊分段法的功率譜估計，下圖為2:1重疊分段的功率譜估計，可見后者估計曲線較為平滑。與上例比較，平均周期圖法功率譜估計具有明顯效果（漲落曲線靠近0dB）。

3.加窗平均周期圖法

加窗平均周期圖法是對分段平均周期圖法的改進。在信號序列x(n)分段后，用非矩形窗口對每一小段信號序列進行預處理，再采用前述分段平均周期圖法進行整個信號序列x(n)的功率譜估計。由窗函數的基本知識（第7章）可知，采用合適的非矩形窗口對信號進行處理可減小“頻譜泄露”，同時可增加頻峰的寬度，從而提高頻譜分辨率。

其中上圖采用無重疊數據分段的加窗平均周期圖法進行功率譜估計，而下圖采用重疊數據分段的加窗平均周期圖法進行功率譜估計，顯然后者是更佳的，信號譜峰加寬，而噪聲譜均在0dB附近，更為平坦（注意采用無重疊數據分段噪聲的最大的下降分貝數大于5dB，而重疊數據分段周期圖法噪聲的最大下降分貝數小于5dB）。

4.??Welch法估計及其MATLAB函數

Welch功率譜密度就是用改進的平均周期圖法來求取隨機信號的功率譜密度估計的。Welch法采用信號重疊分段、加窗函數和FFT算法等計算一個信號序列的自功率譜估計(PSD如上例中的下半部分的求法)和兩個信號序列的互功率譜估計（CSD）。

MATLAB信號處理工具箱函數提供了專門的函數PSD和CSD自動實現Welch法估計，而不需要自己編程。

（1）?函數psd利用Welch法估計一個信號自功率譜密度，函數調用格式為：

[Pxx[,f]]=psd(x[,Nfft,Fs,window,Noverlap,’dflag’])

式中，x為信號序列；Nfft為采用的FFT長度。這一值決定了功率譜估計速度，當Nfft采用2的冪時，程序采用快速算法；Fs為采樣頻率；Window定義窗函數和x分段序列的長度。窗函數長度必須小于或等于Nfft，否則會給出錯誤信息；Noverlap為分段序列重疊的采樣點數（長度），它應小于Nfft;dflag為去除信號趨勢分量的選擇項：’linear’,去除線性趨勢分量，’mean’去除均值分量，’none’不做去除趨勢處理。Pxx為信號x的自功率譜密度估計。f為返回的頻率向量，它和Pxx對應，并且有相同長度。

在psd函數調用格式中，缺省值為：Nfft=min(256,length(x)),Fs=2Hz,?window=hanning(Nfft),noverlap=0.?若x是實序列，函數psd僅計算頻率為正的功率

注意程序前半部分中頻率向量f的創建方法。它與函數psd的輸出Pxx長度的關系如下：若x為實序列，當Nfft為奇數時，f=(0:(Nfft+1)/2-1)/Nfft;當Nfft為偶數時，f=(0:Nfft/2)/Nfft。

函數還有一種缺省返回值的調用格式，用于直接繪制信號序列x的功率譜估計曲線。

函數還可以計算帶有置信區間的功率譜估計，調用格式為：

[Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,Noverlap,p)

式中，p為置信區間，0<=p<=1。

?

由此可知，濾波器輸入白噪聲序列的輸出信號的功率譜或自相關可以確定濾波器的頻率特性。

（2）函數csd利用welch法估計兩個信號的互功率譜密度，函數調用格式為：

????[Pxy[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,’dflag’)

????[Pxy,Pxyc[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,p)

這里，x,y為兩個信號序列；Pxy為x,y的互功率譜估計；其他參數的意義同自功率譜函數psd。

可以看到，兩個白噪聲信號的互功率譜（上圖）雜亂無章，看不出周期成分，大部分功率譜在-5dB以下。然而白噪聲與帶有噪聲的周期信號的功率譜在其周期（頻率為1000Hz）處有一峰值，清楚地表明了周期信號的周期或頻率。因此，利用未知信號與白噪聲信號的互功率譜也可以檢測未知信號中所含有的頻率成分。

?

5 多?窗?口?法

多窗口法(Multitaper?method,簡稱MTM法)利用多個正交窗口（Tapers）獲得各自獨立的近似功率譜估計，然后綜合這些估計得到一個序列的功率譜估計。相對于普通的周期圖法，這種功率譜估計具有更大的自由度，并在估計精度和估計波動方面均有較好的效果。普通的功率譜估計只利用單一窗口，因此在序列始端和末端均會丟失相關信息，而且無法找回。而MTM法估計增加窗口使得丟失的信息盡量減少。

MTM法簡單地采用一個參數：時間帶寬積(Time-bandwidth?product)NW，這個參數用以定義計算功率譜所用窗的數目，為2*NW-1。NW越大，功率譜計算次數越多，時間域分辨率越高，而頻率域分辨率降低，使得功率譜估計的波動減小。隨著NW增大，每次估計中譜泄漏增多，總功率譜估計的偏差增大。對于每一個數據組，通常有一個最優的NW使得在估計偏差和估計波動兩方面求得折中，這需要在程序中反復調試來獲得。

MATLAB信號處理工具箱中函數PMTM就是采用MTM法估計功率譜密度。函數調用格式為：

[Pxx[,f]]=pmtm(x[,nw,Nfft,Fs])

式中，x為信號序列；nw為時間帶寬積，缺省值為4。通常可取2，5/2，3，7/2；Nfft為FFT長度；Fs為采樣頻率。

上面的函數還可以通過無返回值而繪出置信區間，如pmtm(x,nw,Nfft,Fs,’option’,p)繪制帶置信區間的功率譜密度估計曲線，0<=p<=1。

6 最?大?熵?法（Maxmum?entropy?method，?MEM法）

如上所述，周期圖法功率譜估計需要對信號序列“截斷”或加窗處理，其結果是使估計的功率譜密度為信號序列真實譜和窗譜的卷積，導致誤差的產生。

最大熵功率譜估計的目的是最大限度地保留截斷后丟失的“窗口”以外信號的信息，使估計譜的熵最大。主要方法是以已知的自相關序列rxx(0),rxx(1),…,rxx(p)為基礎，外推自相關序列rxx(p+1),rxx(p+2),…，保證信息熵最大。

最大熵功率譜估計法假定隨機過程是平穩高斯過程，可以證明，隨機信號的最大熵譜與AR自回歸（全極點濾波器）模型譜是等價的。

MATLAB信號處理工具箱提供最大熵功率譜估計函數pmem，其調用格式為：

[Pxx,f,a]=pmem(x,p,Nfft,Fs,’xcorr’)

式中，x為輸入信號序列或輸入相關矩陣；p為全極點濾波器階次；a為全極點濾波器模型系數向量；’xcorr’是把x認為是相關矩陣。

?

比較最大熵功率譜估計(MEM)和改進的平均周期圖功率譜估計，可見，MEM法估計的功率譜曲線較光滑。在這一方法中，MEM法選定全極點濾波器的階數取得越大，能夠獲得的窗口外的信息越多，但計算量也越大，需要根據情況折中考慮。

?

7 多信號分類法

MATLAB信號處理工具箱還提供另一種功率譜估計函數pmusic。該函數執行多信號分類法（multiple?signal?classification，?Music法)。將數據自相關矩陣看成由信號自相關矩陣和噪聲自相關矩陣兩部分組成，即數據自相關矩陣R包含有兩個子空間信息：信號子空間和噪聲子空間。這樣，矩陣特征值向量（Eigen?vector）也可分為兩個子空間：信號子空間和噪聲子空間。為了求得功率譜估計，函數pmusic計算信號子空間和噪聲子空間的特征值向量函數，使得在周期信號頻率處函數值最大，功率譜估計出現峰值，而在其他頻率處函數值最小。其調用格式為：

[Pxx,f,a]=pmusic(x,p[[,thresh],Nfft,Fs,window,Noverlap])

式中,x為輸入信號的向量或矩陣；p為信號子空間維數；thresh為閾值，其他參數的意義與函數psd相同。

?

功率譜密度相關方法的MATLAB實現

%% %分段平均周期圖法（Bartlett法） %運用信號不重疊分段估計功率譜 Nsec=256;n=0:sigLength-1;t=n/Fs; %數據點數，分段間隔，時間序列 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(Pxx)-1)*Fs/length(Pxx); %給出功率譜對應的頻率 %%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).^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); %功率譜平均并轉化為dB f=(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; %數據長度，分段數據長度、時間序列 w=hanning(256); %采用的窗口數據 %采用不重疊加窗方法的功率譜估計 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; %第四段加窗振幅譜平方 Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4); %求得平均功率譜，轉換為dB f=(0:length(Pxx)-1)*Fs/length(Pxx); %求得頻率序列 figure subplot(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(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+pxx3+pxx4+pxx5+pxx6+pxx7)/7);%平均功率譜轉換為dB f=(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=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); figure subplot(2,1,1) plot(f,Pxx); xlabel('頻率/Hz');ylabel('功率譜/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(385: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('Power Spectrum (dB)'); title('Averaged Modified Periodogram (half overlap) N=1024'); grid %% %PSD_WELCH方法 %采樣頻率 Nfft=256;n=0:sigLength-1;t=n/Fs; %數據長度、時間序列 window=hanning(256); %選用的窗口 noverlap=128; %分段序列重疊的采樣點數（長度） dflag='none'; %不做趨勢處理 [Pxx,Pxxc,f]=psd(y,Nfft,Fs,window,noverlap,0.95); %功率譜估計,并以0.95的置信度給出置信區間，無返回值是繪制出置信區間 figure plot(f,10*log10(Pxx)); %繪制功率譜 xlabel('頻率/Hz');ylabel('功率譜/dB'); title('PSD—Welch方法'); grid on%% %最大熵法（MEM法） Nfft=256;n=0:sigLength-1;t=n/Fs; %數據長度、分段長度和時間序列 window=hanning(256); %采用窗口 [Pxx1,f]=pmem(x,20,Nfft,Fs); %采用最大熵法，采用濾波器階數14，估計功率譜 figure,subplot(2,1,1),plot(f,10*log10(Pxx1)); %繪制功率譜 xlabel('頻率/Hz');ylabel('功率譜/dB'); title('最大熵法 Order=20原始信號功率譜');grid on [Pxx1,f]=pmem(y0,20,Nfft,Fs); %采用最大熵法，采用濾波器階數14，估計功率譜 subplot(2,1,2),plot(f,10*log10(Pxx1)); %繪制功率譜 xlabel('頻率/Hz');ylabel('功率譜/dB');axis([0 4000 -20 0]) title('最大熵法 Order=20濾波后的信號功率譜');grid on %% %%功率譜密度 %PSD_WELCH方法 Nfft=512;n=0:sigLength-1;t=n/Fs; %數據長度、時間序列 window=hanning(256); %選用的窗口 noverlap=128; %分段序列重疊的采樣點數（長度） dflag='none'; %不做趨勢處理 [Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,noverlap,0.95); %功率譜估計,并以0.95的置信度給出置信區間，無返回值是繪制出置信區間 figure;subplot(211); plot(f,10*log10(Pxx)); %繪制功率譜 xlabel('頻率/Hz');ylabel('功率譜/dB'); grid on;title('PSD—Welch方法的原始信號功率譜') subplot(212) [Pxx,Pxxc,f]=psd(y0,Nfft,Fs,window,noverlap,0.95); %功率譜估計,并以0.95的置信度給出置信區間，無返回值是繪制出置信區間 plot(f,10*log10(Pxx)); %繪制功率譜 xlabel('頻率/Hz');ylabel('功率譜/dB');axis([0 4000 -30 0]) grid on;title('PSD—Welch方法的濾波后的信號功率譜') %% %用多窗口法(MTM) n=0:sigLength-1;t=n/Fs; %數據長度、分段數據長度，時間序列 [Pxx1,f]=pmtm(x,2,Nfft,Fs); %用多窗口法(NW=4)估計功率譜 figure;subplot(2,1,1),plot(f,10*log10(Pxx1)); %繪制功率譜 xlabel('頻率/Hz');ylabel('功率譜/dB'); title('多窗口法(MTM) nw=2原始信號功率譜');grid on [Pxx,f]=pmtm(y0,2,Nfft,Fs); %用多窗口法(NW=2)估計功率譜 subplot(2,1,2),plot(f,10*log10(Pxx)); %繪制功率譜 xlabel('頻率/Hz');ylabel('功率譜/dB');axis([0 4000 -80 -20]) title('多窗口法(MTM) nw=2濾波后的信號功率譜');grid on%% %采用Welch方法估計功率譜 noverlap=128; %重疊數據 dflag='none'; Nfft=1024; figure;subplot(2,1,1) psd(x,Nfft,Fs,window,noverlap,dflag); %采用Welch方法估計功率譜 xlabel('頻率/Hz');ylabel('功率譜/dB') title('Welch方法原始信號功率譜');grid on subplot(2,1,2) psd(y0,Nfft,Fs,window,noverlap,dflag); %采用Welch方法估計功率譜 xlabel('頻率/Hz');ylabel('功率譜/dB'),axis([0 4000 -30 0]) title('Welch方法濾波后的信號功率譜');grid on pmusic(y0,[7,1.1],Nfft,Fs,32,16); %采用多信號分類法估計功率譜 xlabel('頻率/Hz'); ylabel('功率譜/dB') title('通過MUSIC法估計的偽譜')

總結

以上是生活随笔為你收集整理的功率谱密度相关方法MATLAB实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。