本文格式為Word版，下載可任意編輯——Matlab頻譜分析程序Matlab信號(hào)處理工具箱譜估計(jì)專題

頻譜分析

Spectralestimation（譜估計(jì)）的目標(biāo)是基于一個(gè)有限的數(shù)據(jù)集合描述一個(gè)信號(hào)的功率（在頻率上的）分布。功率譜估計(jì)在好多場(chǎng)合下都是有用的，包括對(duì)寬帶噪聲湮沒下的信號(hào)的檢測(cè)。

從數(shù)學(xué)上看，一個(gè)平穩(wěn)隨機(jī)過(guò)程xn的powerspectrum（功率譜）和correlationsequence（相關(guān)序列）通過(guò)discrete-timeFouriertransform（離散時(shí)間傅立葉變換）構(gòu)成聯(lián)系。從normalizedfrequency（歸一化角頻率）角度看，有下式

Sxx????m????R?m?exx2??j?m

注：Sxx????X???1，其中X????limN??Nn??N/2?N/2xnej?n??????。其matlab

近似為X=fft(x,N)/sqrt(N)，在下文中XL?f?就是指matlabfft函數(shù)的計(jì)算結(jié)果了

使用關(guān)系??2?f/fs可以寫成物理頻率f的函數(shù)，其中fs是采樣頻率

?Sxx?f??m????R?m?exx?2?jfm/fs

相關(guān)序列可以從功率譜用IDFT變換求得：

Rxx?m??????sSxx???ej?mSxx?f?e2?jfm/fsd???df

2?fs?fs/2f/2序列xn在整個(gè)Nyquist間隔上的平均功率可以表示為

?sSxx???S?f?d???xxdf2?fs?fs/2f/2Rxx?0??上式中的

???Pxx????Sxx???Sxx?f?以及Pxx?f??2?fs被定義為平穩(wěn)隨機(jī)信號(hào)xn的powerspectraldensity(PSD)（功率譜密度）一個(gè)信號(hào)在頻帶??1,?2?,0??1??2??上的平均功率可以通過(guò)對(duì)PSD在頻帶上積分求出

?2P??1,?2??P???d??????xx1??1Pxx???d?

?2從上式中可以看出Pxx???是一個(gè)信號(hào)在一個(gè)無(wú)窮小頻帶上的功率濃度，這也是為什么它叫做功率譜密度。

PSD的單位是功率（e.g瓦特）每單位頻率。在Pxx???的狀況下，這是瓦特/弧度/抽或只是瓦特/弧度。在Pxx?f?的狀況下單位是瓦特/赫茲。PSD對(duì)頻率的積分得到的單位是瓦特，正如平均功率P??,??所期望的那樣。

12對(duì)實(shí)信號(hào)，PSD是關(guān)于直流信號(hào)對(duì)稱的，所以0????的Pxx???就足夠完整的描述PSD了。然而要獲得整個(gè)Nyquist間隔上的平均功率，有必要引入單邊PSD的概念：

?????0?0Ponesided??????2Pxx???0????信號(hào)在頻帶??1,?2?,0??1??2??上的平均功率可以用單邊PSD求出

?2P??1,?2??P??1onesided???d?

頻譜估計(jì)方法

Matlab信號(hào)處理工具箱提供了三種方法

PSD直接從信號(hào)本身估計(jì)出來(lái)。最簡(jiǎn)單的就是periodogram（周期圖法），一種改進(jìn)

的周期圖法是Welch'smethod。更現(xiàn)代的一種方法是multitapermethod（多椎體法）。

Parametricmethods（參量類方法）

這類方法是假設(shè)信號(hào)是一個(gè)由白噪聲驅(qū)動(dòng)的線性系統(tǒng)的輸出。這類方法的例子是Yule-Walkerautoregressive(AR)method和Burgmethod。這些方法先估計(jì)假設(shè)的產(chǎn)生信號(hào)的線性系統(tǒng)的參數(shù)。這些方法想要對(duì)可用數(shù)據(jù)相對(duì)較少的狀況產(chǎn)生優(yōu)于傳統(tǒng)非參數(shù)方法的結(jié)果。

Subspacemethods（子空間類）

又稱為high-resolutionmethods（高分辯率法）或者super-resolutionmethods（超分辯率方法）基于對(duì)自相關(guān)矩陣的特征分析或者特征值分解產(chǎn)生信號(hào)的頻率分量。代表方法有multiplesignalclassification(MUSIC)method或eigenvector(EV)method。這類方法對(duì)線譜（正弦信號(hào)的譜）最適合，對(duì)檢測(cè)噪聲下的正弦信號(hào)很有效，特別是低信噪比的狀況。

NonparametricMethods非參數(shù)法

下面探討periodogram,modifiedperiodogram,Welch,和multitaper法。同時(shí)也探討CPSD函數(shù)，傳輸函數(shù)估計(jì)和相關(guān)函數(shù)。

Periodogram周期圖法

一個(gè)估計(jì)功率譜的簡(jiǎn)單方法是直接求隨機(jī)過(guò)程抽樣的DFT，然后取結(jié)果的幅度的平方。這樣的方法叫做周期圖法。

一個(gè)長(zhǎng)L的信號(hào)xL?n?的PSD的周期圖估計(jì)是

??f??PxxXL?f?fsL2注：這里XL?f?運(yùn)用的是matlab里面的fft的定義不帶歸一化系數(shù)，所以要除以L其中

XL?f???xL?n?e?2?jfn/fs

n?0L?1實(shí)際對(duì)XL?f?的計(jì)算可以只在有限的頻率點(diǎn)上執(zhí)行并且使用FFT。實(shí)踐上大多數(shù)周期圖法的應(yīng)用都計(jì)算N點(diǎn)PSD估計(jì)

??f??PxxkXL?fk?fsL2，fk?kfs,k?0,1,?,N?1N其中

XL?fk???xL?n?e?2?jkn/N

n?0L?1選擇N是大于L的下一個(gè)2的冪次是明智的，要計(jì)算XL?fk?我們直接對(duì)xL?n?補(bǔ)零到長(zhǎng)度為N。假使L>N，在計(jì)算XL?fk?前，我們必需繞回xL?n?模N。

作為一個(gè)例子，考慮下面1001元素信號(hào)xn，它包含了2個(gè)正弦信號(hào)和噪聲r(shí)andn('state',0);fs=1000;%Samplingfrequencyt=(0:fs)/fs;%OnesecondworthofsamplesA=[12];%Sinusoidamplitudes(rowvector)f=[150;140];%Sinusoidfrequencies(columnvector)xn=A*sin(2*pi*f*t)+0.1*randn(size(t));注意：最終三行說(shuō)明白一個(gè)便利的表示正弦之和的方法，它等價(jià)于：xn=sin(2*pi*150*t)+2*sin(2*pi*140*t)+0.1*randn(size(t));對(duì)這個(gè)PSD的周期圖估計(jì)可以通過(guò)產(chǎn)生一個(gè)周期圖對(duì)象（periodogramobject）來(lái)計(jì)算Hs=spectrum.periodogram('Hamming');估計(jì)的圖形可以用psd函數(shù)顯示。psd(Hs,xn,'Fs',fs,'NFFT',1024,'SpectrumType','twosided')PowerSpectralDensityEstimateviaPeriodogram0-10-20Power/frequency(dB/Hz)-30-40-50-60-70-8000.10.20.30.40.50.6Frequency(kHz)0.70.80.9平均功率通過(guò)用下述求和去近似積分求得[Pxx,F]=psd(Hs,xn,fs,'twosided');Pow=(fs/length(Pxx))*sum(Pxx)Pow=2.5059你還可以用單邊PSD去計(jì)算平均功率[Pxxo,F]=psd(Hs,xn,fs,'onesided');Pow=(fs/(2*length(Pxxo)))*sum(Pxxo)Pow=2.5011周期圖性能

下面從四個(gè)角度探討周期圖法估計(jì)的性能：泄漏，分辯率，偏差和方差。頻譜泄漏

考慮有限長(zhǎng)信號(hào)xL?n?，把它表示成無(wú)限長(zhǎng)序列x?n?乘以一個(gè)有限長(zhǎng)矩形窗wR?n?的乘積的形式經(jīng)常很有用：

xL?n??x?n??wR?n?

由于時(shí)域的乘積等效于頻域的卷積，所以上式的傅立葉變換是

1XL?f??fsfs/2?fs/2?X???WR?f???d?

前文中導(dǎo)出的表達(dá)式

??f??PxxXL?f?fsL2

說(shuō)明卷積對(duì)周期圖有影響。

正弦數(shù)據(jù)的卷積影響最簡(jiǎn)單理解。假設(shè)x?n?是M個(gè)復(fù)正弦的和

Mx?n???Akej?kn

k?1其頻譜是

X?f??fs?Ak??f?fk?

k?1M對(duì)一個(gè)有限長(zhǎng)序列，就變成了

1XL?f??fsfs/2?fs/2?fs?Ak????fk?WR?f???d???AkWR?f?fk?

k?1k?1MM所以在有限長(zhǎng)信號(hào)的頻譜中，Dirac函數(shù)被替換成了形式為WR?f?fk?的項(xiàng)，該項(xiàng)對(duì)應(yīng)于矩形窗的中心在fk的頻率響應(yīng)。

一個(gè)矩形窗的頻率響應(yīng)形狀是一個(gè)sinc信號(hào)，如下所示

矩形窗在物理頻率上的功率譜密度0-10-20PSDdBwatt/Hz-30-40-50-60-70-80-500-400-300-200-1000100frequency/Hz202300400500該圖顯示了一個(gè)主瓣和若干旁瓣，最大旁瓣大約在主瓣下方13.5dB處。這些旁瓣說(shuō)明白頻譜泄漏效應(yīng)。無(wú)限長(zhǎng)信號(hào)的功率嚴(yán)格的集中在離散頻率點(diǎn)fk處，而有限長(zhǎng)信號(hào)在離散頻率點(diǎn)fk附近有連續(xù)的功率。

由于矩形窗越短，它的頻率響應(yīng)對(duì)Dirac沖擊的近似性越差，所以數(shù)據(jù)越短它的頻譜泄漏越明顯?？紤]下面的100個(gè)采樣的序列

randn('state',0)fs=1000;%Samplingfrequencyt=(0:fs/10)/fs;%One-tenthofasecondworthofsamplesA=[12];%Sinusoidamplitudesf=[150;140];%Sinusoidfrequenciesxn=A*sin(2*pi*f*t)+0.1*randn(size(t));Hs=spectrum.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',1024)PowerSpectralDensityEstimateviaPeriodogram0-10-20Power/frequency(dB/Hz)-30-40-50-60-70-8000.050.10.150.20.250.3Frequency(kHz)0.350.40.450.5注意到頻譜泄露只視數(shù)據(jù)長(zhǎng)度而定。周期圖確實(shí)只對(duì)有限數(shù)據(jù)樣本進(jìn)行計(jì)算，但是這和頻譜泄露無(wú)關(guān)。分辯率

分辯率指的是區(qū)分頻譜特征的能力，是分析譜估計(jì)性能的關(guān)鍵概念。

要區(qū)分兩個(gè)在頻率上離得很近的正弦，要求兩個(gè)頻率差大于任何一個(gè)信號(hào)泄漏頻譜的主瓣寬度。主瓣寬度定義為主瓣上峰值功率一半的點(diǎn)間的距離（3dB帶寬）。該寬度近似等于

fs/L

兩個(gè)頻率為f1f2的正弦信號(hào)，可分辯條件是

?f??f1?f2??fsL上例中頻率間隔10Hz，數(shù)據(jù)長(zhǎng)度要大于100抽才能使得周期圖中兩個(gè)頻率可分辯。下圖是只有67個(gè)數(shù)據(jù)長(zhǎng)度的狀況

randn('state',0)fs=1000;%Samplingfrequencyt=(0:fs/15)./fs;%67samplesA=[12];%Sinusoidamplitudesf=[150;140];%Sinusoidfrequenciesxn=A*sin(2*pi*f*t)+0.1*randn(size(t));Hs=spectrum.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',1024)PowerSpectralDensityEstimateviaPeriodogram0-10Power/frequency(dB/Hz)-20-30-40-50-60-7000.050.10.150.20.250.3Frequency(kHz)0.350.40.450.5上述對(duì)分辯率的探討都是在高信噪比的狀況進(jìn)行的，因此沒有考慮噪聲。當(dāng)信噪比低的時(shí)候，譜特征的分辯更難，而且周期圖上會(huì)出現(xiàn)一些噪聲的偽像，如下所示

randn('state',0)fs=1000;%Samplingfrequencyt=(0:fs/10)./fs;%One-tenthofasecondworthofsamplesA=[12];%Sinusoidamplitudesf=[150;140];%Sinusoidfrequenciesxn=A*sin(2*pi*f*t)+2*randn(size(t));Hs=spectrum.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',1024)PowerSpectralDensityEstimateviaPeriodogram-5-10-15Power/frequency(dB/Hz)-20-25-30-35-40-45-50-5500.050.10.150.20.250.3Frequency(kHz)0.350.40.450.5估計(jì)偏差周期圖是對(duì)PSD的有偏估計(jì)。期望值可以是

2f/2?Xf2?L????1sE??P?Wf??d??????xxR?fsL?fsL?fs/2???該式和頻譜泄漏中的XL?f?式相像，除了這里的表達(dá)式用的是平均功率而不是幅度。這示意了周期圖產(chǎn)生的估計(jì)對(duì)應(yīng)于一個(gè)有泄漏的PSD而非真正的PSD。

注意WR?f???本質(zhì)上是一個(gè)三角Bartlett窗（事實(shí)是兩個(gè)矩形脈沖的卷積是三角脈沖。）這導(dǎo)致了最大旁瓣峰值比主瓣峰值低27dB，大致是非平方矩形窗的2倍。周期圖估計(jì)是漸進(jìn)無(wú)偏的。這從早期的一個(gè)觀測(cè)結(jié)果可以明顯看出，隨著記錄數(shù)據(jù)趨于

無(wú)窮大，矩形窗對(duì)頻譜對(duì)Dirac函數(shù)的近似也就越來(lái)越好。然而在某些狀況下，周期圖法估計(jì)很差勁即使數(shù)據(jù)夠長(zhǎng)，這是由于周期圖法的方差，如下所述。周期圖法的方差

22???Xf??sin?2?Lf/fs???L???2?var????Pxx?f??1????fsL???Lsin?2?f/fs????????2L趨于無(wú)窮大，方差也不趨于0。用統(tǒng)計(jì)學(xué)術(shù)語(yǔ)講，該估計(jì)不是無(wú)偏估計(jì)。然而周期圖

在信噪比大的時(shí)候依舊是有用的譜估計(jì)器，特別是數(shù)據(jù)夠長(zhǎng)。

ModifiedPeriodogram修正周期圖法

在fft前先加窗，平滑數(shù)據(jù)的邊緣。可以降低旁瓣的高度。

旁瓣是使用矩形窗產(chǎn)生的陡峭的剪切引入的寄生頻率，對(duì)于非矩形窗，終止點(diǎn)衰減的平滑，所以引入較小的寄生頻率。

但是，非矩形窗增寬了主瓣，因此降低了頻譜分辯率。

函數(shù)periodogram允許指定對(duì)數(shù)據(jù)加的窗，例如默認(rèn)的矩形窗和Hamming窗randn('state',0)fs=1000;%Samplingfrequencyt=(0:fs/10)./fs;%One-tenthofasecondworthofsamplesA=[12];%Sinusoidamplitudesf=[150;140];%Sinusoidfrequenciesxn=A*sin(2*pi*f*t)+0.1*randn(size(t));Hrect=spectrum.periodogram;psd(Hrect,xn,'Fs',fs,'NFFT',1024);PowerSpectralDensityEstimateviaPeriodogram0-10-20Power/frequency(dB/Hz)-30-40-50-60-70-8000.050.10.150.20.250.3Frequency(kHz)0.350.40.450.5Hhamm=spectrum.periodogram('Hamming');psd(Hhamm,xn,'Fs',fs,'NFFT',1024);PowerSpectralDensityEstimateviaPeriodogram0-10-20Power/frequency(dB/Hz)-30-40-50-60-70-8000.050.10.150.20.250.3Frequency(kHz)0.350.40.450.5事實(shí)上加Hamming窗后信號(hào)的主瓣大約是矩形窗主瓣的2倍。對(duì)固定長(zhǎng)度信號(hào)，Hamming窗能達(dá)到的譜估計(jì)分辯率大約是矩形窗分辯率的一半。這種沖突可以在某種程度上被變化窗所解決，例如Kaiser窗。

非矩形窗會(huì)影響信號(hào)的功率，由于一些采樣被減弱了。為了解決這個(gè)問(wèn)題

函數(shù)periodogram將窗歸一化，有平均單位功率。這樣的窗不影響信號(hào)的平均功率。修正周期圖法估計(jì)的PSD是

??f??PxxXL?f?fsLU2其中U是窗歸一化常數(shù)

21L?1U??w?n?Ln?0假使U保證估計(jì)是漸進(jìn)無(wú)偏的。

Welch法

包括：將數(shù)據(jù)序列劃分為不同的段（可以有重疊），對(duì)每段進(jìn)行改進(jìn)周期圖法估計(jì)，再平均。

用spectrum.welch對(duì)象，或pwelch函數(shù)。默認(rèn)狀況下數(shù)據(jù)劃分為4段，50%重疊，應(yīng)用Hamming窗。

取平均的目的是減小方差，重疊會(huì)引入冗余但是加Hamming窗可以部分消除這些冗余，由于窗給邊緣數(shù)據(jù)的權(quán)重比較小。

數(shù)據(jù)段的縮短和非矩形窗的使用使得頻譜分辯率下降。下面的例子展示W(wǎng)elch法的折衷。

首先用周期圖法估計(jì)一個(gè)小信噪比下信號(hào)的PSD：

randn('state',1)fs=1000;%Samplingfrequencyt=(0:0.3*fs)./fs;%301samplesA=[28];%Sinusoidamplitudes(rowvector)f=[150;140];%Sinusoidfrequencies(columnvector)xn=A*sin(2*pi*f*t)+5*randn(size(t));Hs=spectrum.periodogram('rectangular')psd(Hs,xn,'Fs',fs,'NFFT',1024);可以看出由于噪聲太大，150Hz正弦信號(hào)已經(jīng)無(wú)法識(shí)別。

PowerSpectralDensityEstimateviaPeriodogram100Power/frequency(dB/Hz)-10-20-30-40-50-6000.050.10.150.20.250.3Frequency(kHz)0.350.40.450.5Hs=spectrum.welch('rectangular',150,50);psd(Hs,xn,'Fs',fs,'NFFT',512)PowerSpectralDensityEstimateviaWelch105Power/frequency(dB/Hz)0-5-10-15-20-2500.050.10.150.20.250.3Frequency(kHz)0.350.40.450.5可以看出兩個(gè)信號(hào)峰，但是假使進(jìn)一步削減方差，主瓣增寬也使得信號(hào)不可識(shí)別。Hs=spectrum.welch('rectangular',100,75);psd(Hs,xn,'Fs',fs,'NFFT',512);PowerSpectralDensityEstimateviaWelch50Power/frequency(dB/Hz)-5-10-15-2000.050.10.150.20.250.3Frequency(kHz)0.350.40.450.5Welch法的偏差

?EPwelch??1?fsLsUfs/2?fs/2?Pxx???WR?f???d?

221L?1其中Ls是分段數(shù)據(jù)的長(zhǎng)度，U??w?n?是窗歸一化常數(shù)。

Ln?0對(duì)一定長(zhǎng)度的數(shù)據(jù)，Welch法估計(jì)的偏差會(huì)大于周期圖法，由于L?Ls

方差比較難以量化，由于它和分段長(zhǎng)以及實(shí)用的窗都有關(guān)系，但是總的說(shuō)方差反比于使用的段數(shù)。

MultitaperMethod多椎體法

周期圖法估計(jì)可以用濾波器組來(lái)表示。L個(gè)帶通濾波器對(duì)信號(hào)xL?n?進(jìn)行濾波，每個(gè)濾波器的3dB帶寬是fs/L。所有濾波器的幅度響應(yīng)相像于矩形窗的幅度響應(yīng)。周期圖估計(jì)就是對(duì)每個(gè)濾波器輸出信號(hào)功率的計(jì)算，僅僅使用輸出信號(hào)的一個(gè)采樣點(diǎn)計(jì)算輸出信號(hào)功率，而且假設(shè)xL?n?的PSD在每個(gè)濾波器的頻帶上是常數(shù)。

信號(hào)長(zhǎng)度增加，帶通濾波器的帶寬就在減少，近似度就更好。但是有兩個(gè)原因?qū)Υ_切度有影響：1矩形窗對(duì)應(yīng)的帶通濾波器性能很差2每個(gè)帶通濾波器輸出信號(hào)功率的計(jì)算僅僅使用一個(gè)采樣點(diǎn)，這使得估計(jì)很粗糙。

Welch法也可以用濾波器組給出相像的解釋。在Welch法中使用了多個(gè)點(diǎn)來(lái)計(jì)算輸出功率，降低了估計(jì)的方差。另一方面每個(gè)帶通濾波器的帶寬增大了，分辯率下降了。Thompson的多椎體法（MTM）構(gòu)建在上述結(jié)論之上，提供更優(yōu)的PSD估計(jì)。MTM方法沒有使用帶通濾波器（它們本質(zhì)上是矩形窗，宛如周期圖法中一樣），而是使用一組最優(yōu)濾波器計(jì)算估計(jì)值。這些最優(yōu)FIR濾波器是由一組被叫做離散扁平類球體序列（DPSS，也叫做Slepian序列）得到的。

除此之外，MTM方法提供了一個(gè)時(shí)間-帶寬參數(shù)，有了它能在估計(jì)方差和分辯率之間進(jìn)行平衡。該參數(shù)由時(shí)間-帶寬乘積得到，NW，同時(shí)它直接與譜估計(jì)的多椎體數(shù)有關(guān)。總有2*NW-1個(gè)多椎體被用來(lái)形成估計(jì)。這就意味著，隨著NW的提高，會(huì)有越來(lái)越多的功率譜估計(jì)值，估計(jì)方差會(huì)越來(lái)越小。然而，每個(gè)多椎體的帶寬依舊正比于NW，因而NE提高，每個(gè)估計(jì)會(huì)存在更大的泄露，從而整體估計(jì)會(huì)更加浮現(xiàn)有偏。對(duì)每一組數(shù)據(jù)，總有一個(gè)NW值能在估計(jì)偏差和方差見獲得最好的折中。

信號(hào)處理工具箱中實(shí)現(xiàn)MTM方法的函數(shù)是pmtm而實(shí)現(xiàn)該方法的對(duì)象是spectrum.mtm。下面使用spectrum.mtm來(lái)計(jì)算前一個(gè)例子中的PSD：randn('state',0)

fs=1000;%Samplingfrequencyt=(0:fs)/fs;%OnesecondworthofsamplesA=[12];%Sinusoidamplitudesf=[150;140];%Sinusoidfrequenciesxn=A*sin(2*pi*f*t)+0.1*randn(size(t));Hs1=spectrum.mtm(4,'adapt');psd(Hs1,xn,'Fs',fs,'NFFT',1024)ThompsonMultitaperPowerSpectralDensityEstimate-5-10-15Power/frequency(dB/Hz)-20-25-30-35-40-45-50-55050100150202350300Frequency(Hz)350400450500

通過(guò)降低時(shí)間-帶寬積，能夠提高分辯率。Hs2=spectrum.mtm(3/2,'adapt');psd(Hs2,xn,'Fs',fs,'NFFT',1024)ThompsonMultitaperPowerSpectralDensityEstimate0-10-20Power/frequency(dB/Hz)-30-40-50-60-70-80050100150202350300Frequency(Hz)350400450500

注意到兩個(gè)例子中平均功率都被保存：Hs1p=psd(Hs1,xn,'Fs',fs,'NFFT',1024);Pow1=avgpower(Hs1p)Pow1=2.4926Hs2p=psd(Hs2,xn,'Fs',fs,'NFFT',1024);Pow2=avgpower(Hs2p)Pow2=2.4927

這中方法相比Weich方法計(jì)算繁雜度更高，這是計(jì)算離散扁平類球體序列的代價(jià)。對(duì)于長(zhǎng)數(shù)據(jù)序列（10000點(diǎn)以上），尋常計(jì)算一次DPSS序列并將其存為MAT文件更加實(shí)用。Matlab在dpss.mat中提供了dpsssave、dpssload、dpssdir和dpssclear供使用?；プV密度函數(shù)

PSD是互譜密度（CPSD）函數(shù)的一個(gè)特例，CPSD由兩個(gè)信號(hào)xn、yn如下定義：

1Pxy????2?m????R???exy??j?m

宛如相互關(guān)與協(xié)方差的例子，工具箱估計(jì)PSD和CPSD是由于信號(hào)長(zhǎng)度有限。為了使用Welch方法估計(jì)相隔等長(zhǎng)信號(hào)x和y的互功率譜密度，cpsd函數(shù)通過(guò)將x的FFT和y

的FFT再共軛之后相乘的方式得到周期圖。與實(shí)值PSD不同，CPSD是個(gè)復(fù)數(shù)函數(shù)。cpsd宛如pwelch函數(shù)一樣處理信號(hào)的分段和加窗問(wèn)題：Sxy=cpsd(x,y,nwin,noverlap,nfft,fs)傳輸函數(shù)估計(jì)Welch方法的一個(gè)應(yīng)用是非參數(shù)系統(tǒng)的識(shí)別。假設(shè)H是一個(gè)線性時(shí)不變系統(tǒng)，x(n)和y(n)是H的輸入和輸出。則x(n)的功率譜就與x(n)和y(n)的CPSD通過(guò)如下方式相關(guān)聯(lián)：

Pyx????H???Pxx???

x(n)和y(n)的一個(gè)傳輸函數(shù)是：

?????Pyx???HPxx???該方法同時(shí)估計(jì)出幅度和相位信息。tfestimate函數(shù)使用Welch方法計(jì)算CPSD和功率譜，然后得到他們的商作為傳輸函數(shù)的估計(jì)值。tfestimate函數(shù)使用方法和cpsd一致：

將信號(hào)x(n)通過(guò)FIR濾波器，再畫出實(shí)際的幅度響應(yīng)和估計(jì)響應(yīng)如下：h=ones(1,10)/10;%Moving-averagefilteryn=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('ActualTransferFunctionMagnitude');subplot(2,1,2);plot(f,abs(HEST));title('TransferFunctionMagnitudeEstimate');xlabel('Frequency(Hz)');ActualTransferFunctionMagnitude10.50050100150202350300350400450500TransferFunctionMagnitudeEstimate1.510.50050100150202350300Frequency(Hz)350400450500

相干函數(shù)

兩個(gè)信號(hào)幅度平方相干性如下所示：

Cxy????Pxy???2Pxx???Pyy???

該商是一個(gè)0到1之間的實(shí)數(shù)，表征了x(n)和y(n)之間的相干性。

mscohere函數(shù)輸入兩個(gè)序列x和y，計(jì)算其功率譜和CPSD，返回CPSD幅度平方與兩個(gè)功率譜乘積的商。函數(shù)的選項(xiàng)和操作與cpsd和tfestimate相類似。

x和濾波器輸出y的相干函數(shù)如下：mscohere(xn,yn,256,128,256,fs)CoherenceEstimateviaWelch10.90.80.7Magnitude(dB)0.60.50.40.30.20.10050100150202350300Frequency(Hz)350400450500

假使輸入序列長(zhǎng)度nfft，窗長(zhǎng)度window，一個(gè)窗中重疊的數(shù)據(jù)點(diǎn)為numoverlap，這樣的話mscohere只對(duì)一個(gè)樣本操作，函數(shù)返回全1。這是由于相干函數(shù)對(duì)線性獨(dú)立數(shù)據(jù)值為1

ParametricMethods參數(shù)法

參數(shù)法在信號(hào)長(zhǎng)度較短時(shí)能夠獲得比非參數(shù)法更高的分辯率。這類方法使用不同的方式來(lái)估計(jì)頻譜：不是試圖直接從數(shù)據(jù)中估計(jì)PSD，而是將數(shù)據(jù)建模成一個(gè)由白噪聲驅(qū)動(dòng)的線性系統(tǒng)的輸出，并試圖估計(jì)出該系統(tǒng)的參數(shù)。

最常用的線性系統(tǒng)模型是全極點(diǎn)模型，也就是一個(gè)濾波器，它的所有零點(diǎn)都在z平面的原點(diǎn)。這樣一個(gè)濾波器輸入白噪聲后的輸出是一個(gè)自回歸（AR）過(guò)程。正是由于這個(gè)原因，這一類方法被稱作AR方法。

AR方法便于描述譜浮現(xiàn)尖峰的數(shù)據(jù)，即PSD在某些頻點(diǎn)特別大。在好多實(shí)際應(yīng)用中（如語(yǔ)音信號(hào)）數(shù)據(jù)都具有帶尖峰的譜，所以AR模型尋常會(huì)很有用。另外，AR模型具有相對(duì)易于求解的系統(tǒng)線性方程。

信號(hào)處理工具箱提供了以下AR譜估計(jì)方法：?????

Yule-Walker

ARmethod(autocorrelationmethod)BurgmethodCovariancemethodModifiedcovariancemethod

BurgPowerSpectralDensityEstimate0-10Power/frequency(dB/Hz)-20-30-40-50-60050100150202350300Frequency(Hz)350400450500

需要注意的是，隨著Burg法模型階數(shù)的降低，由于正弦信號(hào)初始相位帶來(lái)的頻率偏移逐漸明顯。

Covariance%14thordermodelpsd(Hcov,mtlb(1:64),'Fs',Fs,'NFFT',1024)CovariancePowerSpectralDensityEstimate-30-40Power/frequency(dB/Hz)