頻譜分析完整版

1、Mat l ab信號處理工具箱幫助文檔譜估計專題頻譜分析Spectral estimation （譜估計）的目標是基于一個有限的數(shù)據(jù)集合描 述一個信號的功率（在頻率上的）分布。功率譜估計在很多場合下都是有 用的，包括對寬帶噪聲湮沒下的信號的檢測。從數(shù)學上看，一個平穩(wěn)隨機過程x的power spectrum （功率譜）和 correlat ion sequenceC相關序列）通過 d i screte-time Four ier transform （離散時間傅立葉變換）構成聯(lián)系。從normal ized frequency （歸一化角 頻率）角度看，有下式注：S 0)= X0)2，其中 X (

2、w ) = limN * NN/2 x ejw -k k 。其nn =- N /2mat l ab 近似為 X=fft (x, N) /sqrt (N)，在下文中 X(f)就是指 mat l ab fftxx函數(shù)的計算結(jié)果了使用關系=2兀f / f可以寫成物理頻率f的函數(shù)，其中匕是采樣頻率相關序列可以從功率譜用IDFT變換求得:序列x在整個Nyquist間隔上的平均功率可以表示為上式中的P 0)=史電以及P (f )=史XX2 冗xxfs被定義為平穩(wěn)隨機信號X的power spectral density (PSD)(功率譜 密度)一個信號在頻帶氣,氣,0巴 也丸上的平均功率可以通過對PSD在

3、 頻帶上積分求出從上式中可以看出Px (s)是一個信號在一個無窮小頻帶上的功率濃度， 這也是為什么它叫做功率譜密度。PSD的單位是功率(e.g瓦特)每單位頻率。在p (s)的情況下，這是 瓦特/弧度/抽或只是瓦特/弧度。在P (f )的情況下單位是瓦特/赫茲。PSD 對頻率的積分得到的單位是瓦特，正如平均功率p 所期望的那樣。g對實信號，PSD是關于直流信號對稱的，所以0s k的P (s)就足夠 完整的描述PSD 了。然而要獲得整個Nyqu i st間隔上的平均功率，有必要 引入單邊PSD的概念：信號在頻帶s ,s ,0 s s N，在計算Xl f 前，我們必須繞回Xl n 模N。作為一個例子

4、，考慮下面1001元素信號x，它包含了 2個正弦信號和randn(state,0);fs = 1000;t = (0:fs)/fs;A = 1 2;f = 150;140;噪聲% Sampling frequency% One second worth of samples% Sinusoid amplitudes (row vector)% Sinusoid frequencies (column vector) xn = A*sin(2*pi*f*t) + 0.1*randn(size(t);注意：最后三行表明了一個方便的表示正弦之和的方法，它等價于：xn = sin(2*pi*150*t

5、) + 2*sin(2*pi*140*t) + 0.1*randn(size(t);對這個PSD的周期圖估計可以通過產(chǎn)生一個周期圖對象(periodogram object)來計算Hs = spectrum.periodogram(Hamming);估計的圖形可以用psd函數(shù)顯示。psd(Hs,xn,Fs,fs,NFFT,1024,SpectrumType,twosided)平均功率通過用下述求和去近似積分求得Pxx,F = psd(Hs,xn,fs,twosided); Pow = (fs/length(Pxx) * sum(Pxx) Pow = 2.5059你還可以用單邊PSD去計算平均功

6、率Pxxo,F = psd(Hs,xn,fs,onesided);Pow = (fs/(2*length(Pxxo) * sum(Pxxo)Pow = 2.5011周期圖性能下面從四個角度討論周期圖法估計的性能：泄漏，分辨率，偏差和方差。頻譜泄漏考慮有限長信號x n,把它表示成無限長序列xn乘以一個有限長矩 形窗 R的乘積的形式經(jīng)常很有用：因為時域的乘積等效于頻域的卷積，所以上式的傅立葉變換是前文中導出的表達式說明卷積對周期圖有影響。正弦數(shù)據(jù)的卷積影響最容易理解。假設xn是M個復正弦的和其頻譜是對一個有限長序列，就變成了所以在有限長信號的頻譜中，Dirac函數(shù)被替換成了形式為W (f - f

7、)Rk的項，該項對應于矩形窗的中心在fk的頻率響應。個矩形窗的頻率響應形狀是一個s i nc信號，如下所示該圖顯示了一個主瓣和若干旁瓣，最大旁瓣大約在主瓣下方13. 5dB處。 這些旁瓣說明了頻譜泄漏效應。無限長信號的功率嚴格的集中在離散頻率 點fk處，而有限長信號在離散頻率點fk附近有連續(xù)的功率。因為矩形窗越短，它的頻率響應對Dirac沖擊的近似性越差，所以數(shù) 據(jù)越短它的頻譜泄漏越明顯?？紤]下面的100個采樣的序列randn(state,0) fs = 1000;t = (0:fs/10)/fs;A = 1 2;f = 150;140;% Sampling frequency% One-te

8、nth of a second worth of samples% Sinusoid amplitudes% Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 0.1*randn(size(t);Hs = spectrum.periodogram;psd(Hs,xn,Fs,fs,NFFT,1024)注重到頻譜泄露只視數(shù)據(jù)長度而定。周期圖確實只對有限數(shù)據(jù)樣本進 行計算，但是這和頻譜泄露無關。分辨率分辨率指的是區(qū)分頻譜特征的水平，是分析譜估計性能的關鍵概念。要區(qū)分兩個在頻率上離得很近的正弦，要求兩個頻率差大于任何一個 信號泄漏頻譜的主瓣寬度。主瓣寬度定義為主瓣上

9、峰值功率一半的點間的 距離（3dB帶寬）。該寬度近似等于f L兩個頻率為f 1 f2的正弦信號，可分辨條件是上例中頻率間隔10Hz，數(shù)據(jù)長度要大于100抽才能使得周期圖中兩個頻率可分辨。下圖是只有67個數(shù)據(jù)長度的情況randn(state,0) fs = 1000;t = (0:fs/15)./fs;A = 1 2;f = 150;140;% Sampling frequency% 67 samples% Sinusoid amplitudes% Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 0.1*randn(size(t);Hs=spectrum.p

10、eriodogram;psd(Hs,xn,Fs,fs,NFFT,1024)上述對分辨率的討論都是在高信噪比的情況進行的，因此沒有考慮噪 聲。當信噪比低的時候，譜特征的分辨更難，而且周期圖上會出現(xiàn)一些噪 聲的偽像，如下所示randn(state,0) fs = 1000;t = (0:fs/10)./fs;A = 1 2;f = 150;140;% Sampling frequency% One-tenth of a second worth of samples% Sinusoid amplitudes% Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 2

11、*randn(size(t);Hs=spectrum.periodogram;psd(Hs,xn,Fs,fs,NFFT,1024)估計偏差周期圖是對PSD的有偏估計。期望值可以是該式和頻譜泄漏中的xl （f ）式相似，除了這里的表達式用的是平均功 率而不是幅度。這暗示了周期圖產(chǎn)生的估計對應于一個有泄漏的PSD而非 真正的PSD。注重WR（f-P）F本質(zhì)上是一個三角Bart lett窗（事實是兩個矩形脈沖 的卷積是三角脈沖。）這導致了最大旁瓣峰值比主瓣峰值低27dB，大致是 非平方矩形窗的2倍。周期圖估計是漸進無偏的。這從早期的一個觀察結(jié)果可以明顯看出, 隨著記錄數(shù)據(jù)趨于無窮大，矩形窗對頻譜對D

12、irac函數(shù)的近似也就越來越 好。然而在某些情況下，周期圖法估計很差勁即使數(shù)據(jù)夠長，這是因為周 期圖法的方差，如下所述。周期圖法的方差L趨于無窮大，方差也不趨于0。用統(tǒng)計學術語講，該估計不是無偏估 計。然而周期圖在信噪比大的時候仍然是有用的譜估計器，特別是數(shù)據(jù)夠 長。Mod ified Per i odogram 修正周期圖法在fft前先加窗，平滑數(shù)據(jù)的邊緣?？梢越档团园甑母叨?。旁瓣是使用矩形窗產(chǎn)生的陡峭的剪切引入的寄生頻率，對于非矩形窗, 結(jié)束點衰減的平滑，所以引入較小的寄生頻率。但是，非矩形窗增寬了主瓣，因此降低了頻譜分辨率。函數(shù)periodogram允許指定對數(shù)據(jù)加的窗，例如默認的矩形窗

13、和randn(state,0)fs = 1000;t = (0:fs/10)./fs;A = 1 2;f = 150;140;Hammi ng 窗% Sampling frequency% One-tenth of a second worth of samples% Sinusoid amplitudes% Sinusoid frequencies xn = A*sin(2*pi*f*t) + 0.1*randn(size(t);Hrect = spectrum.periodogram;psd(Hrect,xn,Fs,fs,NFFT,1024);Hhamm = spectrum.period

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

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

16、itudes (row vector)% Sinusoid 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正弦信號已經(jīng)無法識別。Hs = spectrum.welch(rectangular,150,50);psd(Hs,xn,Fs,fs,NFFT,512)可以看出兩個信號峰，但是如果進一步削減方差，主瓣增寬也使得信 號不可識別。Hs = spec

17、trum.welch(rectangular,100,75);psd(Hs,xn,Fs,fs,NFFT,512);We l ch法的偏差其中、是分段數(shù)據(jù)的長度，U = 111 |w(n)|2是窗歸一化常數(shù)。對一定長度的數(shù)據(jù)，Welch法估計的偏差會大于周期圖法，因為L L方差比較難以量化，因為它和分段長以及實用的窗都有關系，但是總 的說方差反比于使用的段數(shù)。Mu ltitaper Method 多椎體法周期圖法估計可以用濾波器組來表示。L個帶通濾波器對信號x屋進 L行濾波，每個濾波器的3dB帶寬是七/ L。所有濾波器的幅度響應相似于矩 形窗的幅度響應。周期圖估計就是對每個濾波器輸出信號功率的計

18、算，僅 僅使用輸出信號的一個采樣點計算輸出信號功率，而且假設x R的PSD在L每個濾波器的頻帶上是常數(shù)。信號長度增加，帶通濾波器的帶寬就在減少，近似度就更好。但是有 兩個原因?qū)_度有影響：1矩形窗對應的帶通濾波器性能很差2每個帶 通濾波器輸出信號功率的計算僅僅使用一個采樣點，這使得估計很粗糙。Welch法也可以用濾波器組給出相似的解釋。在Welch法中使用了多個 點來計算輸出功率，降低了估計的方差。另一方面每個帶通濾波器的帶寬 增大了，分辨率下降了。Thompson的多椎體法（MTM）構建在上述結(jié)論之上，提供更優(yōu)的PSD 估計。MTM方法沒有使用帶通濾波器（它們本質(zhì)上是矩形窗，如同周期圖 法

19、中一樣），而是使用一組最優(yōu)濾波器計算估計值。這些最優(yōu)FIR濾波器 是由一組被叫做離散扁平類球體序列（DPSS，也叫做Slepian序列）得到 的。除此之外，MTM方法提供了一個時間-帶寬參數(shù)，有了它能在估計方差和分辨率之間進行平衡。該參數(shù)由時間-帶寬乘積得到，NW,同時它直接 與譜估計的多椎體數(shù)有關?？傆?*NW-1個多椎體被用來形成估計。這就 意味著，隨著NW的提升，會有越來越多的功率譜估計值，估計方差會越 來越小。然而，每個多椎體的帶寬仍然正比于NW，因而NE提升，每個估 計會存在更大的泄露，從而整體估計會更加呈現(xiàn)有偏。對每一組數(shù)據(jù)，總 有一個NW值能在估計偏差和方差見獲得最好的折中。ra

20、ndn(state,0) fs = 1000; t = (0:fs)/fs; A = 1 2; f = 150；140;% Sampling frequency信號處理工具箱中實現(xiàn)MTM方法的函數(shù)是pmtm而實現(xiàn)該方法的對象是 spectrum. mtmo下面使用spectrum. mtm來計算前 例子中的PSD：% One second worth of samples% Sinusoid amplitudes% Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 0.1*randn(size(t);Hs1 = spectrum.mtm(4,adapt);p

21、sd(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)Pow12. 4926Hs2p = psd(Hs2,xn,Fs,fs,NFFT,1024);Pow2 = avgpower(Hs2p)Pow2 =2. 4927這中方法相比Weich方法計算復雜度更高，這是計算離散扁平類球體 序列的代價。對于長數(shù)據(jù)序

22、列(10000點以上)，通常計算一次DPSS序列 并將其存為MAT文件更加實用。Matlab在dpss.mat中提供了 dpsssave、 dpssload、dpssd i r 和 dpssc l ear 供使用?；プV密度函數(shù)PSD是互譜密度(CPSD)函數(shù)的一個特例，CPSD由兩個信號xn、yn如 下定義：如同互相關與協(xié)方差的例子，工具箱估計PSD和CPSD是因為信號長度 有限。為了使用Welch方法估計相隔等長信號x和y的互功率譜密度，cpsd 函數(shù)通過將x的FFT和y的FFT再共軛之后相乘的方式得到周期圖。與實 值PSD不同，CPSD是個復數(shù)函數(shù)。cpsd如同pwe l ch函數(shù)一樣處理

23、信號的 分段和加窗問題：Sxy = cpsd(x, y, nwin, noverlap, nfft, fs)傳輸函數(shù)估計Welch方法的一個應用是非參數(shù)系統(tǒng)的識別。假設H是一個線性時不變 系統(tǒng)，x(n )和y(n )是H的輸入和輸出。則x(n)的功率譜就與x(n)和y(n)的CPSD通過如下方式相關聯(lián)：x(n)和y(n)的一個傳輸函數(shù)是：該方法同時估計出幅度和相位信息。tfestimate函數(shù)使用Welch方法 計算CPSD和功率譜，然后得到他們的商作為傳輸函數(shù)的估計值。 tfestimate函數(shù)使用方法和cpsd相同：將信號x(n)通過FIR濾波器，再畫出實際的幅度響應和估計響應如下:h =

24、 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);subplot(2,1,2); plot(f,abs(HEST);title(Transfer Function Magnitude Estimate);xlabel(Frequency (Hz);相干函數(shù)兩個信號幅

25、度平方相干性如下所示：該商是一個0到1之間的實數(shù)，表征了 x(n)和y(n )之間的相干性。mscohere函數(shù)輸入兩個序列x和y，計算其功率譜和CPSD，返回CPSD 幅度平方與兩個功率譜乘積的商。函數(shù)的選項和操作與cpsd和tfestimate 相類似。x和濾波器輸出y的相干函數(shù)如下:mscohere(xn, yn, 256, 128, 256, fs)如果輸入序列長度nfft,窗長度window,個窗中重疊的數(shù)據(jù)點為 numoverlap,這樣的話mscohere 只對一個樣本操作，函數(shù)返回全1。這是 因為相干函數(shù)對線性獨立數(shù)據(jù)值為1Parametr i c Methods 參數(shù)法參數(shù)法

26、在信號長度較短時能夠獲得比非參數(shù)法更高的分辨率。這類方 法使用不同的方式來估計頻譜：不是試圖直接從數(shù)據(jù)中估計PSD，而是將 數(shù)據(jù)建模成一個由白噪聲驅(qū)動的線性系統(tǒng)的輸出，并試圖估計出該系統(tǒng)的 參數(shù)。最常用的線性系統(tǒng)模型是全極點模型，也就是一個濾波器，它的所有 零點都在z平面的原點。這樣一個濾波器輸入白噪聲后的輸出是一個自回 歸(AR)過程。正是由于這個原因，這一類方法被稱作AR方法。AR方法便于描述譜呈現(xiàn)尖峰的數(shù)據(jù)，即PSD在某些頻點特別大。在很 多實際應用中(如語音信號)數(shù)據(jù)都具有帶尖峰的譜，所以AR模型通常 會很有用。另外，AR模型具有相對易于求解的系統(tǒng)線性方程。信號處理工具箱提供了下列AR

27、譜估計方法：Yule-WalkerAR method (autocorrelation method)BurgmethodCovar i ancemethodModifiedcovar i ance method所有的AR方法都會給出如下表示的PSD估計:不同的AR方法估計AR參數(shù)ap（k）稍有不同，從而得到不一樣的PSD估計。下表對各種AR方法做了一個陳述總結(jié):Burg協(xié)方差修正協(xié)方差Yule-Walker特點對數(shù)據(jù)不加窗； 在取小一乘意義 上最小化前向后 向預測誤差，限定 AR系數(shù)以滿足 L-D遞歸；對數(shù)據(jù)不加窗；在最小一乘意義 上最小化前向后 向預測誤差；對數(shù)據(jù)不加窗；在最小二乘意義 上

28、最小化前向后 向預測誤差；對數(shù)據(jù)加窗；在最小二乘意義 上最小化前向預 測誤差（也叫自 相關法）；優(yōu)點對短數(shù)據(jù)具有高分辨率；模型總是穩(wěn)定；對短數(shù)據(jù)比Y-W 有更好分辨率（估 計更準確）；能夠從包含p或 更多純正弦信號 的數(shù)據(jù)中提取頻 率；對短數(shù)據(jù)具有高 分辨率；能夠從包含p或 更多純正弦信號 的數(shù)據(jù)中提取頻 率；沒有譜線分裂；對大數(shù)據(jù)性能與其他相當；模型總是穩(wěn)定；缺點峰值位置高度依 賴于初始相位； 在正弦信號包含 噪聲或階數(shù)很高 時可能出現(xiàn)譜線 分裂；正弦信號估計頻 偏；模型可能不穩(wěn)定； 正弦信號估計頻 偏；模型可能不穩(wěn)定； 峰值位置高度依 賴于初始相位； 正弦信號估計較 小頻偏；對于短數(shù)據(jù)性能

29、 不高；正弦信號估計頻偏；非奇異條件階數(shù)必須不大于 輸入幀尺寸一半；階數(shù)必須不大于 輸入幀尺寸的三 分之二；由于估計有偏， 自相關矩陣需要 確保正定；Yule-Walker 法Yule-Walker AR法通過計算信號自相關函數(shù)的有偏估計、求解前向預 測誤差的最小二乘最小化來獲得AR參數(shù)。這就得出了 Yu l e-Wa l ker等式。Yule-Walker AR法結(jié)果與最大熵估計器結(jié)果一致。更多信息參考item2 i n the Selected Bibli ographyo由于自相關函數(shù)的有偏估計的使用，確保了上述自相關矩陣正定。因 此，矩陣可逆且方程一定有解。另外，這樣計算的AR參數(shù)總會

30、產(chǎn)生一個 穩(wěn)定的全極點模型。Yu l e-Wa l ker方程通過Lev i nson算法可以高效的求解。工具箱中的對象spectrum. yu l ear和函數(shù)pyu l ear實現(xiàn)了 Tu l e-Wa l ker 方法。下例比較了一個語音信號通過We l ch法和Yu l e-Wa l ker法的譜：load mtlbHwelch = spectrum.welch(hamming,256,50);psd(Hwelch,mtlb,Fs,Fs,NFFT,1024)Hyulear = spectrum.yulear(14);psd(Hyulear,mtlb,Fs,Fs,NFFT,1024)Yu

31、le-Walker AR法的譜比周期圖法更加平滑。這是因為其內(nèi)在的簡單 全極點模型的緣故。Burg 法Burg AR法譜估計是基于最小化前向后向預測誤差的同時滿足 Levi nson-Durb i n 遞歸(參考 Marp le3, Chapter 7, Proak is6, Sect ion 12. 3. 3)o對比與其它的AR估計方法，Burg法避免了對自相關函數(shù)的計算, 改而直接估計反射系數(shù)。Burg法最首要的優(yōu)勢在于解決含有低噪聲的間隔緊密的正弦信號，并 且對短數(shù)據(jù)的估計，在這種情況下AR功率譜密度估計非常逼近與真值。 另外，Burg法確保產(chǎn)生一個穩(wěn)定AR模型，并且能高效計算。Burg

32、法的精度在階數(shù)高、數(shù)據(jù)記錄長、信噪比高(這會導致線分裂、 或者在譜估計中產(chǎn)生無關峰)的情況下較低Burg法計算的譜密度估計也 易受噪聲正弦信號初始相位導致的頻率偏移(相對于真實頻率)影響。這 一效應在分析短數(shù)據(jù)序列時會被放大。工具箱中的spectrum. burg對象和pburg函數(shù)實現(xiàn)了 Burg法。比較下 對于語音信號通過Burg法和Yule-Walker法得到的譜，在較長信號數(shù)據(jù) 的情況下它們非常相似。load mtlbHburg = spectrum.burg(14); % 14th order model psd(Hburg,mtlb(1:512),Fs,Fs,NFFT,1024)H

33、yulear = spectrum.yulear(14); % 14th order model psd(Hyulear,mtlb(1:512),Fs,Fs,NFFT,1024)比較受噪聲干擾的信號的譜，分別使用Burg法和Welch法計算:randn(state,0) fs = 1000; t = (0:fs)/fs; A = 1 2; f = 150;140;% Sampling frequency% One second worth of samples% Sinusoid amplitudes% Sinusoid frequenciesxn = A*sin(2*pi*f*t) + 0.

34、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法模型階數(shù)的降低，由于正弦信號初始相位 帶來的頻率偏移逐漸明顯。Covar i ance & Mod ified Covar i ance協(xié)方差和修正協(xié)方差法AR譜估計的協(xié)方差算法基于最小化前向預測誤差而產(chǎn)生。而修正協(xié)方 差算法是基于最小化前向和后向預測誤差而產(chǎn)生。工具箱中的 spectrum. cov對象和pcov函數(shù)，以及spectrum. mcov對象和pmcov函數(shù) 實現(xiàn)了各自