1、Matlab 信號處理工具箱 譜估計專題 頻譜分析Spectral estimation（譜估計）的目標是基于一個有限的數據集合描述一個信號的功率（在頻率上的）分布。功率譜估計在很多場合下都是有用的，包括對寬帶噪聲湮沒下的信號的檢測。從數學上看，一個平穩隨機過程的power spectrum（功率譜）和correlation sequence（相關序列）通過discrete-time Fourier transform（離散時間傅立葉變換）構成聯系。從normalized frequency（歸一化角頻率）角度看，有下式注：，其中。其matlab近似為X=fft(x,N)/sqrt(N)，在下

2、文中就是指matlab fft函數的計算結果了使用關系可以寫成物理頻率的函數，其中是采樣頻率相關序列可以從功率譜用IDFT變換求得：序列在整個Nyquist間隔上的平均功率可以表示為上式中的以及被定義為平穩隨機信號的power spectral density (PSD)（功率譜密度）一個信號在頻帶上的平均功率可以通過對PSD在頻帶上積分求出從上式中可以看出是一個信號在一個無窮小頻帶上的功率濃度，這也是為什么它叫做功率譜密度。PSD的單位是功率（e.g 瓦特）每單位頻率。在的情況下，這是瓦特/弧度/抽或只是瓦特/弧度。在的情況下單位是瓦特/赫茲。PSD對頻率的積分得到的單位是瓦特，正如平均功率

3、所期望的那樣。對實信號，PSD是關于直流信號對稱的，所以的就足夠完整的描述PSD了。然而要獲得整個Nyquist間隔上的平均功率，有必要引入單邊PSD的概念：信號在頻帶上的平均功率可以用單邊PSD求出頻譜估計方法Matlab 信號處理工具箱提供了三種方法 PSD直接從信號本身估計出來。最簡單的就是periodogram（周期圖法），一種改進的周期圖法是Welch's method。更現代的一種方法是multitaper method（多椎體法）。Parametric methods （參量類方法）這類方法是假設信號是一個由白噪聲驅動的線性系統的輸出。這類方法的例子是Yule-Walke

4、r autoregressive (AR) method和Burg method。這些方法先估計假設的產生信號的線性系統的參數。這些方法想要對可用數據相對較少的情況產生優于傳統非參數方法的結果。Subspace methods （子空間類）又稱為high-resolution methods（高分辨率法）或者super-resolution methods（超分辨率方法）基于對自相關矩陣的特征分析或者特征值分解產生信號的頻率分量。代表方法有multiple signal classification (MUSIC) method或eigenvector (EV) method。這類方法對線譜（

5、正弦信號的譜）最合適，對檢測噪聲下的正弦信號很有效，特別是低信噪比的情況。Nonparametric Methods非參數法下面討論periodogram, modified periodogram, Welch, 和 multitaper法。同時也討論CPSD函數，傳輸函數估計和相關函數。Periodogram周期圖法一個估計功率譜的簡單方法是直接求隨機過程抽樣的DFT，然后取結果的幅度的平方。這樣的方法叫做周期圖法。一個長L的信號的PSD的周期圖估計是注：這里運用的是matlab里面的fft的定義不帶歸一化系數Matlab FFT函數未做歸一化，所以要除以L其中實際對的計算可以只在有限的頻

6、率點上執行并且使用FFT。實踐上大多數周期圖法的應用都計算N點PSD估計，其中選擇N是大于L的下一個2的冪次是明智的，要計算我們直接對補零到長度為N。假如L>N，在計算前，我們必須繞回模N。作為一個例子，考慮下面1001元素信號，它包含了2個正弦信號和噪聲randn('state',0);fs = 1000; % Sampling frequencyt = (0:fs)/fs; % One second worth of samplesA = 1 2; % Sinusoid amplitudes (row vector)f = 150;140; % Sinusoid fr

7、equencies (column vector)xn = A*sin(2*pi*f*t) + 0.1*randn(size(t);注意：最后三行表明了一個方便的表示正弦之和的方法，它等價于：xn = sin(2*pi*150*t) + 2*sin(2*pi*140*t) + 0.1*randn(size(t); 對這個PSD的周期圖估計可以通過產生一個周期圖對象（periodogram object）來計算Hs = spectrum.periodogram('Hamming');估計的圖形可以用psd函數顯示。psd(Hs,xn,'Fs',fs,'NF

8、FT',1024,'SpectrumType','twosided')平均功率通過用下述求和去近似積分 求得Pxx,F = psd(Hs,xn,fs,'twosided');Pow = (fs/length(Pxx) * sum(Pxx)Pow = 2.5059你還可以用單邊PSD去計算平均功率Pxxo,F = psd(Hs,xn,fs,'onesided');Pow = (fs/(2*length(Pxxo) * sum(Pxxo)Pow = 2.5011周期圖性能下面從四個角度討論周期圖法估計的性能：泄漏，分辨率，偏差

9、和方差。頻譜泄漏考慮有限長信號，把它表示成無限長序列乘以一個有限長矩形窗的乘積的形式經常很有用：因為時域的乘積等效于頻域的卷積，所以上式的傅立葉變換是前文中導出的表達式說明卷積對周期圖有影響。正弦數據的卷積影響最容易理解。假設是M個復正弦的和其頻譜是對一個有限長序列，就變成了所以在有限長信號的頻譜中，Dirac函數被替換成了形式為的項，該項對應于矩形窗的中心在的頻率響應。一個矩形窗的頻率響應形狀是一個sinc信號，如下所示該圖顯示了一個主瓣和若干旁瓣，最大旁瓣大約在主瓣下方13.5dB處。這些旁瓣說明了頻譜泄漏效應。無限長信號的功率嚴格的集中在離散頻率點處，而有限長信號在離散頻率點附近有連續的

10、功率。因為矩形窗越短，它的頻率響應對Dirac沖擊的近似性越差，所以數據越短它的頻譜泄漏越明顯。考慮下面的100個采樣的序列randn('state',0)fs = 1000; % Sampling frequencyt = (0:fs/10)/fs; % One-tenth of a second worth of samplesA = 1 2; % Sinusoid amplitudesf = 150;140; % Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 0.1*randn(size(t);Hs = spectrum.perio

11、dogram;psd(Hs,xn,'Fs',fs,'NFFT',1024)注意到頻譜泄露只視數據長度而定。周期圖確實只對有限數據樣本進行計算，但是這和頻譜泄露無關。分辨率分辨率指的是區分頻譜特征的能力，是分析譜估計性能的關鍵概念。要區分兩個在頻率上離得很近的正弦，要求兩個頻率差大于任何一個信號泄漏頻譜的主瓣寬度。主瓣寬度定義為主瓣上峰值功率一半的點間的距離（3dB帶寬）。該寬度近似等于兩個頻率為的正弦信號，可分辨條件是上例中頻率間隔10Hz，數據長度要大于100抽才能使得周期圖中兩個頻率可分辨。下圖是只有67個數據長度的情況randn('state

12、9;,0)fs = 1000; % Sampling frequencyt = (0:fs/15)./fs; % 67 samplesA = 1 2; % Sinusoid amplitudesf = 150;140; % Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 0.1*randn(size(t);Hs=spectrum.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',1024)上述對分辨率的討論都是在高信噪比的情況進行的，因此沒有考慮噪聲。當信噪比低的時候，譜特征的分辨更難，而且周期圖上會

13、出現一些噪聲的偽像，如下所示randn('state',0)fs = 1000; % Sampling frequencyt = (0:fs/10)./fs; % One-tenth of a second worth of samplesA = 1 2; % Sinusoid amplitudesf = 150;140; % Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 2*randn(size(t);Hs=spectrum.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',102

14、4)估計偏差周期圖是對PSD的有偏估計。期望值可以是該式和頻譜泄漏中的式相似，除了這里的表達式用的是平均功率而不是幅度。這暗示了周期圖產生的估計對應于一個有泄漏的PSD而非真正的PSD。注意本質上是一個三角Bartlett窗（事實是兩個矩形脈沖的卷積是三角脈沖。）這導致了最大旁瓣峰值比主瓣峰值低27dB，大致是非平方矩形窗的2倍。周期圖估計是漸進無偏的。這從早期的一個觀察結果可以明顯看出，隨著記錄數據趨于無窮大，矩形窗對頻譜對Dirac函數的近似也就越來越好。然而在某些情況下，周期圖法估計很差勁即使數據夠長，這是因為周期圖法的方差，如下所述。周期圖法的方差L趨于無窮大，方差也不趨于0。用統計學

15、術語講，該估計不是無偏估計。然而周期圖在信噪比大的時候仍然是有用的譜估計器，特別是數據夠長。Modified Periodogram修正周期圖法在fft前先加窗，平滑數據的邊緣。可以降低旁瓣的高度。旁瓣是使用矩形窗產生的陡峭的剪切引入的寄生頻率，對于非矩形窗，結束點衰減的平滑，所以引入較小的寄生頻率。但是，非矩形窗增寬了主瓣，因此降低了頻譜分辨率。函數periodogram允許指定對數據加的窗，例如默認的矩形窗和Hamming窗randn('state',0)fs = 1000; % Sampling frequencyt = (0:fs/10)./fs; % One-tent

16、h of a second worth of samplesA = 1 2; % Sinusoid amplitudesf = 150;140; % Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 0.1*randn(size(t);Hrect = spectrum.periodogram;psd(Hrect,xn,'Fs',fs,'NFFT',1024);Hhamm = spectrum.periodogram('Hamming');psd(Hhamm,xn,'Fs',fs,'NFF

17、T',1024);事實上加Hamming窗后信號的主瓣大約是矩形窗主瓣的2倍。對固定長度信號，Hamming窗能達到的譜估計分辨率大約是矩形窗分辨率的一半。這種沖突可以在某種程度上被變化窗所解決，例如Kaiser窗。非矩形窗會影響信號的功率，因為一些采樣被削弱了。為了解決這個問題函數periodogram將窗歸一化，有平均單位功率。這樣的窗不影響信號的平均功率。修正周期圖法估計的PSD是其中U是窗歸一化常數假如U保證估計是漸進無偏的。Welch法包括：將數據序列劃分為不同的段（可以有重疊），對每段進行改進周期圖法估計，再平均。用spectrum.welch對象，或pwelch函數。默認

18、情況下數據劃分為4段，50%重疊，應用Hamming窗。取平均的目的是減小方差，重疊會引入冗余但是加Hamming窗可以部分消除這些冗余，因為窗給邊緣數據的權重比較小。數據段的縮短和非矩形窗的使用使得頻譜分辨率下降。下面的例子展示Welch法的折衷。首先用周期圖法估計一個小信噪比下信號的PSD：randn('state',1)fs = 1000; % Sampling frequencyt = (0:0.3*fs)./fs; % 301 samplesA = 2 8; % Sinusoid amplitudes (row vector)f = 150;140; % Sinuso

19、id frequencies (column vector)xn = A*sin(2*pi*f*t) + 5*randn(size(t);Hs = spectrum.periodogram('rectangular')psd(Hs,xn,'Fs',fs,'NFFT',1024);可以看出由于噪聲太大，150Hz正弦信號已經無法識別。Hs = spectrum.welch('rectangular',150,50);psd(Hs,xn,'Fs',fs,'NFFT',512)可以看出兩個信號峰，但是如果

20、進一步削減方差，主瓣增寬也使得信號不可識別。Hs = spectrum.welch('rectangular',100,75);psd(Hs,xn,'Fs',fs,'NFFT',512);Welch法的偏差其中是分段數據的長度，是窗歸一化常數。對一定長度的數據，Welch法估計的偏差會大于周期圖法，因為方差比較難以量化，因為它和分段長以及實用的窗都有關系，但是總的說方差反比于使用的段數。Multitaper Method多椎體法周期圖法估計可以用濾波器組來表示。L個帶通濾波器對信號進行濾波，每個濾波器的3dB帶寬是。所有濾波器的幅度響應相似于矩形

21、窗的幅度響應。周期圖估計就是對每個濾波器輸出信號功率的計算，僅僅使用輸出信號的一個采樣點計算輸出信號功率，而且假設的PSD在每個濾波器的頻帶上是常數。信號長度增加，帶通濾波器的帶寬就在減少，近似度就更好。但是有兩個原因對精確度有影響：1矩形窗對應的帶通濾波器性能很差2每個帶通濾波器輸出信號功率的計算僅僅使用一個采樣點，這使得估計很粗糙。Welch法也可以用濾波器組給出相似的解釋。在Welch法中使用了多個點來計算輸出功率，降低了估計的方差。另一方面每個帶通濾波器的帶寬增大了，分辨率下降了。Thompson的多椎體法（MTM）構建在上述結論之上，提供更優的PSD估計。MTM方法沒有使用帶通濾波器

22、（它們本質上是矩形窗，如同周期圖法中一樣），而是使用一組最優濾波器計算估計值。這些最優FIR濾波器是由一組被叫做離散扁平類球體序列（DPSS，也叫做Slepian序列）得到的。除此之外，MTM方法提供了一個時間-帶寬參數，有了它能在估計方差和分辨率之間進行平衡。該參數由時間-帶寬乘積得到，NW，同時它直接與譜估計的多椎體數有關。總有2*NW-1個多椎體被用來形成估計。這就意味著，隨著NW的提高，會有越來越多的功率譜估計值，估計方差會越來越小。然而，每個多椎體的帶寬仍然正比于NW，因而NE提高，每個估計會存在更大的泄露，從而整體估計會更加呈現有偏。對每一組數據，總有一個NW值能在估計偏差和方差見

23、獲得最好的折中。信號處理工具箱中實現MTM方法的函數是pmtm而實現該方法的對象是spectrum.mtm。下面使用spectrum.mtm來計算前一個例子中的PSD：randn('state',0)fs = 1000; % Sampling frequencyt = (0:fs)/fs; % One second worth of samplesA = 1 2; % Sinusoid amplitudesf = 150;140; % Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 0.1*randn(size(t);Hs1 = spect

24、rum.mtm(4,'adapt');psd(Hs1,xn,'Fs',fs,'NFFT',1024)通過降低時間-帶寬積，能夠提高分辨率。Hs2 = spectrum.mtm(3/2,'adapt');psd(Hs2,xn,'Fs',fs,'NFFT',1024)注意到兩個例子中平均功率都被保留：Hs1p = psd(Hs1,xn,'Fs',fs,'NFFT',1024);Pow1 = avgpower(Hs1p)Pow1 = 2.4926Hs2p = psd(Hs2

25、,xn,'Fs',fs,'NFFT',1024);Pow2 = avgpower(Hs2p)Pow2 = 2.4927這中方法相比Weich方法計算復雜度更高，這是計算離散扁平類球體序列的代價。對于長數據序列（10000點以上），通常計算一次DPSS序列并將其存為MAT文件更加實用。Matlab在dpss.mat中提供了dpsssave、dpssload、dpssdir和 dpssclear供使用。互譜密度函數PSD是互譜密度（CPSD）函數的一個特例，CPSD由兩個信號xn、yn如下定義：如同互相關與協方差的例子，工具箱估計PSD和CPSD是因為信號長度有限。

26、為了使用Welch方法估計相隔等長信號x和y的互功率譜密度，cpsd函數通過將x的FFT和y的FFT再共軛之后相乘的方式得到周期圖。與實值PSD不同，CPSD是個復數函數。cpsd如同pwelch函數一樣處理信號的分段和加窗問題：Sxy = cpsd(x, y, nwin, noverlap, nfft, fs)傳輸函數估計Welch方法的一個應用是非參數系統的識別。假設H是一個線性時不變系統，x(n)和y(n)是H的輸入和輸出。則x(n)的功率譜就與x(n)和y(n)的CPSD通過如下方式相關聯：x(n)和y(n)的一個傳輸函數是：該方法同時估計出幅度和相位信息。tfestimate函數使用

27、Welch方法計算CPSD和功率譜，然后得到他們的商作為傳輸函數的估計值。tfestimate函數使用方法和cpsd相同：將信號x(n)通過FIR濾波器，再畫出實際的幅度響應和估計響應如下：h = ones(1,10)/10; % Moving-average filteryn = filter(h,1,xn);HEST,f = tfestimate(xn,yn,256,128,256,fs);H = freqz(h,1,f,fs);subplot(2,1,1); plot(f,abs(H); title('Actual Transfer Function Magnitude'

28、); subplot(2,1,2); plot(f,abs(HEST);title('Transfer Function Magnitude Estimate'); xlabel('Frequency (Hz)');相干函數兩個信號幅度平方相干性如下所示：該商是一個0到1之間的實數，表征了x(n)和y(n)之間的相干性。mscohere函數輸入兩個序列x和y，計算其功率譜和CPSD，返回CPSD幅度平方與兩個功率譜乘積的商。函數的選項和操作與cpsd和tfestimate相類似。x和濾波器輸出y的相干函數如下：mscohere(xn, yn, 256, 128,

29、 256, fs)如果輸入序列長度nfft，窗長度window，一個窗中重疊的數據點為numoverlap，這樣的話mscohere只對一個樣本操作，函數返回全1。這是因為相干函數對線性獨立數據值為1Parametric Methods參數法參數法在信號長度較短時能夠獲得比非參數法更高的分辨率。這類方法使用不同的方式來估計頻譜：不是試圖直接從數據中估計PSD，而是將數據建模成一個由白噪聲驅動的線性系統的輸出，并試圖估計出該系統的參數。最常用的線性系統模型是全極點模型，也就是一個濾波器，它的所有零點都在z平面的原點。這樣一個濾波器輸入白噪聲后的輸出是一個自回歸（AR）過程。正是由于這個原因，這一

30、類方法被稱作AR方法。AR方法便于描述譜呈現尖峰的數據，即PSD在某些頻點特別大。在很多實際應用中（如語音信號）數據都具有帶尖峰的譜，所以AR模型通常會很有用。另外，AR模型具有相對易于求解的系統線性方程。信號處理工具箱提供了下列AR譜估計方法：l Yule-Walkerl AR method (autocorrelation method)Burgl methodCovariancel methodModifiedl covariance method所有的AR方法都會給出如下表示的PSD估計：不同的AR方法估計AR參數ap(k)稍有不同，從而得到不一樣的PSD估計。下表對各種AR方法做了一

31、個總結：Burg協方差修正協方差Yule-Walker特點對數據不加窗；在最小二乘意義上最小化前向后向預測誤差，限定AR系數以滿足L-D遞歸；對數據不加窗；在最小二乘意義上最小化前向后向預測誤差；對數據不加窗；在最小二乘意義上最小化前向后向預測誤差；對數據加窗；在最小二乘意義上最小化前向預測誤差（也叫自相關法）；優點對短數據具有高分辨率；模型總是穩定；對短數據比Y-W有更好分辨率（估計更準確）；能夠從包含p或更多純正弦信號的數據中提取頻率；對短數據具有高分辨率；能夠從包含p或更多純正弦信號的數據中提取頻率；沒有譜線分裂；對大數據性能與其他相當；模型總是穩定；缺點峰值位置高度依賴于初始相位；在正

32、弦信號包含噪聲或階數很高時可能出現譜線分裂；正弦信號估計頻偏；模型可能不穩定；正弦信號估計頻偏；模型可能不穩定；峰值位置高度依賴于初始相位；正弦信號估計較小頻偏；對于短數據性能不高；正弦信號估計頻偏；非奇異條件階數必須不大于輸入幀尺寸一半；階數必須不大于輸入幀尺寸的三分之二；由于估計有偏，自相關矩陣需要確保正定；Yule-Walker 法Yule-Walker AR法通過計算信號自相關函數的有偏估計、求解前向預測誤差的最小二乘最小化來獲得AR參數。這就得出了Yule-Walker等式。Yule-Walker AR法結果與最大熵估計器結果一致。更多信息參考item 2 in the Select

33、ed Bibliography。由于自相關函數的有偏估計的使用，確保了上述自相關矩陣正定。因此，矩陣可逆且方程一定有解。另外，這樣計算的AR參數總會產生一個穩定的全極點模型。Yule-Walker方程通過Levinson算法可以高效的求解。工具箱中的對象spectrum.yulear和函數pyulear實現了Tule-Walker方法。下例比較了一個語音信號通過Welch法和Yule-Walker法的譜：load mtlbHwelch = spectrum.welch('hamming',256,50);psd(Hwelch,mtlb,'Fs',Fs,'

34、NFFT',1024)Hyulear = spectrum.yulear(14);psd(Hyulear,mtlb,'Fs',Fs,'NFFT',1024)Yule-Walker AR法的譜比周期圖法更加平滑。這是因為其內在的簡單全極點模型的緣故。Burg法Burg AR法譜估計是基于最小化前向后向預測誤差的同時滿足Levinson-Durbin遞歸（參考Marple3，Chapter 7，Proakis6，Section ）。對比與其它的AR估計方法，Burg法避免了對自相關函數的計算，改而直接估計反射系數。Burg法最首要的優勢在于解決含有低噪聲的間

35、隔緊密的正弦信號，并且對短數據的估計，在這種情況下AR功率譜密度估計非常逼近與真值。另外，Burg法確保產生一個穩定AR模型，并且能高效計算。Burg法的精度在階數高、數據記錄長、信噪比高（這會導致線分裂、或者在譜估計中產生無關峰）的情況下較低。Burg法計算的譜密度估計也易受噪聲正弦信號初始相位導致的頻率偏移（相對于真實頻率）影響。這一效應在分析短數據序列時會被放大。工具箱中的spectrum.burg對象和pburg函數實現了Burg法。比較下對于語音信號通過Burg法和Yule-Walker法得到的譜，在較長信號數據的情況下它們非常相似。load mtlbHburg = spectrum

36、.burg(14); % 14th order modelpsd(Hburg,mtlb(1:512),'Fs',Fs,'NFFT',1024)Hyulear = spectrum.yulear(14); % 14th order modelpsd(Hyulear,mtlb(1:512),'Fs',Fs,'NFFT',1024)比較受噪聲干擾的信號的譜，分別使用Burg法和Welch法計算：randn('state',0)fs = 1000; % Sampling frequencyt = (0:fs)/fs; %

37、One second worth of samplesA = 1 2; % Sinusoid amplitudesf = 150;140; % Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 0.1*randn(size(t);Hwelch = spectrum.welch('hamming',256,50);psd(Hwelch,xn,'Fs',fs,'NFFT',1024)Hburg = spectrum.burg(14);psd(Hburg,xn,'Fs',fs,'NFFT',1024)需要注意的是，隨著Burg法模型階數的降低，由于正弦信號初始相位帶來的頻率偏移逐漸明顯。Covariance & Modified Covarian