#基于Matlab的相關(guān)頻譜分析程序教程

1、精選優(yōu)質(zhì)文檔-傾情為你奉上Matlab 信號處理工具箱 譜估計專題 頻譜分析Spectral estimation（譜估計）的目標(biāo)是基于一個有限的數(shù)據(jù)集合描述一個信號的功率（在頻率上的）分布。功率譜估計在很多場合下都是有用的，包括對寬帶噪聲湮沒下的信號的檢測。從數(shù)學(xué)上看，一個平穩(wěn)隨機(jī)過程的power spectrum（功率譜）和correlation sequence（相關(guān)序列）通過discrete-time Fourier transform（離散時間傅立葉變換）構(gòu)成聯(lián)系。從normalized frequency（歸一化角頻率）角度看，有下式注：，其中。其matlab近似為X=fft(x,

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

3、的單位是瓦特，正如平均功率所期望的那樣。對實(shí)信號，PSD是關(guān)于直流信號對稱的，所以的就足夠完整的描述PSD了。然而要獲得整個Nyquist間隔上的平均功率，有必要引入單邊PSD的概念：信號在頻帶上的平均功率可以用單邊PSD求出頻譜估計方法Matlab 信號處理工具箱提供了三種方法Nonparametric methods（非參量類方法）PSD直接從信號本身估計出來。最簡單的就是periodogram（周期圖法），一種改進(jìn)的周期圖法是Welch's method。更現(xiàn)代的一種方法是multitaper method（多椎體法）。Parametric methods （參量類方法）這類方法

4、是假設(shè)信號是一個由白噪聲驅(qū)動的線性系統(tǒng)的輸出。這類方法的例子是Yule-Walker autoregressive (AR) method和Burg method。這些方法先估計假設(shè)的產(chǎn)生信號的線性系統(tǒng)的參數(shù)。這些方法想要對可用數(shù)據(jù)相對較少的情況產(chǎn)生優(yōu)于傳統(tǒng)非參數(shù)方法的結(jié)果。Subspace methods （子空間類）又稱為high-resolution methods（高分辨率法）或者super-resolution methods（超分辨率方法）基于對自相關(guān)矩陣的特征分析或者特征值分解產(chǎn)生信號的頻率分量。代表方法有multiple signal classification (MUSIC

5、) method或eigenvector (EV) method。這類方法對線譜（正弦信號的譜）最合適，對檢測噪聲下的正弦信號很有效，特別是低信噪比的情況。方法描述函數(shù)周期圖PSD 估計spectrum.periodogram, periodogramWelch重疊，加窗的信號段的平均周期圖spectrum.welch, pwelch, cpsd, tfestimate, mscohere多椎體多個正交窗（稱為錐）的組合做譜估計spectrum.mtm, pmtmYule-Walker AR時間序列的估計的自相關(guān)函數(shù)計算自回歸（AR）譜估計spectrum.yulear, pyulearBur

6、g通過最小化線性預(yù)測誤差計算自回歸（AR）譜估計spectrum.burg, pburgCovariance（協(xié)方差）通過最小化前向預(yù)測誤差做時間序列的自回歸（AR）譜估計spectrum.cov, pcov修正協(xié)方差通過最小化前向及后向預(yù)測誤差做時間序列的自回歸（AR）譜估計spectrum.mcov, pmcovMUSIC多重信號分類spectrum.music, pmusic特征向量法虛譜估計spectrum.eigenvector, peigNonparametric Methods非參數(shù)法下面討論periodogram, modified periodogram, Welch, 和

7、multitaper法。同時也討論CPSD函數(shù)，傳輸函數(shù)估計和相關(guān)函數(shù)。Periodogram周期圖法一個估計功率譜的簡單方法是直接求隨機(jī)過程抽樣的DFT，然后取結(jié)果的幅度的平方。這樣的方法叫做周期圖法。一個長L的信號的PSD的周期圖估計是注：這里運(yùn)用的是matlab里面的fft的定義不帶歸一化系數(shù)Matlab FFT函數(shù)未做歸一化，所以要除以L其中實(shí)際對的計算可以只在有限的頻率點(diǎn)上執(zhí)行并且使用FFT。實(shí)踐上大多數(shù)周期圖法的使用都計算N點(diǎn)PSD估計，其中選擇N是大于L的下一個2的冪次是明智的，要計算我們直接對補(bǔ)零到長度為N。假如L>N，在計算前，我們必須繞回模N。作為一個例子，考慮下面1

8、001元素信號，它包含了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 frequencies (column vector)xn = A*sin(2*pi*f*t) + 0.1*randn(size(t);注意：最后三行表明了一個方便的表示正弦之和的方法，它等價于：xn = sin(2*pi*

9、150*t) + 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&#

10、39;);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周期圖性能下面從四個角度討論周期圖法估計的性能：泄漏，分辨率，偏差和方差。頻譜泄漏考慮有限長信號，把它表示成無限長序列乘以一個有限長矩形窗的乘積的形式經(jīng)常很有用：因?yàn)闀r域的乘積等效于頻域的卷積，所以上式的傅立葉變換是前文中導(dǎo)出的表達(dá)式說明卷積對周期圖有影響。正弦數(shù)據(jù)的卷積影響最容易理解

11、。假設(shè)是M個復(fù)正弦的和其頻譜是對一個有限長序列，就變成了所以在有限長信號的頻譜中，Dirac函數(shù)被替換成了形式為的項(xiàng)，該項(xiàng)對應(yīng)于矩形窗的中心在的頻率響應(yīng)。一個矩形窗的頻率響應(yīng)形狀是一個sinc信號，如下所示該圖顯示了一個主瓣和若干旁瓣，最大旁瓣大約在主瓣下方13.5dB處。這些旁瓣說明了頻譜泄漏效應(yīng)。無限長信號的功率嚴(yán)格的集中在離散頻率點(diǎn)處，而有限長信號在離散頻率點(diǎn)附近有連續(xù)的功率。因?yàn)榫匦未霸蕉?，它的頻率響應(yīng)對Dirac沖擊的近似性越差，所以數(shù)據(jù)越短它的頻譜泄漏越明顯?？紤]下面的100個采樣的序列randn('state',0)fs = 1000; % Sampling fr

12、equencyt = (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.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',1024)注意到頻譜泄露只視數(shù)據(jù)長度而定。周期圖確實(shí)只對有限數(shù)據(jù)樣本進(jìn)行計算，但是這和頻譜泄露無關(guān)。分辨率分辨率指的是區(qū)分頻

13、譜特征的能力，是分析譜估計性能的關(guān)鍵概念。要區(qū)分兩個在頻率上離得很近的正弦，要求兩個頻率差大于任何一個信號泄漏頻譜的主瓣寬度。主瓣寬度定義為主瓣上峰值功率一半的點(diǎn)間的距離（3dB帶寬）。該寬度近似等于兩個頻率為的正弦信號，可分辨條件是上例中頻率間隔10Hz，數(shù)據(jù)長度要大于100抽才能使得周期圖中兩個頻率可分辨。下圖是只有67個數(shù)據(jù)長度的情況randn('state',0)fs = 1000; % Sampling frequencyt = (0:fs/15)./fs; % 67 samplesA = 1 2; % Sinusoid amplitudesf = 150;140;

14、% Sinusoid frequenciesxn = 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)信噪比低的時候，譜特征的分辨更難，而且周期圖上會出現(xiàn)一些噪聲的偽像，如下所示randn('state',0)fs = 1000; % Sampling frequencyt = (0:fs/10)./fs; % One-tenth of a secon

15、d 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',1024)估計偏差周期圖是對PSD的有偏估計。期望值可以是該式和頻譜泄漏中的式相似，除了這里的表達(dá)式用的是平均功率而不是幅度。這暗示了周期圖產(chǎn)生的估計對應(yīng)于一個有泄漏的PSD而非真正的PSD。注意本質(zhì)上是一個三角Bartlet

16、t窗（事實(shí)是兩個矩形脈沖的卷積是三角脈沖。）這導(dǎo)致了最大旁瓣峰值比主瓣峰值低27dB，大致是非平方矩形窗的2倍。周期圖估計是漸進(jìn)無偏的。這從早期的一個觀察結(jié)果可以明顯看出，隨著記錄數(shù)據(jù)趨于無窮大，矩形窗對頻譜對Dirac函數(shù)的近似也就越來越好。然而在某些情況下，周期圖法估計很差勁即使數(shù)據(jù)夠長，這是因?yàn)橹芷趫D法的方差，如下所述。周期圖法的方差L趨于無窮大，方差也不趨于0。用統(tǒng)計學(xué)術(shù)語講，該估計不是無偏估計。然而周期圖在信噪比大的時候仍然是有用的譜估計器，特別是數(shù)據(jù)夠長。Modified Periodogram修正周期圖法在fft前先加窗，平滑數(shù)據(jù)的邊緣?？梢越档团园甑母叨?。旁瓣是使用矩形窗產(chǎn)生的

17、陡峭的剪切引入的寄生頻率，對于非矩形窗，結(jié)束點(diǎn)衰減的平滑，所以引入較小的寄生頻率。但是，非矩形窗增寬了主瓣，因此降低了頻譜分辨率。函數(shù)periodogram允許指定對數(shù)據(jù)加的窗，例如默認(rèn)的矩形窗和Hamming窗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

18、*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á)到的譜估計分辨率大約是矩形窗分辨率的一半。這種沖突可以在某種程度上被變化窗所解決，例如Kaiser

19、窗。非矩形窗會影響信號的功率，因?yàn)橐恍┎蓸颖幌魅趿?。為了解決這個問題函數(shù)periodogram將窗歸一化，有平均單位功率。這樣的窗不影響信號的平均功率。修正周期圖法估計的PSD是其中U是窗歸一化常數(shù)假如U保證估計是漸進(jìn)無偏的。Welch法包括：將數(shù)據(jù)序列劃分為不同的段（可以有重疊），對每段進(jìn)行改進(jìn)周期圖法估計，再平均。用spectrum.welch對象，或pwelch函數(shù)。默認(rèn)情況下數(shù)據(jù)劃分為4段，50%重疊，使用Hamming窗。取平均的目的是減小方差，重疊會引入冗余但是加Hamming窗可以部分消除這些冗余，因?yàn)榇敖o邊緣數(shù)據(jù)的權(quán)重比較小。數(shù)據(jù)段的縮短和非矩形窗的使用使得頻譜分辨率下降。下面

20、的例子展示W(wǎng)elch法的折衷。首先用周期圖法估計一個小信噪比下信號的PSD：randn('state',1)fs = 1000; % Sampling frequencyt = (0:0.3*fs)./fs; % 301 samplesA = 2 8; % Sinusoid amplitudes (row vector)f = 150;140; % Sinusoid frequencies (column vector)xn = A*sin(2*pi*f*t) + 5*randn(size(t);Hs = spectrum.periodogram('rectangula

21、r')psd(Hs,xn,'Fs',fs,'NFFT',1024);可以看出由于噪聲太大，150Hz正弦信號已經(jīng)無法識別。Hs = spectrum.welch('rectangular',150,50);psd(Hs,xn,'Fs',fs,'NFFT',512)可以看出兩個信號峰，但是如果進(jìn)一步削減方差，主瓣增寬也使得信號不可識別。Hs = spectrum.welch('rectangular',100,75);psd(Hs,xn,'Fs',fs,'NFFT

22、9;,512);Welch法的偏差其中是分段數(shù)據(jù)的長度，是窗歸一化常數(shù)。對一定長度的數(shù)據(jù)，Welch法估計的偏差會大于周期圖法，因?yàn)榉讲畋容^難以量化，因?yàn)樗头侄伍L以及實(shí)用的窗都有關(guān)系，但是總的說方差反比于使用的段數(shù)。Multitaper Method多椎體法周期圖法估計可以用濾波器組來表示。L個帶通濾波器對信號進(jìn)行濾波，每個濾波器的3dB帶寬是。所有濾波器的幅度響應(yīng)相似于矩形窗的幅度響應(yīng)。周期圖估計就是對每個濾波器輸出信號功率的計算，僅僅使用輸出信號的一個采樣點(diǎn)計算輸出信號功率，而且假設(shè)的PSD在每個濾波器的頻帶上是常數(shù)。信號長度增加，帶通濾波器的帶寬就在減少，近似度就更好。但是有兩個原因?qū)?/p>

23、精確度有影響：1矩形窗對應(yīng)的帶通濾波器性能很差2每個帶通濾波器輸出信號功率的計算僅僅使用一個采樣點(diǎn)，這使得估計很粗糙。Welch法也可以用濾波器組給出相似的解釋。在Welch法中使用了多個點(diǎn)來計算輸出功率，降低了估計的方差。另一方面每個帶通濾波器的帶寬增大了，分辨率下降了。Thompson的多椎體法（MTM）構(gòu)建在上述結(jié)論之上，提供更優(yōu)的PSD估計。MTM方法沒有使用帶通濾波器（它們本質(zhì)上是矩形窗，如同周期圖法中一樣），而是使用一組最優(yōu)濾波器計算估計值。這些最優(yōu)FIR濾波器是由一組被叫做離散扁平類球體序列（DPSS，也叫做Slepian序列）得到的。除此之外，MTM方法提供了一個時間-帶寬參數(shù)

24、，有了它能在估計方差和分辨率之間進(jìn)行平衡。該參數(shù)由時間-帶寬乘積得到，NW，同時它直接和譜估計的多椎體數(shù)有關(guān)?？傆?*NW-1個多椎體被用來形成估計。這就意味著，隨著NW的提高，會有越來越多的功率譜估計值，估計方差會越來越小。然而，每個多椎體的帶寬仍然正比于NW，因而NE提高，每個估計會存在更大的泄露，從而整體估計會更加呈現(xiàn)有偏。對每一組數(shù)據(jù)，總有一個NW值能在估計偏差和方差見獲得最好的折中。信號處理工具箱中實(shí)現(xiàn)MTM方法的函數(shù)是pmtm而實(shí)現(xiàn)該方法的對象是spectrum.mtm。下面使用spectrum.mtm來計算前一個例子中的PSD：randn('state',0)fs

25、 = 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 = spectrum.mtm(4,'adapt');psd(Hs1,xn,'Fs',fs,'NFFT',1024)通過降低時間-帶寬積，能夠提高分辨率。Hs2 = spectrum.mt

26、m(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,xn,'Fs',fs,'NFFT',1024);Pow2 = avgpower(Hs2p)Pow2 = 2.4927這中方法相比Weich方法計算復(fù)雜度更高，這是計算離散扁平類球體序列的

27、代價。對于長數(shù)據(jù)序列（10000點(diǎn)以上），通常計算一次DPSS序列并將其存為MAT文件更加實(shí)用。Matlab在dpss.mat中提供了dpsssave、dpssload、dpssdir和 dpssclear供使用?；プV密度函數(shù)PSD是互譜密度（CPSD）函數(shù)的一個特例，CPSD由兩個信號xn、yn如下定義：如同互相關(guān)和協(xié)方差的例子，工具箱估計PSD和CPSD是因?yàn)樾盘栭L度有限。為了使用Welch方法估計相隔等長信號x和y的互功率譜密度，cpsd函數(shù)通過將x的FFT和y的FFT再共軛之后相乘的方式得到周期圖。和實(shí)值PSD不同，CPSD是個復(fù)數(shù)函數(shù)。cpsd如同pwelch函數(shù)一樣處理信號的分段和

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

29、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('Frequen

30、cy (Hz)');相干函數(shù)兩個信號幅度平方相干性如下所示：該商是一個0到1之間的實(shí)數(shù)，表征了x(n)和y(n)之間的相干性。mscohere函數(shù)輸入兩個序列x和y，計算其功率譜和CPSD，返回CPSD幅度平方和兩個功率譜乘積的商。函數(shù)的選項(xiàng)和操作和cpsd和tfestimate相類似。x和濾波器輸出y的相干函數(shù)如下：mscohere(xn, yn, 256, 128, 256, fs)如果輸入序列長度nfft，窗長度window，一個窗中重疊的數(shù)據(jù)點(diǎn)為numoverlap，這樣的話mscohere只對一個樣本操作，函數(shù)返回全1。這是因?yàn)橄喔珊瘮?shù)對線性獨(dú)立數(shù)據(jù)值為1Parametric

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

32、提供了下列AR譜估計方法：l Yule-Walkerl AR method (autocorrelation method)Burgl methodCovariancel methodModifiedl covariance method所有的AR方法都會給出如下表示的PSD估計：不同的AR方法估計AR參數(shù)ap(k)稍有不同，從而得到不一樣的PSD估計。下表對各種AR方法做了一個總結(jié)：Burg協(xié)方差修正協(xié)方差Yule-Walker特點(diǎn)對數(shù)據(jù)不加窗；在最小二乘意義上最小化前向后向預(yù)測誤差，限定AR系數(shù)以滿足L-D遞歸；對數(shù)據(jù)不加窗；在最小二乘意義上最小化前向后向預(yù)測誤差；對數(shù)據(jù)不加窗；在最小二乘

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

34、大于輸入幀尺寸一半；階數(shù)必須不大于輸入幀尺寸的三分之二；由于估計有偏，自相關(guān)矩陣需要確保正定；Yule-Walker 法Yule-Walker AR法通過計算信號自相關(guān)函數(shù)的有偏估計、求解前向預(yù)測誤差的最小二乘最小化來獲得AR參數(shù)。這就得出了Yule-Walker等式。Yule-Walker AR法結(jié)果和最大熵估計器結(jié)果一致。更多信息參考item 2 in the Selected Bibliography。由于自相關(guān)函數(shù)的有偏估計的使用，確保了上述自相關(guān)矩陣正定。因此，矩陣可逆且方程一定有解。另外，這樣計算的AR參數(shù)總會產(chǎn)生一個穩(wěn)定的全極點(diǎn)模型。Yule-Walker方程通過Levinson

35、算法可以高效的求解。工具箱中的對象spectrum.yulear和函數(shù)pyulear實(shí)現(xiàn)了Tule-Walker方法。下例比較了一個語音信號通過Welch法和Yule-Walker法的譜：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)Yule-Walker

36、AR法的譜比周期圖法更加平滑。這是因?yàn)槠鋬?nèi)在的簡單全極點(diǎn)模型的緣故。Burg法Burg AR法譜估計是基于最小化前向后向預(yù)測誤差的同時滿足Levinson-Durbin遞歸（參考Marple3，Chapter 7，Proakis6，Section 12.3.3）。對比和其它的AR估計方法，Burg法避免了對自相關(guān)函數(shù)的計算，改而直接估計反射系數(shù)。Burg法最首要的優(yōu)勢在于解決含有低噪聲的間隔緊密的正弦信號，并且對短數(shù)據(jù)的估計，在這種情況下AR功率譜密度估計非常逼近和真值。另外，Burg法確保產(chǎn)生一個穩(wěn)定AR模型，并且能高效計算。Burg法的精度在階數(shù)高、數(shù)據(jù)記錄長、信噪比高（這會導(dǎo)致線分裂、或

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

38、m.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; % 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法模型階數(shù)的降低，由于正弦信號初始相位帶來的