Matlab信號處理工具箱 幫助文檔 譜估計專題翻譯：無名網(wǎng)友 &Lyra頻譜分析Spectralestimation（譜估計）的目標(biāo)是基于一個有限的數(shù)據(jù)集合描述一個信號的功率（在頻率上的）分布。功率譜估計在很多場合下都是有用的，包括對寬帶噪聲湮沒下的信號的檢測。從數(shù)學(xué)上看，一個平穩(wěn)隨機(jī)過程 xn的 powerspectrum（功率譜）和correlationsequence（相關(guān)序列）通過 discrete-timeFouriertransform（離散時間傅立葉變換） 構(gòu)成聯(lián)系。角度看，有下式（離散時間傅立葉變換） 構(gòu)成聯(lián)系。角度看，有下式從normalizedfrequency（歸一化角頻率）Sxx Rxxmejmm注： Sxx 注： Sxx X2，其中matlab近似為X=fft(x,N)/sqrt(N)的計算結(jié)果了1N/2Xlim xnejn 。其NNnN/2n，在下文中 XLf就是指matlabfft函數(shù)使用關(guān)系 2f/fs可以寫成物理頻率 f的函數(shù)，其中 fs是采樣頻率Sxxf Rxxme2jfm/fsm相關(guān)序列可以從功率譜用 IDFT變換求得：RxxmSxx ejmdfsRxxmSxx ejmdfs/2 fs序列 xn在整個 Nyquist間隔上的平均功率可以表示為Rxx0Sxx dfRxx0Sxx dfs/2Sxxf

f

fs/2 sdf上式中的PxxSxx 以及Pxx f Sxx fxxs被定義為平穩(wěn)隨機(jī)信號 xn的powerspectraldensity（PSD）（功率譜密度）一個信號在頻帶 1,2,0 1 2上的平均功率可以通過對 PSD在頻帶上積分求出21P1,2 Pxx d Pxx d12從上式中可以看出 Pxx 是一個信號在一個無窮小頻帶上的功率濃度， 這也是為什么它叫做功率譜密度。PSD的單位是功率（ e.g瓦特）每單位頻率。在 Pxx 的情況下，這是瓦特/弧度/抽或只是瓦特 /弧度。在 Pxxf的情況下單位是瓦特 /赫茲。PSD對頻率的積分得到的單位是瓦特，正如平均功率 P,所期望的那樣。對實(shí)信號， PSD是關(guān)于直流信號對稱的， 所以0的Pxx 就足夠完整的描述 PSD了。然而要獲得整個 Nyquist間隔上的平均功率，有必要引入單邊PSD的概念：00PonesidedTOC\o"1-5"\h\z2Pxx 0信號在頻帶 1,2,0 1 2上的平均功率可以用單邊 PSD求出2P1,2 Ponesided d1頻譜估計方法 cpsdCpsdMatlab信號處理工具箱提供了三種方法Nonparametricmethods（非參量類方法）PSD直接從信號本身估計出來。 最簡單的就是 periodogram（周期圖法） ，一種改進(jìn)的周期圖法是 Welch'smethod。更現(xiàn)代的一種方法是 multitapermethod（多椎體法） 。Parametricmethods（參量類方法）這類方法是假設(shè)信號是一個由白噪聲驅(qū)動的線性系統(tǒng)的輸出。 這類方法的例子是Yule-Walkerautoregressive（AR）method和Burgmethod。這些方法先估計假設(shè)的產(chǎn)生信號的線性系統(tǒng)的參數(shù)。 這些方法想要對可用數(shù)據(jù)相對較少的情況產(chǎn)生優(yōu)于傳統(tǒng)非參數(shù)方法的結(jié)果。

Subspacemethods又稱為high-resolutionmethods（高分辨率法）或者 super-resolutionmethods（超分辨率方法）基于對自相關(guān)矩陣的特征分析或者特征值分解產(chǎn)生信號的頻率分量。代表方法有multiplesignalclassification（MUSIC）method或eigenvector（EV）method。這類方法對線譜（正弦信號的譜）最合適，對檢測噪聲下的正弦信號很有效，特別是低信噪比的情況。方法描述函數(shù)周期圖PSD估計spectrum.periodogram,periodogramWelch重疊，加窗的信號段的平均周期圖spectrum.welch,pwelch,cpsd,tfestimate,mscohere多椎體多個正交窗（稱為錐）的組合做譜估計spectrum.mtm,pmtmYule-WalkerAR時間序列的估計的自相關(guān)函數(shù)計算自回歸（AR）譜估計spectrum.yulear,pyulearBurg通過最小化線性預(yù)測誤差計算自回歸（AR）譜估計spectrum.burg,pburgCovariance（協(xié)方差）通過最小化前向預(yù)測誤差做時間序列的自回歸（ AR）譜估計spectrum.cov,pcov修正協(xié)方差通過最小化前向及后向預(yù)測誤差做時間序列的自回歸（ AR）譜估計spectrum.mcov,pmcovMUSIC多重信號分類spectrum.music,pmusic特征向量法虛譜估計spectrum.eigenvector,peigNonparametricMethods非參數(shù)法下面討論periodogram,modifiedperiodogram,Welch,和multitaper法。同時也討論 CPSD函數(shù)，傳輸函數(shù)估計和相關(guān)函數(shù)。Periodogram周期圖法DFT，DFT，然后取結(jié)果的幅一個長L的信號 x一個長L的信號 xLn的PSD的周期圖估計是P?xx fXLf2

fsL注：這里 XLf運(yùn)用的是 matlab里面的fft的定義不帶歸一化系數(shù)，所以PowerSpectralDensityEstimateviaPeriodogram0PowerSpectralDensityEstimateviaPeriodogram0要除以L其中L1XLfxLne2jfn/fn0實(shí)際對 XLf的計算可以只在有限的頻率點(diǎn)上執(zhí)行并且使用 FFT。實(shí)踐上大多數(shù)周期圖法的應(yīng)用都計算 N點(diǎn)PSD估計2XLfk kfP?xxfk ，fk ,k0,1,,N1fsL N其中L1XLfk xLne2jkn/Nn0選擇N是大于L的下一個 2的冪次是明智的，要計算 XLfk我們直接對xLn補(bǔ)零到長度為 N。假如L>N，在計算 XLfk前，我們必須繞回 xLn模N。randn('state',0);fs=1000;t=(0:fs)/fs;A=[12];frandn('state',0);fs=1000;t=(0:fs)/fs;A=[12];f=[150;140];%Samplingfrequency%Onesecondworthofsamples%Sinusoidamplitudes(rowvector)%Sinusoidfrequencies(columnvector)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的周期圖估計可以通過產(chǎn)生一個周期圖對象( periodogramobject)來計算Hs=spectrum.periodogram('Hamming');估計的圖形可以用 psd函數(shù)顯示。psd(Hs,xn,'Fs',fs,'NFFT',1024,'SpectrumType','twosided')-10-20-30-40-50-60-70-800 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9Frequency(kHz)求得)zH/Bd(ycneuqerf/rewo[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周期圖性能下面從四個角度討論周期圖法估計的性能：泄漏，分辨率，偏差和方差。頻譜泄漏考慮有限長信號 xLn，把它表示成無限長序列 xn乘以一個有限長矩形窗wRn的乘積的形式經(jīng)常很有用：xLnxnwRnfs/21sXLf1X WR fdfsfs/2P?xx f2XLffsL說明卷積對周期圖有影響。正弦數(shù)據(jù)的卷積影響最容易理解。假設(shè) xn是M個復(fù)正弦的和Mxn Akejknk1其頻譜是MXf fs Ak ffkk1對一個有限長序列，就變成了fs/2 M MXLf1 fs Ak fkWRfd AkWRf fkfsfs/2 k1 k1的項，所以在有限長信號的頻譜中， Dirac函數(shù)被替換成了形式為 WRf fk的項，該項對應(yīng)于矩形窗的中心在 fk的頻率響應(yīng)。一個矩形窗的頻率響應(yīng)形狀是一個 sinc信號，如下所示PowerSpectralDensityEstimateviaPeriodogram0PowerSpectralDensityEstimateviaPeriodogram0PowerSpectralDensityEstimateviaPeriodogram0PowerSpectralDensityEstimateviaPeriodogram000zH/tawBdDSP0矩形窗在物理頻率上的功率譜密度-500 -400 -300 -200 -100 0 100 200 300 400 500frequency/Hz00-200zH/tawBdDSP0矩形窗在物理頻率上的功率譜密度-500 -400 -300 -200 -100 0 100 200 300 400 500frequency/Hz00-200-3-405-0-6該圖顯示了一個主瓣和若干旁瓣，最大旁瓣大約在主瓣下方 13.5dB處。這些旁瓣說明了頻譜泄漏效應(yīng)。無限長信號的功率嚴(yán)格的集中在離散頻率點(diǎn) fk處，fk附近有連續(xù)的功率。因為矩形窗越短，它的頻率響應(yīng)對 Dirac沖擊的近似性越差，所以數(shù)據(jù)越短randn('state',0)fs=1000;t=randn('state',0)fs=1000;t=(0:fs/10)/fs;A=[12];f=[150;140];%Samplingfrequency%One-tenthofasecondworthofsamples%Sinusoidamplitudes%Sinusoidfrequenciesxn=A*sin(2*pi*f*t)+0.1*randn(size(t));Hs=spectrum.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',1024)

)zH/Bd(ycneuqerf/rewo-800-10-20-30-40-50-60-700.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5)zH/Bd(ycneuqerf/rewo-800-10-20-30-40-50-60-700.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5Frequency(kHz)注意到頻譜泄露只視數(shù)據(jù)長度而定。周期圖確實(shí)只對有限數(shù)據(jù)樣本進(jìn)行計分辨率指的是區(qū)分頻譜特征的能力，是分析譜估計性能的關(guān)鍵概念。要區(qū)分兩個在頻率上離得很近的正弦， 要求兩個頻率差大于任何一個信號泄3dB。該寬度近似等于 fs/L兩個頻率為 f1f2的正弦信號，可分辨條件是ff1f2 Ls上例中頻率間隔 10Hz，數(shù)據(jù)長度要大于 100抽才能使得周期圖中兩個頻率randn('state',0)fs=1000;trandn('state',0)fs=1000;t=(0:fs/15)./fs;A=[12];f=[150;140];%Samplingfrequency%67samples%Sinusoidamplitudes%Sinusoidfrequenciesxn=A*sin(2*pi*f*t)+0.1*randn(size(t));Hs=spectrum.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',1024))zH/Bd(ycneuqerf/rewoPPowerSpectralDensityEstimateviaPeriodogram0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5Frequency(kHz))zH/Bd(ycneuqerf/rewoPPowerSpectralDensityEstimateviaPeriodogram0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5Frequency(kHz)上述對分辨率的討論都是在高信噪比的情況進(jìn)行的，因此沒有考慮噪聲。當(dāng)信噪比低的時候， 譜特征的分辨更難， 而且周期圖上會出現(xiàn)一些噪聲的偽像， 如下所示randn('state',0)fsrandn('state',0)fs=1000;t=(0:fs/10)./fs;samplesA=[12];f=[150;140];xn=A*sin(2*pi*f*t)+2*randn(size(t));Hs=spectrum.periodogram;%Samplingfrequency%One-tenthofasecondworthof%Sinusoidamplitudes%Sinusoidfrequenciespsd(Hs,xn,'Fs',fs,'NFFT',1024)PowerSpectralDensityEstimateviaPeriodogram-PowerSpectralDensityEstimateviaPeriodogram-#-45-50-55-60-65-70-75-80-85ModifiedCovariancePowerSpectralDensityEstimate0 0.511.5 2Frequency(kHz)2.53 3.5)zH/Bd(ycneuqerf/rewoPSubspaceMethods子空間法MUSIC方法和特征向量分析法spectrum.music對象和pmusic函數(shù)，以及 spectrum.eigenvector對peig函數(shù)提供了兩種相關(guān)的譜分析方法：spectrum.music對象和pmusic函數(shù)提供 Schmidt提出的MUSIC算法。spectrum.eigenvector對象和peig函數(shù)提供 Johnson提出的EV算法。這兩種算法均基于對自相關(guān)矩陣的特征分析，用于對頻率的估計。這種譜分考慮存在于白噪聲中的一些復(fù)正弦信號。 可以把該系統(tǒng)的自相關(guān)矩陣 R寫作S與噪聲自相關(guān)矩陣 W之和：RSW信號自相關(guān)矩陣的特征向量與信號噪聲子空間有緊密聯(lián)系。S的特征向量 v信號自相關(guān)矩陣的特征向量與信號噪聲子空間有緊密聯(lián)系。像信號向量一樣擴(kuò)展成為信號子空間。如果系統(tǒng)包含 M個復(fù)正弦信號，自相關(guān)矩陣階數(shù)為 p，則特征向量 vM1至vp1擴(kuò)展成自相關(guān)矩陣的噪聲子空間。頻率估計函數(shù)特征分析方法通過計算信號和噪聲子空間向量的某些函數(shù)來生成其頻率估計值。MUSIC和EV技術(shù)選擇一個函數(shù)，它在一個輸入正弦信號頻率點(diǎn)上趨于無窮（分母趨于 0）。使用數(shù)字技術(shù)，得到的估計值在感興趣的頻點(diǎn)上具有尖銳的峰值，這就意味著向量中可能沒用無窮大點(diǎn)。MUSIC估計由下面方程所示：Pf 1 1PMUSICf N N 2H2eHf vkvkHef vkefkp1 kp1此處N是特征向量的維數(shù)， ef是復(fù)正弦信號向量：ef1expj2fexpj2f2expj2f4 ??expj2fn1v表示輸入信號則相關(guān)矩陣的特征向量， vk是第k個特征向量； H代表共軛轉(zhuǎn)置。求和中的特征向量對應(yīng)了最小的特征值并張成噪聲空間 （p是信號子空間維度）。表達(dá)式 vkHef等價于一個傅里葉變換（向量 ef由復(fù)指數(shù)組成） 。這一形式對于數(shù)值計算有用，因為 FFT能夠?qū)γ恳粋€ vk進(jìn)行計算，然后幅度平方再被求和。PEVfPEVfN2vkHef/kkp1工具箱中的 pmusic函數(shù)和peig函數(shù)采用 SVD（奇值分解）對信號計算，采用eig函數(shù)來分析自相關(guān)矩陣，并將特征向量歸為信號子空間和噪聲子空間。當(dāng)使用SVD的時候，pmusic和peig并未顯式地計算相關(guān)矩陣， 但是奇異值卻是特征值。先說一下 pwelch函數(shù)吧，相信你也查過了，我就直接把其他人的說法粘貼過來了：[Pxx,f]=pwelch(x,window,noverlap,nfft,fs)Welch方法是一種修正周期圖功率譜密度估計方法， 它通過選取的窗口對數(shù)據(jù)進(jìn)行加窗處理，先分段求功率譜之后再進(jìn)行平均。其中窗口的長度 N表示每次處理的分段數(shù)據(jù)長度， Noverlap是指相鄰兩段數(shù)據(jù)之間的重疊部分長度。 N越大得到的功率譜分辨率越高 (越準(zhǔn)確 )，但方差加大 (及功率譜曲線不太平滑 )；N越小，結(jié)果的方差會變小，但功率譜分辨率較低 (估計結(jié)果不太準(zhǔn)確 )。將信號分為多段，每段之間可以有 overlapping,也可以沒有。每一段加窗每一段做譜分析求平均。關(guān)于信號到底要分多少段？每一段長度多少？我的理解，這個問題的回答，需要考慮頻域分辨率。 當(dāng)然，每一段時間跨度越大，頻域分辨率就越高。如果整個信號不分段， 那么頻域分辨率最高。 當(dāng)然分段也是為了考慮噪聲的影響。 分得段越多，噪聲越小。 有時候，如果一個信號本身不太長， 那么如果有 overlapping,顯然可以比不 overlapping的時候的頻域分辨率高。 總結(jié)一下，每一段的長度，和overlapping的長度，需要平衡噪音，頻域分辨率。這個說法網(wǎng)上很多，直接粘過來了。意思就是 pwelch的輸出的 Pxx是信號的功率譜，而具體怎么個求法呢？下面就用FFT的結(jié)果來求功率譜：fs=10e6;ts=1/fs;N=4000;t=(0:N-1)*ts;f0=1e5;st=cos(2*pi*f0*t)+1;st_fft=fft(st);st_fft1=st_fft(1:end/2+1);psdx=abs(st_fft1).^2/fs/N; %功率譜密度為能量譜密度除時間，模值的平方即為能量譜psdx(2:end-1)=2*psdx(2:end-1); %第一個點(diǎn)是 0頻，最后一個點(diǎn)是第一段多出的點(diǎn)，這兩個點(diǎn)都不是對稱的，因此乘 2的時候把這兩個點(diǎn)去掉figure;plot(psdx);title('fft方法求的功率譜 ');%figure;plot(fftshift(abs(st_fft)));[Pxx,f]=pwelch(st,rectwin(lengt