




版權(quán)說(shuō)明:本文檔由用戶(hù)提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請(qǐng)進(jìn)行舉報(bào)或認(rèn)領(lǐng)
文檔簡(jiǎn)介
1、Matlab 信號(hào)處理工具箱 幫助文檔 譜估計(jì)專(zhuān)題 翻譯:無(wú)名網(wǎng)友 & Lyra頻譜分析Spectral estimation(譜估計(jì))的目標(biāo)是基于一個(gè)有限的數(shù)據(jù)集合描述一個(gè)信號(hào)的功率(在頻率上的)分布。功率譜估計(jì)在很多場(chǎng)合下都是有用的,包括對(duì)寬帶噪聲湮沒(méi)下的信號(hào)的檢測(cè)。從數(shù)學(xué)上看,一個(gè)平穩(wěn)隨機(jī)過(guò)程的power spectrum(功率譜)和correlation sequence(相關(guān)序列)通過(guò)discrete-time Fourier transform(離散時(shí)間傅立葉變換)構(gòu)成聯(lián)系。從normalized frequency(歸一化角頻率)角度看,有下式注:,其中。其matlab近
2、似為X=fft(x,N)/sqrt(N),在下文中就是指matlab fft函數(shù)的計(jì)算結(jié)果了使用關(guān)系可以寫(xiě)成物理頻率的函數(shù),其中是采樣頻率相關(guān)序列可以從功率譜用IDFT變換求得:序列在整個(gè)Nyquist間隔上的平均功率可以表示為上式中的以及被定義為平穩(wěn)隨機(jī)信號(hào)的power spectral density (PSD)(功率譜密度)一個(gè)信號(hào)在頻帶上的平均功率可以通過(guò)對(duì)PSD在頻帶上積分求出從上式中可以看出是一個(gè)信號(hào)在一個(gè)無(wú)窮小頻帶上的功率濃度,這也是為什么它叫做功率譜密度。PSD的單位是功率(e.g 瓦特)每單位頻率。在的情況下,這是瓦特/弧度/抽或只是瓦特/弧度。在的情況下單位是瓦特/赫茲。P
3、SD對(duì)頻率的積分得到的單位是瓦特,正如平均功率所期望的那樣。對(duì)實(shí)信號(hào),PSD是關(guān)于直流信號(hào)對(duì)稱(chēng)的,所以的就足夠完整的描述PSD了。然而要獲得整個(gè)Nyquist間隔上的平均功率,有必要引入單邊PSD的概念:信號(hào)在頻帶上的平均功率可以用單邊PSD求出頻譜估計(jì)方法cpsdCpsdMatlab 信號(hào)處理工具箱提供了三種方法Nonparametric methods(非參量類(lèi)方法)PSD直接從信號(hào)本身估計(jì)出來(lái)。最簡(jiǎn)單的就是periodogram(周期圖法),一種改進(jìn)的周期圖法是Welch's method。更現(xiàn)代的一種方法是multitaper method(多椎體法)。Parametric m
4、ethods (參量類(lèi)方法)這類(lèi)方法是假設(shè)信號(hào)是一個(gè)由白噪聲驅(qū)動(dòng)的線(xiàn)性系統(tǒng)的輸出。這類(lèi)方法的例子是Yule-Walker autoregressive (AR) method和Burg method。這些方法先估計(jì)假設(shè)的產(chǎn)生信號(hào)的線(xiàn)性系統(tǒng)的參數(shù)。這些方法想要對(duì)可用數(shù)據(jù)相對(duì)較少的情況產(chǎn)生優(yōu)于傳統(tǒng)非參數(shù)方法的結(jié)果。Subspace methods (子空間類(lèi))又稱(chēng)為high-resolution methods(高分辨率法)或者super-resolution methods(超分辨率方法)基于對(duì)自相關(guān)矩陣的特征分析或者特征值分解產(chǎn)生信號(hào)的頻率分量。代表方法有multiple signal cla
5、ssification (MUSIC) method或eigenvector (EV) method。這類(lèi)方法對(duì)線(xiàn)譜(正弦信號(hào)的譜)最合適,對(duì)檢測(cè)噪聲下的正弦信號(hào)很有效,特別是低信噪比的情況。方法描述函數(shù)周期圖PSD 估計(jì)spectrum.periodogram, periodogramWelch重疊,加窗的信號(hào)段的平均周期圖spectrum.welch, pwelch, cpsd, tfestimate, mscohere多椎體多個(gè)正交窗(稱(chēng)為錐)的組合做譜估計(jì)spectrum.mtm, pmtmYule-Walker AR時(shí)間序列的估計(jì)的自相關(guān)函數(shù)計(jì)算自回歸(AR)譜估計(jì)spectrum.
6、yulear, pyulearBurg通過(guò)最小化線(xiàn)性預(yù)測(cè)誤差計(jì)算自回歸(AR)譜估計(jì)spectrum.burg, pburgCovariance(協(xié)方差)通過(guò)最小化前向預(yù)測(cè)誤差做時(shí)間序列的自回歸(AR)譜估計(jì)spectrum.cov, pcov修正協(xié)方差通過(guò)最小化前向及后向預(yù)測(cè)誤差做時(shí)間序列的自回歸(AR)譜估計(jì)spectrum.mcov, pmcovMUSIC多重信號(hào)分類(lèi)spectrum.music, pmusic特征向量法虛譜估計(jì)spectrum.eigenvector, peigNonparametric Methods非參數(shù)法下面討論periodogram, modified peri
7、odogram, 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)的PSD的周期圖估計(jì)是注:這里運(yùn)用的是matlab里面的fft的定義不帶歸一化系數(shù)Matlab FFT函數(shù)未做歸一化,所以要除以L(fǎng)其中實(shí)際對(duì)的計(jì)算可以只在有限的頻率點(diǎn)上執(zhí)行并且使用FFT。實(shí)踐上大多數(shù)周期圖法的應(yīng)用都計(jì)算N點(diǎn)PSD估計(jì),其中選擇N是大于L的下一個(gè)2的冪次是明智的,要計(jì)算我們直接對(duì)補(bǔ)零到長(zhǎng)度為N。假如L>N,在計(jì)算前,我們必
8、須繞回模N。作為一個(gè)例子,考慮下面1001元素信號(hào),它包含了2個(gè)正弦信號(hào)和噪聲r(shí)andn('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);注意:最后三行表明了一個(gè)方便的表示正弦之和的方法,它
9、等價(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ì)象(periodogram object)來(lái)計(jì)算Hs = spectrum.periodogram('Hamming');估計(jì)的圖形可以用psd函數(shù)顯示。psd(Hs,xn,'Fs',fs,'NFFT',1024,'SpectrumType','twosided')平均功率通過(guò)用下述求和去近似積分 求得Pxx,F = psd(Hs,xn,
10、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),把它表示成無(wú)限長(zhǎng)序列乘以一個(gè)有限長(zhǎng)矩形窗的乘積的形式經(jīng)常很有用:因?yàn)闀r(shí)域的乘積等效于頻域的卷積,所以上式的傅立葉變換是前文中導(dǎo)出的表達(dá)式說(shuō)明卷積對(duì)周期圖
11、有影響。正弦數(shù)據(jù)的卷積影響最容易理解。假設(shè)是M個(gè)復(fù)正弦的和其頻譜是對(duì)一個(gè)有限長(zhǎng)序列,就變成了所以在有限長(zhǎng)信號(hào)的頻譜中,Dirac函數(shù)被替換成了形式為的項(xiàng),該項(xiàng)對(duì)應(yīng)于矩形窗的中心在的頻率響應(yīng)。一個(gè)矩形窗的頻率響應(yīng)形狀是一個(gè)sinc信號(hào),如下所示該圖顯示了一個(gè)主瓣和若干旁瓣,最大旁瓣大約在主瓣下方13.5dB處。這些旁瓣說(shuō)明了頻譜泄漏效應(yīng)。無(wú)限長(zhǎng)信號(hào)的功率嚴(yán)格的集中在離散頻率點(diǎn)處,而有限長(zhǎng)信號(hào)在離散頻率點(diǎn)附近有連續(xù)的功率。因?yàn)榫匦未霸蕉?,它的頻率響應(yīng)對(duì)Dirac沖擊的近似性越差,所以數(shù)據(jù)越短它的頻譜泄漏越明顯。考慮下面的100個(gè)采樣的序列randn('state',0)fs = 1
12、000; % 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*f*t) + 0.1*randn(size(t);Hs = spectrum.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',1024)注意到頻譜泄露只視數(shù)據(jù)長(zhǎng)度而定。周期圖確實(shí)只對(duì)有限數(shù)據(jù)樣本進(jìn)行計(jì)算,但是這和頻
13、譜泄露無(wú)關(guān)。分辨率分辨率指的是區(qū)分頻譜特征的能力,是分析譜估計(jì)性能的關(guān)鍵概念。要區(qū)分兩個(gè)在頻率上離得很近的正弦,要求兩個(gè)頻率差大于任何一個(gè)信號(hào)泄漏頻譜的主瓣寬度。主瓣寬度定義為主瓣上峰值功率一半的點(diǎn)間的距離(3dB帶寬)。該寬度近似等于兩個(gè)頻率為的正弦信號(hào),可分辨條件是上例中頻率間隔10Hz,數(shù)據(jù)長(zhǎng)度要大于100抽才能使得周期圖中兩個(gè)頻率可分辨。下圖是只有67個(gè)數(shù)據(jù)長(zhǎng)度的情況randn('state',0)fs = 1000; % Sampling frequencyt = (0:fs/15)./fs; % 67 samplesA = 1 2; % Sinusoid ampli
14、tudesf = 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)上述對(duì)分辨率的討論都是在高信噪比的情況進(jìn)行的,因此沒(méi)有考慮噪聲。當(dāng)信噪比低的時(shí)候,譜特征的分辨更難,而且周期圖上會(huì)出現(xiàn)一些噪聲的偽像,如下所示randn('state',0)fs = 1000; % Sampling frequencyt = (0:fs/10)./fs; % On
15、e-tenth of a second 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)估計(jì)偏差周期圖是對(duì)PSD的有偏估計(jì)。期望值可以是該式和頻譜泄漏中的式相似,除了這里的表達(dá)式用的是平均功率而不是幅度。這暗示了周期圖產(chǎn)生的估計(jì)對(duì)應(yīng)于一個(gè)有泄漏的PSD而非真正的PSD
16、。注意本質(zhì)上是一個(gè)三角Bartlett窗(事實(shí)是兩個(gè)矩形脈沖的卷積是三角脈沖。)這導(dǎo)致了最大旁瓣峰值比主瓣峰值低27dB,大致是非平方矩形窗的2倍。周期圖估計(jì)是漸進(jìn)無(wú)偏的。這從早期的一個(gè)觀察結(jié)果可以明顯看出,隨著記錄數(shù)據(jù)趨于無(wú)窮大,矩形窗對(duì)頻譜對(duì)Dirac函數(shù)的近似也就越來(lái)越好。然而在某些情況下,周期圖法估計(jì)很差勁即使數(shù)據(jù)夠長(zhǎng),這是因?yàn)橹芷趫D法的方差,如下所述。周期圖法的方差L趨于無(wú)窮大,方差也不趨于0。用統(tǒng)計(jì)學(xué)術(shù)語(yǔ)講,該估計(jì)不是無(wú)偏估計(jì)。然而周期圖在信噪比大的時(shí)候仍然是有用的譜估計(jì)器,特別是數(shù)據(jù)夠長(zhǎng)。Modified Periodogram修正周期圖法在fft前先加窗,平滑數(shù)據(jù)的邊緣。可以降
17、低旁瓣的高度。旁瓣是使用矩形窗產(chǎn)生的陡峭的剪切引入的寄生頻率,對(duì)于非矩形窗,結(jié)束點(diǎn)衰減的平滑,所以引入較小的寄生頻率。但是,非矩形窗增寬了主瓣,因此降低了頻譜分辨率。函數(shù)periodogram允許指定對(duì)數(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 frequenc
18、iesxn = 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窗后信號(hào)的主瓣大約是矩形窗主瓣的2倍。對(duì)固定長(zhǎng)度信號(hào),Hamming窗能達(dá)到的譜估計(jì)分辨率大約是矩形窗分辨率的一半。這種沖突可以在某種程
19、度上被變化窗所解決,例如Kaiser窗。非矩形窗會(huì)影響信號(hào)的功率,因?yàn)橐恍┎蓸颖幌魅趿?。為了解決這個(gè)問(wèn)題函數(shù)periodogram將窗歸一化,有平均單位功率。這樣的窗不影響信號(hào)的平均功率。修正周期圖法估計(jì)的PSD是其中U是窗歸一化常數(shù)假如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窗可以部分消除這些冗余,因?yàn)榇敖o邊緣數(shù)據(jù)的權(quán)重比較小。數(shù)據(jù)段的縮短和非
20、矩形窗的使用使得頻譜分辨率下降。下面的例子展示W(wǎng)elch法的折衷。首先用周期圖法估計(jì)一個(gè)小信噪比下信號(hào)的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.periodogr
21、am('rectangular')psd(Hs,xn,'Fs',fs,'NFFT',1024);可以看出由于噪聲太大,150Hz正弦信號(hào)已經(jīng)無(wú)法識(shí)別。Hs = spectrum.welch('rectangular',150,50);psd(Hs,xn,'Fs',fs,'NFFT',512)可以看出兩個(gè)信號(hào)峰,但是如果進(jìn)一步削減方差,主瓣增寬也使得信號(hào)不可識(shí)別。Hs = spectrum.welch('rectangular',100,75);psd(Hs,xn,'Fs
22、9;,fs,'NFFT',512);Welch法的偏差其中是分段數(shù)據(jù)的長(zhǎng)度,是窗歸一化常數(shù)。對(duì)一定長(zhǎng)度的數(shù)據(jù),Welch法估計(jì)的偏差會(huì)大于周期圖法,因?yàn)榉讲畋容^難以量化,因?yàn)樗头侄伍L(zhǎng)以及實(shí)用的窗都有關(guān)系,但是總的說(shuō)方差反比于使用的段數(shù)。Multitaper Method多椎體法周期圖法估計(jì)可以用濾波器組來(lái)表示。L個(gè)帶通濾波器對(duì)信號(hào)進(jìn)行濾波,每個(gè)濾波器的3dB帶寬是。所有濾波器的幅度響應(yīng)相似于矩形窗的幅度響應(yīng)。周期圖估計(jì)就是對(duì)每個(gè)濾波器輸出信號(hào)功率的計(jì)算,僅僅使用輸出信號(hào)的一個(gè)采樣點(diǎn)計(jì)算輸出信號(hào)功率,而且假設(shè)的PSD在每個(gè)濾波器的頻帶上是常數(shù)。信號(hào)長(zhǎng)度增加,帶通濾波器的帶寬就在
23、減少,近似度就更好。但是有兩個(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方法沒(méi)有使用帶通濾波器(它們本質(zhì)上是矩形窗,如同周期圖法中一樣),而是使用一組最優(yōu)濾波器計(jì)算估計(jì)值。這些最優(yōu)FIR濾波器是由一組被叫做離散扁平類(lèi)球體序列(DPSS,也叫做Slepian序列)得到的。除此之外
24、,MTM方法提供了一個(gè)時(shí)間-帶寬參數(shù),有了它能在估計(jì)方差和分辨率之間進(jìn)行平衡。該參數(shù)由時(shí)間-帶寬乘積得到,NW,同時(shí)它直接與譜估計(jì)的多椎體數(shù)有關(guān)??傆?*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ì)偏差和方差見(jiàn)獲得最好的折中。信號(hào)處理工具箱中實(shí)現(xiàn)MTM方法的函數(shù)是pmtm而實(shí)現(xiàn)該方法的對(duì)象是spectrum.mtm。下面使用spectrum.mtm來(lái)計(jì)算前一個(gè)例子中的PSD:randn(
25、39;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);Hs1 = spectrum.mtm(4,'adapt');psd(Hs1,xn,'Fs',fs,'NFFT',1024)通過(guò)降低時(shí)間-帶寬積,能夠提高分辨率
26、。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ù)雜
27、度更高,這是計(jì)算離散扁平類(lèi)球體序列的代價(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如下定義:如同互相關(guān)與協(xié)方差的例子,工具箱估計(jì)PSD和CPSD是因?yàn)樾盘?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如同
28、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è)線(xiàn)性時(shí)不變系統(tǒng),x(n)和y(n)是H的輸入和輸出。則x(n)的功率譜就與x(n)和y(n)的CPSD通過(guò)如下方式相關(guān)聯(lián):x(n)和y(n)的一個(gè)傳輸函數(shù)是:該方法同時(shí)估計(jì)出幅度和相位信息。tfestimate函數(shù)使用Welch方法計(jì)算CPSD和功率譜,然后得到他們的商作為傳輸函數(shù)的估計(jì)值。tfestimate函數(shù)使用方法和cpsd相同:將信號(hào)x(n)通過(guò)FIR濾波器,再畫(huà)出實(shí)際的幅度響應(yīng)和估計(jì)響應(yīng)
29、如下:h = 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'); x
30、label('Frequency (Hz)');相干函數(shù)兩個(gè)信號(hào)幅度平方相干性如下所示:該商是一個(gè)0到1之間的實(shí)數(shù),表征了x(n)和y(n)之間的相干性。mscohere函數(shù)輸入兩個(gè)序列x和y,計(jì)算其功率譜和CPSD,返回CPSD幅度平方與兩個(gè)功率譜乘積的商。函數(shù)的選項(xiàng)和操作與cpsd和tfestimate相類(lèi)似。x和濾波器輸出y的相干函數(shù)如下:mscohere(xn, yn, 256, 128, 256, fs)如果輸入序列長(zhǎng)度nfft,窗長(zhǎng)度window,一個(gè)窗中重疊的數(shù)據(jù)點(diǎn)為numoverlap,這樣的話(huà)mscohere只對(duì)一個(gè)樣本操作,函數(shù)返回全1。這是因?yàn)橄喔珊瘮?shù)對(duì)線(xiàn)
31、性獨(dú)立數(shù)據(jù)值為1Parametric Methods參數(shù)法參數(shù)法在信號(hào)長(zhǎng)度較短時(shí)能夠獲得比非參數(shù)法更高的分辨率。這類(lèi)方法使用不同的方式來(lái)估計(jì)頻譜:不是試圖直接從數(shù)據(jù)中估計(jì)PSD,而是將數(shù)據(jù)建模成一個(gè)由白噪聲驅(qū)動(dòng)的線(xiàn)性系統(tǒng)的輸出,并試圖估計(jì)出該系統(tǒng)的參數(shù)。最常用的線(xiàn)性系統(tǒng)模型是全極點(diǎn)模型,也就是一個(gè)濾波器,它的所有零點(diǎn)都在z平面的原點(diǎn)。這樣一個(gè)濾波器輸入白噪聲后的輸出是一個(gè)自回歸(AR)過(guò)程。正是由于這個(gè)原因,這一類(lèi)方法被稱(chēng)作AR方法。AR方法便于描述譜呈現(xiàn)尖峰的數(shù)據(jù),即PSD在某些頻點(diǎn)特別大。在很多實(shí)際應(yīng)用中(如語(yǔ)音信號(hào))數(shù)據(jù)都具有帶尖峰的譜,所以AR模型通常會(huì)很有用。另外,AR模型具有相對(duì)易
32、于求解的系統(tǒng)線(xiàn)性方程。信號(hào)處理工具箱提供了下列AR譜估計(jì)方法:l Yule-Walkerl AR method (autocorrelation method)Burgl methodCovariancel methodModifiedl covariance method所有的AR方法都會(huì)給出如下表示的PSD估計(jì):不同的AR方法估計(jì)AR參數(shù)ap(k)稍有不同,從而得到不一樣的PSD估計(jì)。下表對(duì)各種AR方法做了一個(gè)總結(jié):Burg協(xié)方差修正協(xié)方差Yule-Walker特點(diǎn)對(duì)數(shù)據(jù)不加窗;在最小二乘意義上最小化前向后向預(yù)測(cè)誤差,限定AR系數(shù)以滿(mǎn)足L-D遞歸;對(duì)數(shù)據(jù)不加窗;在最小二乘意義上最小化前向后
33、向預(yù)測(cè)誤差;對(duì)數(shù)據(jù)不加窗;在最小二乘意義上最小化前向后向預(yù)測(cè)誤差;對(duì)數(shù)據(jù)加窗;在最小二乘意義上最小化前向預(yù)測(cè)誤差(也叫自相關(guān)法);優(yōu)點(diǎn)對(duì)短數(shù)據(jù)具有高分辨率;模型總是穩(wěn)定;對(duì)短數(shù)據(jù)比Y-W有更好分辨率(估計(jì)更準(zhǔn)確);能夠從包含p或更多純正弦信號(hào)的數(shù)據(jù)中提取頻率;對(duì)短數(shù)據(jù)具有高分辨率;能夠從包含p或更多純正弦信號(hào)的數(shù)據(jù)中提取頻率;沒(méi)有譜線(xiàn)分裂;對(duì)大數(shù)據(jù)性能與其他相當(dāng);模型總是穩(wěn)定;缺點(diǎn)峰值位置高度依賴(lài)于初始相位;在正弦信號(hào)包含噪聲或階數(shù)很高時(shí)可能出現(xiàn)譜線(xiàn)分裂;正弦信號(hào)估計(jì)頻偏;模型可能不穩(wěn)定;正弦信號(hào)估計(jì)頻偏;模型可能不穩(wěn)定;峰值位置高度依賴(lài)于初始相位;正弦信號(hào)估計(jì)較小頻偏;對(duì)于短數(shù)據(jù)性能不高;正
34、弦信號(hào)估計(jì)頻偏;非奇異條件階數(shù)必須不大于輸入幀尺寸一半;階數(shù)必須不大于輸入幀尺寸的三分之二;由于估計(jì)有偏,自相關(guān)矩陣需要確保正定;Yule-Walker 法Yule-Walker AR法通過(guò)計(jì)算信號(hào)自相關(guān)函數(shù)的有偏估計(jì)、求解前向預(yù)測(cè)誤差的最小二乘最小化來(lái)獲得AR參數(shù)。這就得出了Yule-Walker等式。Yule-Walker AR法結(jié)果與最大熵估計(jì)器結(jié)果一致。更多信息參考item 2 in the Selected Bibliography。由于自相關(guān)函數(shù)的有偏估計(jì)的使用,確保了上述自相關(guān)矩陣正定。因此,矩陣可逆且方程一定有解。另外,這樣計(jì)算的AR參數(shù)總會(huì)產(chǎn)生一個(gè)穩(wěn)定的全極點(diǎn)模型。Yule-
35、Walker方程通過(guò)Levinson算法可以高效的求解。工具箱中的對(duì)象spectrum.yulear和函數(shù)pyulear實(shí)現(xiàn)了Tule-Walker方法。下例比較了一個(gè)語(yǔ)音信號(hào)通過(guò)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'
36、,1024)Yule-Walker AR法的譜比周期圖法更加平滑。這是因?yàn)槠鋬?nèi)在的簡(jiǎn)單全極點(diǎn)模型的緣故。Burg法Burg AR法譜估計(jì)是基于最小化前向后向預(yù)測(cè)誤差的同時(shí)滿(mǎn)足Levinson-Durbin遞歸(參考Marple3,Chapter 7,Proakis6,Section )。對(duì)比與其它的AR估計(jì)方法,Burg法避免了對(duì)自相關(guān)函數(shù)的計(jì)算,改而直接估計(jì)反射系數(shù)。Burg法最首要的優(yōu)勢(shì)在于解決含有低噪聲的間隔緊密的正弦信號(hào),并且對(duì)短數(shù)據(jù)的估計(jì),在這種情況下AR功率譜密度估計(jì)非常逼近與真值。另外,Burg法確保產(chǎn)生一個(gè)穩(wěn)定AR模型,并且能高效計(jì)算。Burg法的精度在階數(shù)高、數(shù)據(jù)記錄長(zhǎng)、信噪
37、比高(這會(huì)導(dǎo)致線(xiàn)分裂、或者在譜估計(jì)中產(chǎn)生無(wú)關(guān)峰)的情況下較低。Burg法計(jì)算的譜密度估計(jì)也易受噪聲正弦信號(hào)初始相位導(dǎo)致的頻率偏移(相對(duì)于真實(shí)頻率)影響。這一效應(yīng)在分析短數(shù)據(jù)序列時(shí)會(huì)被放大。工具箱中的spectrum.burg對(duì)象和pburg函數(shù)實(shí)現(xiàn)了Burg法。比較下對(duì)于語(yǔ)音信號(hào)通過(guò)Burg法和Yule-Walker法得到的譜,在較長(zhǎng)信號(hào)數(shù)據(jù)的情況下它們非常相似。load mtlbHburg = spectrum.burg(14); % 14th order modelpsd(Hburg,mtlb(1:512),'Fs',Fs,'NFFT',1024)Hyule
38、ar = spectrum.yulear(14); % 14th order modelpsd(Hyulear,mtlb(1:512),'Fs',Fs,'NFFT',1024)比較受噪聲干擾的信號(hào)的譜,分別使用Burg法和Welch法計(jì)算:randn('state',0)fs = 1000; % Sampling frequencyt = (0:fs)/fs; % One second worth of samplesA = 1 2; % Sinusoid amplitudesf = 150;140; % Sinusoid frequencies
39、xn = 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)需要注意的是,隨著B(niǎo)urg法模型階數(shù)的降低,由于正弦信號(hào)初始相位帶來(lái)的頻率偏移逐漸明顯。Covariance & Modified Covariance協(xié)方差和
40、修正協(xié)方差法AR譜估計(jì)的協(xié)方差算法基于最小化前向預(yù)測(cè)誤差而產(chǎn)生。而修正協(xié)方差算法是基于最小化前向和后向預(yù)測(cè)誤差而產(chǎn)生。工具箱中的spectrum.cov對(duì)象和pcov函數(shù),以及spectrum.mcov對(duì)象和pmcov函數(shù)實(shí)現(xiàn)了各自算法。比較語(yǔ)音信號(hào)分別用協(xié)方差算法和修正協(xié)方差算法得到的譜,它們幾乎完全一致,即使是對(duì)于短信號(hào)而言。load mtlbHcov = spectrum.cov(14); % 14th order modelpsd(Hcov,mtlb(1:64),'Fs',Fs,'NFFT',1024)Hmcov = spectrum.mcov(14);
41、 % 14th order modelpsd(Hmcov,mtlb(1:64),'Fs',Fs,'NFFT',1024)Subspace Methods子空間法MUSIC方法和特征向量分析法spectrum.music對(duì)象和pmusic函數(shù),以及spectrum.eigenvector對(duì)象和peig函數(shù)提供了兩種相關(guān)的譜分析方法:l spectrum.music對(duì)象和pmusic函數(shù)提供Schmidt提出的MUSIC算法。l spectrum.eigenvector對(duì)象和peig函數(shù)提供Johnson提出的EV算法。這兩種算法均基于對(duì)自相關(guān)矩陣的特征分析,用于對(duì)
42、頻率的估計(jì)。這種譜分析將數(shù)據(jù)相關(guān)矩陣的信息分為信號(hào)子空間或者噪聲子空間。特征分析方法簡(jiǎn)介考慮存在于白噪聲中的一些復(fù)正弦信號(hào)??梢园言撓到y(tǒng)的自相關(guān)矩陣R寫(xiě)作時(shí)信號(hào)自相關(guān)矩陣S與噪聲自相關(guān)矩陣W之和:信號(hào)自相關(guān)矩陣的特征向量與信號(hào)噪聲子空間有緊密聯(lián)系。S的特征向量v像信號(hào)向量一樣擴(kuò)展成為信號(hào)子空間。如果系統(tǒng)包含M個(gè)復(fù)正弦信號(hào),自相關(guān)矩陣階數(shù)為p,則特征向量至擴(kuò)展成自相關(guān)矩陣的噪聲子空間。頻率估計(jì)函數(shù)特征分析方法通過(guò)計(jì)算信號(hào)和噪聲子空間向量的某些函數(shù)來(lái)生成其頻率估計(jì)值。MUSIC和EV技術(shù)選擇一個(gè)函數(shù),它在一個(gè)輸入正弦信號(hào)頻率點(diǎn)上趨于無(wú)窮(分母趨于0)。使用數(shù)字技術(shù),得到的估計(jì)值在感興趣的頻點(diǎn)上具有
43、尖銳的峰值,這就意味著向量中可能沒(méi)用無(wú)窮大點(diǎn)。MUSIC估計(jì)由下面方程所示:此處N是特征向量的維數(shù),是復(fù)正弦信號(hào)向量:v表示輸入信號(hào)則相關(guān)矩陣的特征向量,vk是第k個(gè)特征向量;H代表共軛轉(zhuǎn)置。求和中的特征向量對(duì)應(yīng)了最小的特征值并張成噪聲空間(p是信號(hào)子空間維度)。表達(dá)式等價(jià)于一個(gè)傅里葉變換(向量由復(fù)指數(shù)組成)。這一形式對(duì)于數(shù)值計(jì)算有用,因?yàn)镕FT能夠?qū)γ恳粋€(gè)進(jìn)行計(jì)算,然后幅度平方再被求和。EV算法通過(guò)自相關(guān)矩陣的特征值對(duì)求和進(jìn)行加權(quán):工具箱中的pmusic函數(shù)和peig函數(shù)采用SVD(奇值分解)對(duì)信號(hào)計(jì)算,采用eig函數(shù)來(lái)分析自相關(guān)矩陣,并將特征向量歸為信號(hào)子空間和噪聲子空間。當(dāng)使用SVD的時(shí)
44、候,pmusic和peig并未顯式地計(jì)算相關(guān)矩陣,但是奇異值卻是特征值。先說(shuō)一下pwelch函數(shù)吧,相信你也查過(guò)了,我就直接把其他人的說(shuō)法粘貼過(guò)來(lái)了:Pxx,f=pwelch(x,window,noverlap,nfft,fs)Welch方法是一種修正周期圖功率譜密度估計(jì)方法,它通過(guò)選取的窗口對(duì)數(shù)據(jù)進(jìn)行加窗處理,先分段求功率譜之后再進(jìn)行平均。其中窗口的長(zhǎng)度N表示每次處理的分段數(shù)據(jù)長(zhǎng)度,Noverlap是指相鄰兩段數(shù)據(jù)之間的重疊部分長(zhǎng)度。N越大得到的功率譜分辨率越高(越準(zhǔn)確),但方差加大(及功率譜曲線(xiàn)不太平滑);N越小,結(jié)果的方差會(huì)變小,但功率譜分辨率較低(估計(jì)結(jié)果不太準(zhǔn)確)。1. 將信號(hào)分為多段,每段之間可以有overlapping,也可以沒(méi)有。2. 每一段
溫馨提示
- 1. 本站所有資源如無(wú)特殊說(shuō)明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請(qǐng)下載最新的WinRAR軟件解壓。
- 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請(qǐng)聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶(hù)所有。
- 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁(yè)內(nèi)容里面會(huì)有圖紙預(yù)覽,若沒(méi)有圖紙預(yù)覽就沒(méi)有圖紙。
- 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
- 5. 人人文庫(kù)網(wǎng)僅提供信息存儲(chǔ)空間,僅對(duì)用戶(hù)上傳內(nèi)容的表現(xiàn)方式做保護(hù)處理,對(duì)用戶(hù)上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對(duì)任何下載內(nèi)容負(fù)責(zé)。
- 6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容,請(qǐng)與我們聯(lián)系,我們立即糾正。
- 7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時(shí)也不承擔(dān)用戶(hù)因使用這些下載資源對(duì)自己和他人造成任何形式的傷害或損失。
最新文檔
- 山茶油訂購(gòu)合同協(xié)議書(shū)
- 道路共建協(xié)議書(shū)
- 邵昕離婚協(xié)議書(shū)
- 美伊退出協(xié)議書(shū)
- 酒店受傷協(xié)議書(shū)
- 酒店分成協(xié)議書(shū)
- 維護(hù)婚姻協(xié)議書(shū)
- 注冊(cè)類(lèi)人員引薦協(xié)議書(shū)
- 良田流轉(zhuǎn)協(xié)議書(shū)
- 礦石交易協(xié)議書(shū)
- 初中 初一 心理健康 生活中的小確幸 課件
- 輸液泵/微量注射泵使用技術(shù)操作考核評(píng)分標(biāo)準(zhǔn)
- 《微生物與免疫學(xué)》期末考試復(fù)習(xí)題及參考答案
- 梁若瑜著-十二宮六七二象書(shū)增注版
- 安全文明環(huán)保施工現(xiàn)場(chǎng)綜合規(guī)劃和詳細(xì)措施
- 《第二單元 遼宋夏金元時(shí)期:民族關(guān)系發(fā)展和社會(huì)變化》單元梳理
- 外研版三年級(jí)英語(yǔ)下冊(cè)全冊(cè)教材分析解讀
- 巴蜀文化(課堂PPT)課件
- 質(zhì)量部組織架構(gòu)
- 電氣裝置安裝工程接地裝置施工及驗(yàn)收規(guī)范——50169-2006
- 水電站自動(dòng)化運(yùn)行專(zhuān)業(yè)術(shù)語(yǔ)
評(píng)論
0/150
提交評(píng)論