Matlab頻譜分析程序文件

可編輯版/Matlab信號處理工具箱譜估計(jì)專題頻譜分析Spectralestimation〔譜估計(jì)的目標(biāo)是基于一個(gè)有限的數(shù)據(jù)集合描述一個(gè)信號的功率〔在頻率上的分布。功率譜估計(jì)在很多場合下都是有用的,包括對寬帶噪聲湮沒下的信號的檢測。從數(shù)學(xué)上看,一個(gè)平穩(wěn)隨機(jī)過程的powerspectrum〔功率譜和correlationsequence〔相關(guān)序列通過discrete-timeFouriertransform〔離散時(shí)間傅立葉變換構(gòu)成聯(lián)系。從normalizedfrequency〔歸一化角頻率角度看,有下式注：,其中。其matlab近似為X=fft<x,N>/sqrt<N>,在下文中就是指matlabfft函數(shù)的計(jì)算結(jié)果了使用關(guān)系可以寫成物理頻率的函數(shù),其中是采樣頻率相關(guān)序列可以從功率譜用IDFT變換求得：序列在整個(gè)Nyquist間隔上的平均功率可以表示為上式中的以及被定義為平穩(wěn)隨機(jī)信號的powerspectraldensity<PSD>〔功率譜密度一個(gè)信號在頻帶上的平均功率可以通過對PSD在頻帶上積分求出從上式中可以看出是一個(gè)信號在一個(gè)無窮小頻帶上的功率濃度,這也是為什么它叫做功率譜密度。PSD的單位是功率〔e.g瓦特每單位頻率。在的情況下,這是瓦特/弧度/抽或只是瓦特/弧度。在的情況下單位是瓦特/赫茲。PSD對頻率的積分得到的單位是瓦特,正如平均功率所期望的那樣。對實(shí)信號,PSD是關(guān)于直流信號對稱的,所以的就足夠完整的描述PSD了。然而要獲得整個(gè)Nyquist間隔上的平均功率,有必要引入單邊PSD的概念：信號在頻帶上的平均功率可以用單邊PSD求出頻譜估計(jì)方法Matlab信號處理工具箱提供了三種方法PSD直接從信號本身估計(jì)出來。最簡單的就是periodogram〔周期圖法,一種改進(jìn)的周期圖法是Welch'smethod。更現(xiàn)代的一種方法是multitapermethod〔多椎體法。Parametricmethods〔參量類方法這類方法是假設(shè)信號是一個(gè)由白噪聲驅(qū)動(dòng)的線性系統(tǒng)的輸出。這類方法的例子是Yule-Walkerautoregressive<AR>method和Burgmethod。這些方法先估計(jì)假設(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。這類方法對線譜〔正弦信號的譜最合適,對檢測噪聲下的正弦信號很有效,特別是低信噪比的情況。NonparametricMethods非參數(shù)法下面討論periodogram,modifiedperiodogram,Welch,和multitaper法。同時(shí)也討論CPSD函數(shù),傳輸函數(shù)估計(jì)和相關(guān)函數(shù)。Periodogram周期圖法一個(gè)估計(jì)功率譜的簡單方法是直接求隨機(jī)過程抽樣的DFT,然后取結(jié)果的幅度的平方。這樣的方法叫做周期圖法。一個(gè)長L的信號的PSD的周期圖估計(jì)是注：這里運(yùn)用的是matlab里面的fft的定義不帶歸一化系數(shù)MatlabFFT函數(shù)未做歸一化,所以要除以LMatlabFFT函數(shù)未做歸一化其中實(shí)際對的計(jì)算可以只在有限的頻率點(diǎn)上執(zhí)行并且使用FFT。實(shí)踐上大多數(shù)周期圖法的應(yīng)用都計(jì)算N點(diǎn)PSD估計(jì),其中選擇N是大于L的下一個(gè)2的冪次是明智的,要計(jì)算我們直接對補(bǔ)零到長度為N。假如L>N,在計(jì)算前,我們必須繞回模N。作為一個(gè)例子,考慮下面1001元素信號,它包含了2個(gè)正弦信號和噪聲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>>;注意：最后三行表明了一個(gè)方便的表示正弦之和的方法,它等價(jià)于：xn=sin<2*pi*150*t>+2*sin<2*pi*140*t>+0.1*randn<size<t>>;對這個(gè)PSD的周期圖估計(jì)可以通過產(chǎn)生一個(gè)周期圖對象〔periodogramobject來計(jì)算Hs=spectrum.periodogram<'Hamming'>;估計(jì)的圖形可以用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去計(jì)算平均功率[Pxxo,F]=psd<Hs,xn,fs,'onesided'>;Pow=<fs/<2*length<Pxxo>>>*sum<Pxxo>Pow=2.5011周期圖性能下面從四個(gè)角度討論周期圖法估計(jì)的性能：泄漏,分辨率,偏差和方差。頻譜泄漏考慮有限長信號,把它表示成無限長序列乘以一個(gè)有限長矩形窗的乘積的形式經(jīng)常很有用：因?yàn)闀r(shí)域的乘積等效于頻域的卷積,所以上式的傅立葉變換是前文中導(dǎo)出的表達(dá)式說明卷積對周期圖有影響。正弦數(shù)據(jù)的卷積影響最容易理解。假設(shè)是M個(gè)復(fù)正弦的和其頻譜是對一個(gè)有限長序列,就變成了所以在有限長信號的頻譜中,Dirac函數(shù)被替換成了形式為的項(xiàng),該項(xiàng)對應(yīng)于矩形窗的中心在的頻率響應(yīng)。一個(gè)矩形窗的頻率響應(yīng)形狀是一個(gè)sinc信號,如下所示該圖顯示了一個(gè)主瓣和若干旁瓣,最大旁瓣大約在主瓣下方13.5dB處。這些旁瓣說明了頻譜泄漏效應(yīng)。無限長信號的功率嚴(yán)格的集中在離散頻率點(diǎn)處,而有限長信號在離散頻率點(diǎn)附近有連續(xù)的功率。因?yàn)榫匦未霸蕉?它的頻率響應(yīng)對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>注意到頻譜泄露只視數(shù)據(jù)長度而定。周期圖確實(shí)只對有限數(shù)據(jù)樣本進(jìn)行計(jì)算,但是這和頻譜泄露無關(guān)。分辨率分辨率指的是區(qū)分頻譜特征的能力,是分析譜估計(jì)性能的關(guān)鍵概念。要區(qū)分兩個(gè)在頻率上離得很近的正弦,要求兩個(gè)頻率差大于任何一個(gè)信號泄漏頻譜的主瓣寬度。主瓣寬度定義為主瓣上峰值功率一半的點(diǎn)間的距離〔3dB帶寬。該寬度近似等于兩個(gè)頻率為的正弦信號,可分辨條件是上例中頻率間隔10Hz,數(shù)據(jù)長度要大于100抽才能使得周期圖中兩個(gè)頻率可分辨。下圖是只有67個(gè)數(shù)據(jù)長度的情況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>上述對分辨率的討論都是在高信噪比的情況進(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>估計(jì)偏差周期圖是對PSD的有偏估計(jì)。期望值可以是該式和頻譜泄漏中的式相似,除了這里的表達(dá)式用的是平均功率而不是幅度。這暗示了周期圖產(chǎn)生的估計(jì)對應(yīng)于一個(gè)有泄漏的PSD而非真正的PSD。注意本質(zhì)上是一個(gè)三角Bartlett窗〔事實(shí)是兩個(gè)矩形脈沖的卷積是三角脈沖。這導(dǎo)致了最大旁瓣峰值比主瓣峰值低27dB,大致是非平方矩形窗的2倍。周期圖估計(jì)是漸進(jìn)無偏的。這從早期的一個(gè)觀察結(jié)果可以明顯看出,隨著記錄數(shù)據(jù)趨于無窮大,矩形窗對頻譜對Dirac函數(shù)的近似也就越來越好。然而在某些情況下,周期圖法估計(jì)很差勁即使數(shù)據(jù)夠長,這是因?yàn)橹芷趫D法的方差,如下所述。周期圖法的方差L趨于無窮大,方差也不趨于0。用統(tǒng)計(jì)學(xué)術(shù)語講,該估計(jì)不是無偏估計(jì)。然而周期圖在信噪比大的時(shí)候仍然是有用的譜估計(jì)器,特別是數(shù)據(jù)夠長。ModifiedPeriodogram修正周期圖法在fft前先加窗,平滑數(shù)據(jù)的邊緣?？梢越档团园甑母叨取Ｅ园晔鞘褂镁匦未爱a(chǎn)生的陡峭的剪切引入的寄生頻率,對于非矩形窗,結(jié)束點(diǎn)衰減的平滑,所以引入較小的寄生頻率。但是,非矩形窗增寬了主瓣,因此降低了頻譜分辨率。函數(shù)periodogram允許指定對數(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>;Hhamm=spectrum.periodogram<'Hamming'>;psd<Hhamm,xn,'Fs',fs,'NFFT',1024>;事實(shí)上加Hamming窗后信號的主瓣大約是矩形窗主瓣的2倍。對固定長度信號,Hamming窗能達(dá)到的譜估計(jì)分辨率大約是矩形窗分辨率的一半。這種沖突可以在某種程度上被變化窗所解決,例如Kaiser窗。非矩形窗會(huì)影響信號的功率,因?yàn)橐恍┎蓸颖幌魅趿?。為了解決這個(gè)問題函數(shù)periodogram將窗歸一化,有平均單位功率。這樣的窗不影響信號的平均功率。修正周期圖法估計(jì)的PSD是其中U是窗歸一化常數(shù)假如U保證估計(jì)是漸進(jìn)無偏的。Welch法包括：將數(shù)據(jù)序列劃分為不同的段〔可以有重疊,對每段進(jìn)行改進(jìn)周期圖法估計(jì),再平均。用spectrum.welch對象,或pwelch函數(shù)。默認(rèn)情況下數(shù)據(jù)劃分為4段,50%重疊,應(yīng)用Hamming窗。取平均的目的是減小方差,重疊會(huì)引入冗余但是加Hamming窗可以部分消除這些冗余,因?yàn)榇敖o邊緣數(shù)據(jù)的權(quán)重比較小。數(shù)據(jù)段的縮短和非矩形窗的使用使得頻譜分辨率下降。下面的例子展示W(wǎng)elch法的折衷。首先用周期圖法估計(jì)一個(gè)小信噪比下信號的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正弦信號已經(jīng)無法識別。Hs=spectrum.welch<'rectangular',150,50>;psd<Hs,xn,'Fs',fs,'NFFT',512>可以看出兩個(gè)信號峰,但是如果進(jìn)一步削減方差,主瓣增寬也使得信號不可識別。Hs=spectrum.welch<'rectangular',100,75>;psd<Hs,xn,'Fs',fs,'NFFT',512>;Welch法的偏差其中是分段數(shù)據(jù)的長度,是窗歸一化常數(shù)。對一定長度的數(shù)據(jù),Welch法估計(jì)的偏差會(huì)大于周期圖法,因?yàn)榉讲畋容^難以量化,因?yàn)樗头侄伍L以及實(shí)用的窗都有關(guān)系,但是總的說方差反比于使用的段數(shù)。MultitaperMethod多椎體法周期圖法估計(jì)可以用濾波器組來表示。L個(gè)帶通濾波器對信號進(jìn)行濾波,每個(gè)濾波器的3dB帶寬是。所有濾波器的幅度響應(yīng)相似于矩形窗的幅度響應(yīng)。周期圖估計(jì)就是對每個(gè)濾波器輸出信號功率的計(jì)算,僅僅使用輸出信號的一個(gè)采樣點(diǎn)計(jì)算輸出信號功率,而且假設(shè)的PSD在每個(gè)濾波器的頻帶上是常數(shù)。信號長度增加,帶通濾波器的帶寬就在減少,近似度就更好。但是有兩個(gè)原因?qū)_度有影響：1矩形窗對應(yīng)的帶通濾波器性能很差2每個(gè)帶通濾波器輸出信號功率的計(jì)算僅僅使用一個(gè)采樣點(diǎn),這使得估計(jì)很粗糙。Welch法也可以用濾波器組給出相似的解釋。在Welch法中使用了多個(gè)點(diǎn)來計(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è)多椎體被用來形成估計(jì)。這就意味著,隨著NW的提高,會(huì)有越來越多的功率譜估計(jì)值,估計(jì)方差會(huì)越來越小。然而,每個(gè)多椎體的帶寬仍然正比于NW,因而NE提高,每個(gè)估計(jì)會(huì)存在更大的泄露,從而整體估計(jì)會(huì)更加呈現(xiàn)有偏。對每一組數(shù)據(jù),總有一個(gè)NW值能在估計(jì)偏差和方差見獲得最好的折中。信號處理工具箱中實(shí)現(xiàn)MTM方法的函數(shù)是pmtm而實(shí)現(xiàn)該方法的對象是spectrum.mtm。下面使用spectrum.mtm來計(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>通過降低時(shí)間-帶寬積,能夠提高分辨率。Hs2=spectrum.mtm<3/2,'adapt'>;psd<Hs2,xn,'Fs',fs,'NFFT',1024>注意到兩個(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ì)算復(fù)雜度更高,這是計(jì)算離散扁平類球體序列的代價(jià)。對于長數(shù)據(jù)序列〔10000點(diǎn)以上,通常計(jì)算一次DPSS序列并將其存為MAT文件更加實(shí)用。Matlab在dpss.mat中提供了dpsssave、dpssload、dpssdir和dpssclear供使用?；プV密度函數(shù)PSD是互譜密度〔CPSD函數(shù)的一個(gè)特例,CPSD由兩個(gè)信號xn、yn如下定義：如同互相關(guān)與協(xié)方差的例子,工具箱估計(jì)PSD和CPSD是因?yàn)樾盘栭L度有限。為了使用Welch方法估計(jì)相隔等長信號x和y的互功率譜密度,cpsd函數(shù)通過將x的FFT和y的FFT再共軛之后相乘的方式得到周期圖。與實(shí)值PSD不同,CPSD是個(gè)復(fù)數(shù)函數(shù)。cpsd如同pwelch函數(shù)一樣處理信號的分段和加窗問題：Sxy=cpsd<x,y,nwin,noverlap,nfft,fs>傳輸函數(shù)估計(jì)Welch方法的一個(gè)應(yīng)用是非參數(shù)系統(tǒng)的識別。假設(shè)H是一個(gè)線性時(shí)不變系統(tǒng),x<n>和y<n>是H的輸入和輸出。則x<n>的功率譜就與x<n>和y<n>的CPSD通過如下方式相關(guān)聯(lián)：x<n>和y<n>的一個(gè)傳輸函數(shù)是：該方法同時(shí)估計(jì)出幅度和相位信息。tfestimate函數(shù)使用Welch方法計(jì)算CPSD和功率譜,然后得到他們的商作為傳輸函數(shù)的估計(jì)值。tfestimate函數(shù)使用方法和cpsd相同：將信號x<n>通過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>'>;相干函數(shù)兩個(gè)信號幅度平方相干性如下所示：該商是一個(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>如果輸入序列長度nfft,窗長度window,一個(gè)窗中重疊的數(shù)據(jù)點(diǎn)為numoverlap,這樣的話mscohere只對一個(gè)樣本操作,函數(shù)返回全1。這是因?yàn)橄喔珊瘮?shù)對線性獨(dú)立數(shù)據(jù)值為1ParametricMethods參數(shù)法參數(shù)法在信號長度較短時(shí)能夠獲得比非參數(shù)法更高的分辨率。這類方法使用不同的方式來估計(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過程。正是由于這個(gè)原因,這一類方法被稱作AR方法。AR方法便于描述譜呈現(xiàn)尖峰的數(shù)據(jù),即PSD在某些頻點(diǎn)特別大。在很多實(shí)際應(yīng)用中〔如語音信號數(shù)據(jù)都具有帶尖峰的譜,所以AR模型通常會(huì)很有用。另外,AR模型具有相對易于求解的系統(tǒng)線性方程。信號處理工具箱提供了下列AR譜估計(jì)方法：Yule-WalkerARmethod<autocorrelationmethod>BurgmethodCovariancemethodModifiedcovariancemethod所有的AR方法都會(huì)給出如下表示的PSD估計(jì)：不同的AR方法估計(jì)AR參數(shù)ap<k>稍有不同,從而得到不一樣的PSD估計(jì)。下表對各種AR方法做了一個(gè)總結(jié)：Burg協(xié)方差修正協(xié)方差Yule-Walker特點(diǎn)對數(shù)據(jù)不加窗；在最小二乘意義上最小化前向后向預(yù)測誤差,限定AR系數(shù)以滿足L-D遞歸；對數(shù)據(jù)不加窗；在最小二乘意義上最小化前向后向預(yù)測誤差；對數(shù)據(jù)不加窗；在最小二乘意義上最小化前向后向預(yù)測誤差；對數(shù)據(jù)加窗；在最小二乘意義上最小化前向預(yù)測誤差〔也叫自相關(guān)法；優(yōu)點(diǎn)對短數(shù)據(jù)具有高分辨率；模型總是穩(wěn)定；對短數(shù)據(jù)比Y-W有更好分辨率〔估計(jì)更準(zhǔn)確；能夠從包含p或更多純正弦信號的數(shù)據(jù)中提取頻率；對短數(shù)據(jù)具有高分辨率；能夠從包含p或更多純正弦信號的數(shù)據(jù)中提取頻率；沒有譜線分裂；對大數(shù)據(jù)性能與其他相當(dāng)；模型總是穩(wěn)定；缺點(diǎn)峰值位置高度依賴于初始相位；在正弦信號包含噪聲或階數(shù)很高時(shí)可能出現(xiàn)譜線分裂；正弦信號估計(jì)頻偏；模型可能不穩(wěn)定；正弦信號估計(jì)頻偏；模型可能不穩(wěn)定；峰值位置高度依賴于初始相位；正弦信號估計(jì)較小頻偏；對于短數(shù)據(jù)性能不高；正弦信號估計(jì)頻偏；非奇異條件階數(shù)必須不大于輸入幀尺寸一半；階數(shù)必須不大于輸入幀尺寸的三分之二；由于估計(jì)有偏,自相關(guān)矩陣需要確保正定；Yule-Walker法Yule-WalkerAR法通過計(jì)算信號自相關(guān)函數(shù)的有偏估計(jì)、求解前向預(yù)測誤差的最小二乘最小化來獲得AR參數(shù)。這就得出了Yule-Walker等式。Yule-WalkerAR法結(jié)果與最大熵估計(jì)器結(jié)果一致。更多信息參考item[2]intheSelectedBibliography。由于自相關(guān)函數(shù)的有偏估計(jì)的使用,確保了上述自相關(guān)矩陣正定。因此,矩陣可逆且方程一定有解。另外,這樣計(jì)算的AR參數(shù)總會(huì)產(chǎn)生一個(gè)穩(wěn)定的全極點(diǎn)模型。Yule-Walker方程通過Levinson算法可以高效的求解。工具箱中的對象spectrum.yulear和函數(shù)pyulear實(shí)現(xiàn)了Tule-Walker方法。下例比較了一個(gè)語音信號通過Welch法和Yule-Walker法的譜：loadmtlbHwelch=spectrum.welch<'hamming',256,50>;psd<Hwelch,mtlb,'Fs',Fs,'NFFT',1024>Hyulear=spectrum.yulear<14>;psd<Hyulear,mtlb,'Fs',Fs,'NFFT',1024>Yule-WalkerAR法的譜比周期圖法更加平滑。這是因?yàn)槠鋬?nèi)在的簡單全極點(diǎn)模型的緣故。Burg法BurgAR法譜估計(jì)是基于最小化前向后向預(yù)測誤差的同時(shí)滿足Levinson-Durbin遞歸〔參考Marple[3],Chapter7,Proakis[6],Section。對比與其它的AR估計(jì)方法,Burg法避免了對自相關(guān)函數(shù)的計(jì)算,改而直接估計(jì)反射系數(shù)。Burg法最首要的優(yōu)勢在于解決含有低噪聲的間隔緊密的正弦信號,并且對短數(shù)據(jù)的估計(jì),在這種情況下AR功率譜密度估計(jì)非常逼近與真值。另外,Burg法確保產(chǎn)生一個(gè)穩(wěn)定AR模型,并且能高效計(jì)算。Burg法的精度在階數(shù)高、數(shù)據(jù)記錄長、信噪比高〔這會(huì)導(dǎo)致線分裂、或者在譜估計(jì)中產(chǎn)生無關(guān)峰的情況下較低。Burg法計(jì)算的譜密度估計(jì)也易受噪聲正弦信號初始相位導(dǎo)致的頻率偏移〔相對于真實(shí)頻率影響。這一效應(yīng)在分析短數(shù)據(jù)序列時(shí)會(huì)被放大。工具箱中的spectrum.burg對象和pburg函數(shù)實(shí)現(xiàn)了Burg法。比較下對于語音信號通過Burg法和Yule-Walker法得到的譜,在較長信號數(shù)據(jù)的情況下它們非常相似。loadmtlbHburg=spectrum.burg<14>;%14thordermodelpsd<Hburg,mtlb<1:512>,'Fs',Fs,'NFFT',1024>Hyulear=spectrum.yulear<14>;%14thordermodelpsd<Hyulear,mtlb<1:512>,'Fs',Fs,'NFFT',1024>比較受噪聲干擾的信號的譜,分別使用Burg法和Welch法計(jì)算：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>>;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ù)的降低,由于正弦信號初始相位帶來的頻率偏移逐漸明顯。Covariance&ModifiedCovari