對功率譜估計常用方法的探討及應(yīng)用分析

DSP課程設(shè)計對功率譜估計常用方法的探討及應(yīng)用分析進行傅里葉變換在頻域中研究信號，是研究確定性信號最簡單且有效的手段，但在現(xiàn)代信號分析中，對于常見的隨機信號，不可能用清楚的數(shù)學(xué)關(guān)系式來描述，其傅里葉變換更不存在，轉(zhuǎn)而可以利用給定的N個樣本數(shù)據(jù)估計一個平穩(wěn)隨機信號的功率譜密度。功率譜估計是數(shù)字信號處理的重要研究內(nèi)容之一。功率譜估計可以分為經(jīng)典功率譜估計和現(xiàn)代功率譜估計。本文介紹了各種經(jīng)典功率譜估計方法，不僅從理論上對各種方法的譜估計質(zhì)量進行了分析比較，而且通過Matlab進行了仿真。在對經(jīng)典譜估計進行討論之后，還分析了現(xiàn)代譜估計即參數(shù)譜估計方法，通過觀測數(shù)據(jù)估計參數(shù)模型再按照求參數(shù)模型輸出功率的方法估計信號功率譜。現(xiàn)代譜估計的內(nèi)容極其豐富，設(shè)計的學(xué)科及應(yīng)用的領(lǐng)域都相當(dāng)廣泛，至今每年都有大量的科研成果出來。在本文的最后利用現(xiàn)代譜估計的方法討論了功率譜方法在噪聲源信號識別中的應(yīng)用。文章還給出了常見譜估計方法的比較，便于深刻理解各種方法的特點，從而在實際工作中做出合理的選擇。1.功率譜方法的發(fā)展功率譜估計是隨機信號處理的重要內(nèi)容，其技術(shù)淵源很長，而且在過去的40余年中獲得了飛速的發(fā)展。涉及到信號與系統(tǒng)、隨機信號分析、概率統(tǒng)計、矩陣代數(shù)等一系列的基礎(chǔ)學(xué)科，廣泛應(yīng)用于人民的日常生活及軍事、工業(yè)、農(nóng)業(yè)活動中，是一個具有強大生命力的研究領(lǐng)域。本文將簡要回顧一下功率譜估計的發(fā)展歷程，對常用的一些方法進行總結(jié)。功率譜的估計方法有很多，主要有經(jīng)典譜估計和現(xiàn)代譜估計。經(jīng)典譜估計又可以分成兩種：一種是BT法，也叫間接法；另一種是直接法又稱周期圖法?，F(xiàn)代譜估計的方法又大致可分為參數(shù)模型譜估計和非參數(shù)模型譜估計，前者有AR模型、MA模型、ARMA模型、PRONY模型等，后者有最小方差方法、多分量的MUSIC方法等。1.1功率譜研究的發(fā)展過程功率譜估計是數(shù)字信號處理的主要內(nèi)容之一，主要研究信號在頻域中的各種特征，目的是根據(jù)有限數(shù)據(jù)在頻域內(nèi)提取被淹沒在噪聲中的有用信號。下面對譜估計的發(fā)展過程做簡要回顧：英國科學(xué)家牛頓最早給出了“譜”的概念。后來，1822年，法國工程師傅立葉提出了著名的傅立葉諧波分析理論。該理論至今依然是進行信號分析和信號處理的理論基礎(chǔ)。傅立葉級數(shù)提出后，首先在人們觀測自然界中的周期現(xiàn)象時得到應(yīng)用。19世紀(jì)末，Schuster提出用傅立葉級數(shù)的幅度平方作為函數(shù)中功率的度量，并將其命名為“周期圖”（periodogram）。這是經(jīng)典譜估計的最早提法，這種提法至今仍然被沿用，只不過現(xiàn)在是用快速傅立葉變換（FFT）來計算離散傅立葉變換（DFT），用DFT的幅度平方作為信號中功率的度量。1958年，R,Blackman和J.Tukey首先提出BT法，并命名為布萊克曼-杜基譜估計器（簡稱BT譜估計器）。這種方法是先按照有限個觀測數(shù)據(jù)估計自相關(guān)函數(shù)，再對其求傅里葉變換得到功率譜。在1965年FFT未出現(xiàn)以前，BT法一直是最常用的方法。周期圖較差的方差性能促使人們研究另外的分析方法。1927年，Yule提出用線性回歸方程來模擬一個時間序列。Yule的工作實際上成了現(xiàn)代譜估計中最重要的方法——參數(shù)模型法譜估計的基礎(chǔ)。Walker利用Yule的分析方法研究了衰減正弦時間序列，得出Yule-Walker方程，可以說，Yule和Walker都是開拓自回歸模型的先鋒。1930年，著名控制理論專家Wiener在他的著作中首次精確定義了一個隨機過程的自相關(guān)函數(shù)及功率譜密度，并把譜分析建立在隨機過程統(tǒng)計特征的基礎(chǔ)上，艮L“功率譜密度是隨機過程二階統(tǒng)計量自相關(guān)函數(shù)的傅立葉變換”，這就是Wiener—Khintchine定理。該定理把功率譜密度定義為頻率的連續(xù)函數(shù)，而不再像以前定義為離散的諧波頻率的函數(shù)。1949年，Tukey根據(jù)Wiener—Khintchine定理提出了對有限長數(shù)據(jù)進行譜估計的自相關(guān)法，即利用有限長數(shù)據(jù)估計自相關(guān)函數(shù)，再對該自相關(guān)函數(shù)球傅立葉變換，從而得到譜的估計。1958年，Blackman和Tukey在出版的有關(guān)經(jīng)典譜估計的專著中討論了自相關(guān)譜估計法，所以自相關(guān)法又叫BT法。周期圖法和自相關(guān)法都可用快速傅立葉變換算法來實現(xiàn)，且物理概念明確，因而仍是目前較常用的譜估計方法。1948年，Bartlett首次提出了用自回歸模型系數(shù)計算功率譜。自回歸模型和線性預(yù)測都用到了1911年提出的Toeplitz矩陣結(jié)構(gòu)，Levinson曾根據(jù)該矩陣的特點于1947年提出了解Yule-Walker的快速計算方法。這些工作為現(xiàn)代譜估計的發(fā)展打下了良好的理論基礎(chǔ)。1965年，Cooley和Tukey提出的FFT算法，也促進了譜估計的迅速發(fā)展。現(xiàn)代譜估計的提出主要是針對經(jīng)典譜估計（周期圖和自相關(guān)法）的分辨率和方差性能不好的問題。1967年，Burg提出的最大嫡譜估計，即是朝著高分辨率譜估計所作的最有意義的努力。雖然，Bartlett在1948年，Parzem于1957年都曾經(jīng)建議用自回歸模型做譜估計，但在Burg的論文發(fā)表之前，都沒有引起注意。現(xiàn)代譜估計的內(nèi)容極其豐富，涉及的學(xué)科及應(yīng)用領(lǐng)域也相當(dāng)廣泛，至今，每年都有大量的論文出現(xiàn)。目前尚難對現(xiàn)代譜估計的方法作出準(zhǔn)確的分類。從現(xiàn)代譜估計的方法上，大致可以分為參數(shù)模型譜估計和非參數(shù)模型譜估計，前者有AR模型、MA模型、ARMA模型、PRONY模型等；后者有最小方差方法、多分量的MUSIC方法等。非參數(shù)模型譜估計的特點是其模型不是用有限參數(shù)來描述，而直接由相關(guān)函數(shù)序列得到，這種方法能提高低信噪比時的譜分辨率。參數(shù)模型譜估計是先根據(jù)過程的先驗信息或者一些假定，建立一個數(shù)學(xué)模型來表示所給定采樣數(shù)據(jù)的過程，或者選擇一個較好的近似實際模型，而后利用采樣數(shù)據(jù)序列或者自相關(guān)序列，估計該模型的參數(shù)，最后把參數(shù)代入到該模型對應(yīng)的理論功率譜表達式，得到所需要的譜估計。目前大量的論文集中在模型參數(shù)的求解上，以求得到速度更快、更穩(wěn)健、統(tǒng)計性能更好的算法。1.2功率譜估計方法提出在通信系統(tǒng)中，往往需要研究具有目中統(tǒng)計特性的隨機信號。由于隨機信號是一類持續(xù)時間無限長，具有無限大能量的功率信號，它不滿足傅里葉變換條件，而且也不存在解析表達式，因此就不能夠應(yīng)用確定信號的頻譜計算方法去分析隨機信號的頻譜。然而，雖然隨機信號的頻譜不存在，但其相關(guān)函數(shù)是可以確定的。如果隨機信號是平穩(wěn)的，那么其相關(guān)函數(shù)的傅里葉變換就是它的功率譜密度函數(shù)，簡稱功率譜。功率譜反映了單位頻帶內(nèi)隨機信號的一個樣本信號來對該隨機過程的功率譜密度函數(shù)做出估計。1.3功率譜估計應(yīng)用及用途功率譜估計有著極其廣泛的應(yīng)用，不僅在認(rèn)識一個隨機信號時，需要估計它的功率譜。它還被廣泛地應(yīng)用于各種信號處理中。在信號處理的許多場所，要求預(yù)先知道信號的功率譜密度(或自相關(guān)函數(shù))。例如，在最佳線性過濾問題中，要設(shè)計一個維納濾波器就首先要求知道(或估計出)信號與噪聲的功率譜密度(或自相關(guān)函數(shù))。根據(jù)信號與噪聲的功率譜(或4(m))才能設(shè)計出能夠盡量不失真的重現(xiàn)信號，而把XX噪聲最大限度抑制的維納濾波器。常常利用功率譜估計來得到線性系統(tǒng)的參數(shù)估計。例如，當(dāng)我們要了解某一系統(tǒng)的幅頻特性H(o)時，可用一白色噪聲sS通過該系統(tǒng)。再從該系統(tǒng)的輸出樣本y(n)估計功率譜密度P(o)。由于白色噪聲的PSD(用pjo)表示)為一常數(shù)即p①(o)=b\，于是有：P(o)=o*[H(o)2(1-1)故通過估計輸出信號的PSD，可以估計出系統(tǒng)的頻率特性H(o)|(模特性)。從寬帶噪聲中檢測窄帶信號。這是功率譜估計在信號處理中的一個重要用途。但是這要求功率譜估計有足夠好的頻率的分辨率，否則就不一定能夠清楚地檢測出來。所謂譜估計的分辨率可以粗略地定義為能夠分辨出的二個分立的譜分量間的最小頻率間隙(距)。提高譜估計的分辨率已成為目前譜估計研究中的一個重要方向功率譜估計就是通過信號的相關(guān)性估計出接受到信號的功率隨頻率的變化關(guān)系，實際用途有濾波，信號識別(分析出信號的頻率)，信號分離，系統(tǒng)辨識等。譜估計技術(shù)是現(xiàn)代信號處理的一個重要部分，還包括空間譜估計，高階譜估計等。維納濾波、卡爾曼濾波，可用于自適應(yīng)濾波，信號波形預(yù)測等(火控系統(tǒng)中的飛機航跡預(yù)判)。2經(jīng)典譜估計2.1周期圖法周期圖法又稱直接法。它是從隨機信號x(n)中截取N長的一段，把它視為能量有限x(n)真實功率譜七(球)的估計七(幻w)的抽樣.周期圖這一概念早在1899年就提出了，但由于點數(shù)N一般比較大，該方法的計算量過大而在當(dāng)時無法使用。只是1965年FFT出現(xiàn)后，此法才變成譜估計的一個常用方法。周期圖法包含了下列二條假設(shè)：認(rèn)為隨機序列是廣義平穩(wěn)且各態(tài)遍歷的，可以用其一個樣本x(n)中的一段,(n)來估計該隨機序列的功率譜。這當(dāng)然必然帶來誤差。由于對x(n)采用DFT，就默認(rèn)x(n)在時域是周期的，以及x(k)在頻域NNN是周期的。這種方法把隨機序列樣本x(n)看成是截得一段Xn(n)的周期延拓，這也就是周期圖法這個名字的來歷。與相關(guān)法相比，相關(guān)法在求相相關(guān)函數(shù)%(m)時將Xn(n)以外是數(shù)據(jù)全都看成零，因此相關(guān)法認(rèn)為除Xn(n)外x(n)是全零序列，這種處理方法顯然與周期圖法不一樣。N但是，當(dāng)相關(guān)法被引入基于FFT的快速相關(guān)后，相關(guān)法和周期圖法開始融合。比我們發(fā)現(xiàn)：如果相關(guān)法中M=N，不加延遲窗，那么就和充(N-1)個零的周期圖法一樣了。簡單地可以這樣說：周期圖法是M=N時相關(guān)法的特例。因此相關(guān)法和周期圖法可結(jié)合使用。2.2相關(guān)法譜估計(BT)這種方法以相關(guān)函數(shù)為媒介來計算功率譜，所以又叫間接法。它是1958年由Blackman和Tukey提出。這種方法的具體步驟是：第一步：從無限長隨機序列x(n)中截取長度N的有限長序列列Xn(n)第二步：由N長序列Xn(n)求(2M-1)點的自相關(guān)函數(shù)Rx(m)序列。即Rx(m)=—^1x(n)x(n+m)(2-1)n=0這里，m=-(M-1)...,-1,0,1...,M-1，M^N，Rx(m)是雙邊序列，但是由自相關(guān)函數(shù)的偶對稱性式，只要求出m=0,。。。,M-1的傅里葉變換，另一半也就知道了。第三步：由相關(guān)函數(shù)的傅式變換求功率譜。即S(ejw)=覽R(m)e-jwm(2-2)尤Xm=-(M-1)以上過程中經(jīng)歷了兩次截斷，一次是將x(n)截成N長，稱為加數(shù)據(jù)窗，一次是將想x(n)截成(2M-1)長，稱為加延遲窗。因此所得的功率譜僅是近似值，也叫譜估計，式中的七(ejw)代表估值。一般取M<<N，因為只有當(dāng)M較小時，序列傅式變換的點數(shù)才破小，功率譜的計算量才不至于大到難以實現(xiàn)，而且譜估計質(zhì)量也較好。因此，在FFT問世之前，相關(guān)法是最常用的譜估計方法。當(dāng)FFT問世后，情況有所變化。因為截斷后的七(n)可視作能量信號，由相關(guān)卷積定理可得(2-3)Rx(m)=N[x(m)*x(-m)](2-3)這就將相關(guān)化為線性卷積，而線性卷積又可以用快速卷積來實現(xiàn)。我們可對上式兩邊取(2N-1)點DFT，則有(2-4)R,(m)=1[x(k)*x(K)=1\X(K)|2(2-4)xN2N-12N-1N2N-1于是將時域卷積變?yōu)轭l域乘積，用快速相關(guān)求Rx(m)的完整方案如下:1.對N長xN(n)的充(N-1)個零，成為(2N-1)長的。2求(2N-1)點的FFT，得X(K)=籠-2x(n)W-mk。2N-12N-12N-1N=0求￡|X霜1(K)2。由DFT性質(zhì)，乂咨1(n)是純實的，％\k)滿足共軛偶對稱，而N|x1(K)|2一定是實偶的，且以(2N-1)為周期。求(2N-1)點的IFFT：Rx(m)=2；-1￥NX2N-1(幻2Wm(2-5)k=-(N-1)這里-1!X2n1(K)2是實偶的，m=-(N-1)...0...N-1。本來IFFT求和范圍是0至2N-2,由于TX2N-1(K)2的實偶性與周期性，求和范圍改為-(N-1)至(N-1)不影響計算結(jié)果。同理可將m的范圍改為-(N-1)￥(N-1)。上述的快速相關(guān)中，充零的目的是為了能用圓周卷積代替線性卷積，以便進一步采用快速卷積算法?？焖傧嚓P(guān)輸出是-(N-1)至(N-1)的2N-1點，加W(m)窗M后截取的是-(M-1)至(M-1)的頻段，最后作(2M-1)點FFT，得Sx(k)。我們注意到：如果數(shù)據(jù)點數(shù)與自相關(guān)序列點數(shù)相同即M=N，^U(2N-1)點的IFFT后緊跟一個(2N-1)點的FFT，利用Rx(m)的對稱性，F(xiàn)T運算框的計算式變?yōu)镾(K)=支'Rx(m)W-mkX(2-6)m=-(N-1)由于N=M并假設(shè)窗形狀是矩形的，第二次Wm(m)的截斷就不需要了。比較式(3-5)和式(3-6)，S(k)=FFT[R《m)],Rx(m)=IFFT[N|x2州K)|2]正反傅氏變換可以抵消，直接得(2-7)S*)=NX2n-1(K)2(2-7)為了實行基2FFT，也可將(2N-1)點換成2N點，這樣做不影響結(jié)果的正確性。2.3巴特利特(Bartlett)平均周期圖的方法首先讓我們來看一下為什么周期圖經(jīng)過某種平均(或平滑)后會使它的方差當(dāng)N-8時趨于零，達到一致估計的目的。如果氣，x2,…，七是不相關(guān)的隨機變量，每一個具有期望值p，方差a2，則可以證明它們的數(shù)學(xué)平均如=(x1+x2+...+x)lL的期望值等于p，數(shù)學(xué)平均的方差等于a2/L，即：12'E[x]=Le\x+xHFx]=L.加=pVar\x]=e[X—E(x))2LE[x2]-(e\x)=—E[x+x+Fx)^—P2VarL212L>—p2E^x2+x2+?.?+x2)+左LEj=1i=1母j工ZeLxLEeIM]>—p2IjjIj=1ij=1i豐j=卬(L—1川=L曰2-卬2所以Varix]=—&C21+L2口2-m2**1所以Varix]=—&C21+L2口2-m2**1-日2=——&!11x2■—Lp2L2i-i=1L2Li=1iL21(切氣］)2+J2I(E[x])2}一L1b2=—Lb2=—L2L(2-8)由(2-8)可見，L個平均的方差比每個隨機變量的單獨方差小L倍。當(dāng)LT3Varix〕T0，可達到一致譜估計的目的。因而降低估計量的方差的一種有效方法是將若十個獨立估計值進行平均。把這種方法應(yīng)用于譜估計通常歸功于Bartlett。Bartlett平均周期圖的方法是將序列x(n)(0<n<N-1)分段求周期圖再平均。設(shè)將x(n)分成L段，每段有M個樣本，因而N=LM，第i段樣本序列可寫成xi(n)=x(n+iM—M)0<n<M-1,1<i<L第i段的周期圖為M1xj(n)e-j^nn=0如果m>M,4xx(m)很小，則可假定各段的周期圖Im(①)是互相獨立的。對功率譜密度的概念的討論，譜估計可定義為L段周期圖的平均，即p(①)=1￡ii(①)i=1(2-9)于是它的期望值為E^P(①)11￡e^Ii(①)LE^Ii(①J_xxLMM(o)1二j”P(6)W(ej(?-6))d62兀一”xxBEPxx1-八j”P(6)2兀M一兀xxsinK?-e)M..2]

sin(L-6)/2(2-10)這里M=N/L，因此Bartlett估計的期望值是真實譜PJo)與三角窗函數(shù)的卷積。由于三角窗函數(shù)不等于8函數(shù)，所以Bartlett估計也是有偏估計即Bias。0，但

當(dāng)N—s時，Bias—0。由于我們假定各段周期圖是相互獨立的，所以可按式(3-8)得到下式:Vark(Vark(rn)11Var\lxxLM(①)](2-11)由此可見，隨著L的增加Vai^P(w)k下降的，當(dāng)L—s時，Var^P(w)］—0。因xxxx此Bartlett估計是一致估計。比較式(2-10)的eP(w)］與式(2-1)的血(°)］可見在二種情況的估計量的期望值xxN都是真值P(W)與窗口函數(shù)Wg)的卷積形式，但后者將前者WB中N改為M，今式(2-10)中w(eg)的主瓣寬度大于后式中的主瓣寬度，因而Bias。(W)］，而主瓣愈寬顯然使分辨率愈差。因此Bias可用來說明譜的N偏倚△xxBiasM=N/L<N。因而使wb(j)主瓣的寬度增大。由于主瓣的寬度愈窄愈接近5函數(shù)，分辨率，Bias愈大說明譜分辨率愈差。一個固定的記錄長度N，周期圖分段的數(shù)目L愈大將使方差愈小，但M也愈小，因而使Bias今式(2-10)中w(eg)的主瓣寬度大于后式中的主瓣寬度，因而Bias。(W)］，而主瓣愈寬顯然使分辨率愈差。因此Bias可用來說明譜的N偏倚△xxBias由此可見Bartlett法使譜估計的方差減小是用增加Bias以及降低譜分辨率的代價換來的。實際上，當(dāng)N一定時，Bias與Var的互換性是譜估計的一個固有特性。為了說明經(jīng)平均后的周期圖作為功率譜估計的實際效果，設(shè)有一零均值高斯分布的隨機過程，其功率譜密度為■?”)■?”)(1-az-1)(1-az)a=1—*'0.1這一功率譜密度是由一零值高斯分布單位方差的噪聲序列通過一個其H(z)=1/(1-az-1)的濾波器形成的。為了簡單設(shè)選用N=8的矩形窗函數(shù)。其周期圖的期望值(用E［18(w)］表示)與真值匕(w)均表示在圖2.4中。它說明N=8的周期圖可以得到的Bias的情況。圖2-2表示xM=8分4段與16段二種情況平均后的周期圖。顯然L=4的方差比L=16的大。將L=16的曲線與圖2.4的曲線比較可見在這種估計中誤差的大部分起因于Bias而不是方差(因為二種情況均為M=8，Bias相仿，誤差也相仿)。M=16,L=2及8的周期圖表示在圖2-3。圖2-3中L不同造成的影響也是明

顯的。但是這二個曲線的起伏都很大，因此有理由認(rèn)為誤差主要起因于方差。比較M=8與M=16的周期圖可見，M愈大（即§2.3中的N愈大），將使周期圖的起伏愈增快的結(jié)論，在這里也同樣成立。比較圖2-2與圖2-3發(fā)現(xiàn)在這個例子中最好的選擇是應(yīng)用L=16,M=8的估計而不是L=8、M=16的估計，即寧可減小方差，犧牲Bias。在實際中，當(dāng)然功率譜密度的真值是不知道的。但是譜的窗口函數(shù)以及關(guān)于功率譜密度的某些信息往往是預(yù)函先知道的。通過改變M和L以及利用預(yù)先知道的情況，通常可以很好地進行選擇。平均周期圖的方法特別適合于應(yīng)用FFT算法。因此在FFT出現(xiàn)以后這個方法比下面將要討論的利用窗口函數(shù)的處理法用得更多。而在FFT出現(xiàn)以前主要是用窗口函數(shù)處理法平滑周期圖。圖2-1P(o圖2-1P(o)與E[18(0)]的特性圖2-2平滑后的周期圖（每段取8個數(shù)據(jù)）圖2-3平均后的周期圖（每段取16個數(shù)據(jù)）2.4Welch法Welch提出了對Bartlett法的修正使更適合于用FFT進行計算。他主要提出二方面的修正，其一是選擇適當(dāng)?shù)拇昂瘮?shù)o⑴，并在周期圖計算前直接加進法，這樣得到的每一段的周期圖為電1xi(n)o(n)e-on=0這里U=_L*1W2(n)為歸一化因子，而Bartlett法每段的周期圖為n=0/、1N-1,、2I⑺(電1xi(n)o(n)e-on=0這樣加窗函數(shù)的優(yōu)點是無論什么樣的窗函數(shù)均可使譜估計非負。其二是在分段時，可使各段之間有重迭，這樣將會使方差減小(當(dāng)N與M一定時)。他建議可以重迭50%。3經(jīng)典譜估計方法編程與分析3.1直接法直接法又稱周期圖法，它是把隨機序列x(n)的N個觀測數(shù)據(jù)視為一能量有限的序列，直接計算x(n)的離散傅立葉變換，得X(k),然后再取其幅值的平方，并除以N，作為序列x(n)真實功率譜的估計。Matlab代碼示例：clear;Fs=1000;%采樣頻率n=0:1/Fs:1;%產(chǎn)生含有噪聲的序列改變n的取值范圍觀察圖形的變換xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));window=boxcar(length(xn));%矩形窗nfft=1024;[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法plot(f,10*log10(Pxx));

圖3-1周期圖法中N=1000功率譜-1UIIIIIIII050100160200260300360400460600圖3-2周期圖中N=100000功率譜圖3-3周期圖當(dāng)N=100功率譜隨著采樣點數(shù)的增加，該股集是漸進無騙的。從圖中可以看出，采用周期突發(fā)估計得出的功率譜很不平滑，相應(yīng)的估計協(xié)方差比較大。而且采用增加采樣點的辦法也不能吃周期圖變得更加平滑，這是周期圖法的缺點。周期圖法得出的估計譜方差特性不好：當(dāng)數(shù)據(jù)長度N太大時，撲線的起伏加劇；N太小時譜的分辨率又不好。對其改進的主要方法有二種，即平均和平滑，平均就是將截取的數(shù)據(jù)段七(n)再分成L個小段，分別計算功率譜后取功率譜的平均，這種方法使估計的方差減少，但偏差加大，分辨率下降。平滑是用一個適當(dāng)?shù)拇昂瘮?shù)W(em)與算出的功率譜5{ejw)進行卷積，使譜線平滑。這種方法得出的譜估計是無偏的，方差也小，但X分辨率下降。3.2間接法間接法先由序列x(n)估計出自相關(guān)函數(shù)R(n)，然后對R(n)進行傅立葉變換，便得到x(n)的功率譜估計。Matlab代碼示例：clear;Fs=1000;%采樣頻率n=0:1/Fs:1;%產(chǎn)生含有噪聲的序列xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));nfft=1024;cxn=xcorr(xn,'unbiased');%計算序列的自相關(guān)函數(shù)CXk=fft(cxn,nfft);Pxx=abs(CXk);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));plot(k,plot_Pxx);圖3-4直接發(fā)功率譜圖當(dāng)M=N-1時,BT法與周期圖法估計出的功率譜是一樣的;當(dāng)M<N-1時,BT法的偏差大于周期圖法，在窗函數(shù)滿足一定條件時是漸進無偏估計;方差小于周期圖的方差;分辨率比周期圖法低,與窗函數(shù)的選擇有關(guān)。BT法的缺點在于當(dāng)M一N時，R(m)的方差很大,使譜估計質(zhì)量下降;由R(m)得到的$(明不一定為正值，從而可能失去功率譜的物理意義。3.3改進的直接法：對于直接法的功率譜估計，當(dāng)數(shù)據(jù)長度N太大時，譜曲線起伏加劇，若N太小，譜的分辨率又不好，因此需要改進。3.3.1Bartlett法對周期圖法改進的思想是將信號分段進行估計，然后再將這些估計結(jié)果進行平均，從而減小估計的協(xié)方差，是功率譜圖變得比直接法更平滑。增加分段數(shù)可以進一步減低估計的協(xié)方差，然而若每段中的數(shù)據(jù)點數(shù)太少，就會使估計的頻率分辨率下降很多。從樣本信號序列總點數(shù)一定的條件下，可以采用使分段相互重疊的方法來增加分段數(shù)，從而保持每段信號點數(shù)不變，這樣就在保證頻率分辨率的前提下進一步降低估計協(xié)方差。而一般的Bartlett法是通過降低分辨率來降低其方差的。Bartlett平均周期圖的方法是將N點的有限長序列x(n)分段求周期圖再平均。Matlab代碼示例：clear；Fs=1000;n=0:1/Fs:1;xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));nfft=1024;window=boxcar(length(n));%矩形窗noverlap=0;%數(shù)據(jù)無重疊p=0.9;%置信概率[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));plot_Pxxc=10*log10(Pxxc(index+1));figure(1)plot(k,plot_Pxx);pause;figure(2)plot(k,[plot_Pxxplot_Pxx-plot_Pxxcplot_Pxx+plot_Pxxc]);圖3-5bartlett法估計功率譜圖現(xiàn)在比較常用的改進方法是Welch法，又叫加權(quán)交疊平均法，簡記為WOSA法，這種方法以加窗（加權(quán)）求取平滑，以分段重疊求得平均，因此集平均與平滑的優(yōu)點于一體，同時也不可避免帶有兩者的缺點，因此歸根到底是一種折中。其主要步驟是：（1）將N長的數(shù)據(jù)段分成L個小段，每小段M點，相鄰小段間交疊M/2點（即2:1分段）。因為L（M/2）+M/2=N，所以段數(shù)(3-1),N-M/2L=M/2(3-1)（2）對各小段加同樣的平滑窗由3）后起傅氏變換（3-2）X（ejw）=勿x（n）ax（n）e-jwn,i=1,...L,n=0（3-2）（3）用下式求各小段功率譜的平均（3-3）Si（ejw）=1工上|X.（ejw）|2i=1■IAJ-1U=——▽oj2（?）這里，財口代表窗函數(shù)平均功率，MU是M長窗函數(shù)的能量。（3-3）Matlab代碼示例：clear;Fs=1000;n=0:1/Fs:1;xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));nfft=1024;window=boxcar(100);%矩形窗window1=hamming(100);%海明窗window2=blackman(100);%blackman窗noverlap=20;range='half;[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);plot_Pxx=10*log10(Pxx);plot_Pxx1=10*log10(Pxx1);plot_Pxx2=10*log10(Pxx2);figure(1)plot(f,plot_Pxx);pause;figure(2)plot(f,plot_Pxx1);pause;figure(3)0圖3-8blackman窗35IIIIIIIIIWelch法對Bartlett^進亍了兩方面的修正，一是選擇適當(dāng)?shù)拇昂瘮?shù)w(n),并再周期圖計算前直接加進去，加窗的優(yōu)點是無論什么樣的窗函數(shù)均可使譜估計非負。二是在分段時，可使各段之間有重疊，這樣會使方差減小。不同窗函數(shù)的Welch譜估計在選擇窗函數(shù)時，一般有如下要求：窗口寬度M要遠小于樣本序列長度N，以排除不可靠的自相關(guān)值；當(dāng)平穩(wěn)信號為實過程時，為保證平滑周期圖和真是功率譜也是實偶函數(shù)，平滑窗函數(shù)必須是實偶對稱的；平滑窗函數(shù)應(yīng)當(dāng)在m=0出游峰值，并且m隨絕對值增加而單調(diào)下降，使可靠的自相關(guān)值有較大的權(quán)值；功率譜是頻率的非負函數(shù)且周期圖是非負的，因而要求窗函數(shù)的Fourier變換是非負的。仍以上述平穩(wěn)隨機信號為例，采樣頻率、采樣點數(shù)、FFT點數(shù)、窗長度及重疊數(shù)據(jù)不變，窗函數(shù)采用矩形窗、Blackman窗、Hamming窗，仿真結(jié)果如圖所示。由仿真結(jié)果可以看出：使用不同的窗函數(shù)譜估計的質(zhì)量是不一樣的，矩形窗的主瓣較窄，分辨率較好，但方差較大，噪聲水平較高；而Blackman窗和Hamming窗的主瓣較寬，分辨率較低，但方差較小，噪聲水平較低。因此，在進行譜分析時選擇何種窗函數(shù)，要視具體情況而定。如果強調(diào)高分辨率，能精確讀出主瓣頻率，而不關(guān)心幅度的精度，例如測量震動物體的自震頻率，可以選用主瓣寬度比較窄的矩形窗；對受到強干擾的窄帶信號若干擾靠近信號，則可選用旁瓣幅度較小的窗函數(shù)，若離開通帶較遠，則可選用漸近線衰減速度比較快的窗函數(shù)?？傊?，要針對不同的信號和不同的處理目的來選擇合適的窗函數(shù)，這樣才能得到良好的效果。3.4經(jīng)典譜估計方法的比較總結(jié)周期圖法和BT法的特點是離散性大，曲線粗糙，方差較大，但是分辨率較高；Bartlett法和Welch法的收斂性較好，曲線平滑，方差較小，但是功率譜主瓣較寬，分辨率低，這是由于對隨機序列加窗截斷所引起的Gibbs效應(yīng)造成的；(3)與Bartlett法相比，Welch法的估計曲線比較粗糙，但是分辨率較好，原因是Welch法中對數(shù)據(jù)進行截斷時加的是Hanning窗，而在Bartlett法中使用的是矩形窗，相對于矩形窗，Hanning窗的主瓣包含更多的能量，因而使功率譜的主瓣較窄，分辨率較高。由上述理論分析及仿真實驗可知，Welch法采用加窗交疊求功率譜，可以有效減小方差和偏差，一般情況下能接近一致估計的要求，因而得到廣泛應(yīng)用。同時還可以發(fā)現(xiàn)，對信號加不同的窗函數(shù)，譜估計的質(zhì)量是不同的。在經(jīng)典譜估計中，無論是周期圖法還是其改進的方法，都存在著頻率分辨率低、方差性能不好的問題，原因是譜估計時需要對數(shù)據(jù)加窗截斷，用有限個數(shù)據(jù)或其自相關(guān)函數(shù)來估計無限個數(shù)據(jù)的功率譜，這其實是假定了窗以外的數(shù)據(jù)或自相關(guān)函數(shù)全為零，這種假定是不符合實際的，正是由于這些不符合實際的假設(shè)造成了經(jīng)典譜估計分辨率較差。另外，經(jīng)典譜估計的功率譜定義中既無求均值運算又無求極限運算，因而使得譜估計的方差性能較差，當(dāng)數(shù)據(jù)很短時，這個問題更為突出。如何選取最佳窗函數(shù)、提高頻譜分辨率，如何在短數(shù)據(jù)情況下提高信號譜估計質(zhì)量，還需要進一步研究。4.現(xiàn)代譜估計現(xiàn)代功率譜估計即參數(shù)譜估計方法是通過觀測數(shù)據(jù)估計參數(shù)模型再按照求參數(shù)模型輸出功率的方法估計信號功率譜，主要是針對經(jīng)典譜估計的分辨率低和方差性能不好等問題提出的。常用模型有ARMA模型、AR模型、MA模型,AR模型的正則方程是一組線性方程，而MA和ARMA模型是非線性方程。由于AR模型具有一系列良好的性能，因此被研究最多也得到最廣泛的應(yīng)用。本節(jié)將較為詳細的討論AR模型，并對MA和ARMA模型譜估計方法做簡要的討論。4.1AR模型AR模型又稱為自回歸模型，是一個全極點的模型,可用如下差分方程來表示：x(n)=—￡ax(n-k)+w(n)k=1(4-1)其中，p是AR模型的階數(shù)，『k}=l，2，…，p是p階AR模型的參數(shù).將該模型記為AR(p)，它的系統(tǒng)轉(zhuǎn)移函數(shù)為

旦=1序一1+眼孕(4-2)kk=1在功率譜估計中，若觀測的數(shù)據(jù)x(n)是平穩(wěn)隨機過程，則該系統(tǒng)的輸入w(n)也可認(rèn)為是平穩(wěn)的，因而根據(jù)線性系統(tǒng)對平穩(wěn)隨機信號的響應(yīng)理論可得觀測數(shù)據(jù)的功率譜為(4-3)P(w)=b2ih(ejw)|2=ZWxw11+2^ae~jkwk(4-2)(4-3)k=1由式可知，利用ar模型進行功率譜估計的實質(zhì)是求解模型系數(shù)｛a｝和b2的問題.將式(1)兩端乘以x(n-m)求平均(數(shù)學(xué)期望)，可以求得觀測數(shù)據(jù)的AR(p)模型參數(shù)與自相關(guān)函數(shù)的關(guān)系式為R(m)R(m)={XX-^LaR(m—k)k=1XX—￡aR(m—k)+bk=1kRX(—m)m>02m=0wm<0(4-4)可見，階AR模型輸出的相關(guān)函數(shù)具有遞推的性質(zhì)，因而選用AR模型進行譜估計只需較少的觀測數(shù)據(jù)將式(4)寫成矩陣形式得-R-R(0)R(1)...R(P)Tb2R(1):R(0)R(p—1):a1:=w0LR(p)R(p—1)R(0)a1-p」L0」(4-5)進而按式(3)求得信號功率譜的估上式就是著名的Yule-Walker(Y—W)方程.它表明，只要已知觀測數(shù)據(jù)的自相關(guān)數(shù)，就能求出AR模型參數(shù)｛ak｝進而按式(3)求得信號功率譜的估該模型的現(xiàn)在輸出值是它本身過去

它們之間有著非常密切的關(guān)系，即它另外，從AR模型的差分方程式可知值的回歸，這與預(yù)測器存在著一定的相似性，們的系統(tǒng)函數(shù)互為倒數(shù)，也就是說預(yù)測誤差濾波器A(z)就是AR信號模型H(z)p的逆濾波器.因此通過預(yù)測誤差濾波器優(yōu)化設(shè)計使預(yù)測均方誤差最小就可求得AR信號模型的最優(yōu)參數(shù)，即P階線性預(yù)測器的預(yù)測系數(shù)｛a(k)｝等于p階AR模型的系數(shù)｛ak｝，其最小均方預(yù)測誤差E等于自噪聲方差°：。該模型的現(xiàn)在輸出值是它本身過去

它們之間有著非常密切的關(guān)系，即它因此，根據(jù)上述的Y-W方程以及AR”模型與預(yù)測誤差濾波器之間的關(guān)系，就可提

取AR模型參數(shù).目前主要有三種：Levinson-Durbin算法、Burg算法和Marple算法。這三種方法都可以用時間平均代替集合平均的最小平方準(zhǔn)則推導(dǎo)得到。在本文中，我們主要采用TBurg法來估計信號的模型參數(shù)，因此主要介紹一下Burg法。Levinson-Durbin算法L-D遞推算法是在滿足前向預(yù)測均方誤差最小的前提下，先求得觀測數(shù)據(jù)的自相關(guān)函數(shù)，然后利用Y-W方程的遞推性質(zhì)求得模型參數(shù)，進而根據(jù)式(3)求得功率譜的估值.它是模型階次逐次加大的一種算法，即先計算階次m=l時的預(yù)測系數(shù){a(k)}=a⑴和b、再計算m=2時的勺(1)，Of)和b2按此依次計算到階次m=p時的a⑴，止遞推.遞推公式為：{a(k)}=a⑴和b、到階次m=p時的a⑴，止遞推.遞推公式為：…a(p)及b2當(dāng)b2滿足精度要求時即可停pR(m)+2aa(m)=—(k)R(m—k)m—kk=^Em—1(4-6)a(k)=a(k)+a(m)a(m—k)m其中m—1m—1k=1,2,m-1(4-7)2wm1—a(m)21E=r(0)k=1H「1-a(k)(4-8)設(shè)觀察到的N個數(shù)據(jù)為x(0),x(1)，…，設(shè)觀察到的N個數(shù)據(jù)為x(0),x(1)，…，x(N—1)，則b2=R(0)=_w0Nn=0(2)計算反射系數(shù)：p—22e(n)eb(n—1)-1)}—22-—tmm-1]2ef(n)J+m—1n=mb(n)m—1Burg算法的基本思想是直接從觀測的數(shù)據(jù)利用線性預(yù)測器的前向和后向預(yù)測的總均方誤差之和為最小的準(zhǔn)則來估計反射系數(shù)，進而通過L-D算法的遞推公式求出AR模型優(yōu)化參數(shù).具體算法如下：⑴取m=1,初始化:TOC\o"1-5"\h\z€f(n)=eb(n)=x(n),n=0,1N-1⑴取m=1,初始化:001N-12x2(n)(3)計算濾波器系數(shù)及預(yù)測誤差功率：a(m)=pma(k)=a(k)+pa(m—k)k=1,2,m-1EK"(1-P?匕(4)遞推高一階前、后向預(yù)測誤差：ef(n)=ef(n)+Peb(〃-1)mm-1mm-1eb(n)=eb(n-1)+Pef(n)把m更新為m+1,重復(fù)②?④直到。2滿足要求。Marple算法Marple算法又稱為不受約束的最小二乘法，它的主要思路是為了擺脫因采用遞推運算對確定預(yù)測系數(shù)的約束，讓每一預(yù)測系數(shù)(模型參數(shù))的確定直接與前、后向預(yù)測的總的平方誤差最小(最小二乘法)聯(lián)系起來.即令總的平方誤差pm-1m-1n—p無1Xa(k)x(n-k)2+X1Xa(k)x(n+k-p)2a(0)—1n—ppLk—0」n—ppLk—0」p最小.由上式可見，總的平方誤差8是系數(shù)a(k)的函數(shù).若把8對各預(yù)測系數(shù)a(k)而非單一地對a(p)=p求導(dǎo)數(shù)，并令其為零，就可以得到一組線性方程.解此方程組所得的a(k)就是在最小平方誤差準(zhǔn)則下的最優(yōu)預(yù)測系數(shù).但由于方程組系數(shù)矩陣不是Toeplitz型，所以不能利用L—D算法求解.為了減少運算量，Marple提出一種格型結(jié)構(gòu)的高效遞推算法。4.1.1AR模型的穩(wěn)定性及階的確定AR(p)模型穩(wěn)定的充分必要條件是H(z)的極點(即A(z)的根)都在單位圓內(nèi)。穩(wěn)定的AR(p)模型將具有以下的性質(zhì)：H(z)的全部極點或A(z)的所有根都在單位圓內(nèi)。自相關(guān)矩陣是正定的。激勵信號的方差(能量)隨階次的增加而遞。反射系數(shù)的模恒小于1。但是在實際應(yīng)用中，levinson算法的已知數(shù)據(jù)(自相關(guān)值)是由來估計的，有限字長效應(yīng)有可能造成大的誤差，致使估計出來的AR(p)參數(shù)所構(gòu)成的A(z)的根跑到單圓上或單位圓外，從而使模型失去穩(wěn)定。在遞推計算的過程中如果出現(xiàn)這種情況，將致，或|即停止遞推計算。通常事先并不知道AR模型的階。階選得太低，功率譜受到的平滑太厲害，平滑后的譜已經(jīng)分辨不出真實譜中的兩個峰了。階選的太高，固然會提高譜估計的分辨率，但同時會產(chǎn)生虛假的譜峰或譜的細節(jié)。一種簡單而直觀的確定AR模型的階的方法，是不斷增加模型的階，同時觀察預(yù)測誤差功率，當(dāng)其下降到最小時，對應(yīng)的階便可選定為模型的階。但是預(yù)測誤差功率(或AR模型激勵源的方差)是隨著階次的增加而單調(diào)下降的，因此，很難確定降到什么程度才合適。另一方面，應(yīng)該注意到，隨著階次增加，模型參數(shù)的數(shù)目亦增多了，譜估計的方差會變大(表現(xiàn)在虛假譜峰的出現(xiàn))。因此，不能簡單地依靠觀察預(yù)測誤差功率的下降來確定模型的階。與此相應(yīng)的另一種簡單方法是觀察各階模型預(yù)測誤差序列的周期圖，當(dāng)它很接近于平坦(白色譜)時即對應(yīng)最佳的階。4.2MA模型及功率譜估計MA模型是一個全零點模型，其功率譜不易體現(xiàn)信號中的峰值，即分辨率較低。從譜估計的角度來看，MA模型譜估計等效于經(jīng)典譜估計中的自相關(guān)法，若單純?yōu)榱藢σ欢斡邢揲L數(shù)據(jù)做譜估計，就沒必要求解MA模型了。但在系統(tǒng)分析于辨識以及ARMA譜估計中都要用到MA模型，因此仍有必要討論MA系數(shù)的求解方法。目前提出的方法主要有：1，譜分解法；2用高階的AR模型來近似MA模型；3,最大似然法。以第二種方法為例討論MA模型參數(shù)的提取。i有N點數(shù)據(jù)x(n)建立一個p階的AR模型，p>>q，可用AR模型參數(shù)的計算方法求出p階的AR系數(shù)ap(k)，k=1,2...p；^利用a(k)，k=1,2...p進行線性預(yù)測，等效于一個q階的AR模型，再p一次利用AR系數(shù)的求解方法得到b(k)，k=1,2，…，q。進而通過兩次求AR系數(shù)可得MA模型的系數(shù)。4.3ARMA模型理論上可以證明，在一定條件下，ARMA和MA模型都可以用一個階次足夠大的AR模型來近似，所以在實際使用中我們用一個階次比較高的AR模型代替ARMA模型來進行功率譜估計，避免了求解非線性方程組帶來的困難。5現(xiàn)代譜估計的應(yīng)用一基于AR模型的噪聲源識別5.1.1研究背景在噪聲治理實踐中，常常需要鑒別多源噪聲中各聲源對總體噪聲的貢獻情況，以便有針對性地進行治理。但是，在大多數(shù)情況下。不允許或無法使噪聲源設(shè)備單獨運行。例如，油爐燃燒時必須供水，水泵噪聲將與油爐燃燒噪聲并存；發(fā)電廠中所有設(shè)備都不得隨意停機等等。因此，迫切需要能夠進行在線分析各聲源發(fā)聲情況的聲源鑒別系統(tǒng)。本文就是針對上述問題提出的。使用現(xiàn)代信號分析方法處理采集到的儀器設(shè)備在運行過程中產(chǎn)生的噪聲信號，實現(xiàn)對噪聲源的識別。噪聲已被公認(rèn)為是與水質(zhì)污染，大氣污染并列的三大污染之一。在我們生活的環(huán)境中每個人都會受到各種噪聲的干擾和影響。特別是在交通，工業(yè)生產(chǎn)和日常生活中，噪聲更是無處不在。據(jù)全國統(tǒng)計，在反映環(huán)境污染的投訴中，關(guān)于噪聲污染的人民來信和來訪的件數(shù)逐年增加，已從1991年的2.78萬件增加至1995年的3.90萬件，增加了40%以上；而反映噪聲污染問題的投訴占環(huán)境污染投訴的信訪比例則從1991年25%增加到1995年的35.6%，五年中增加10個百分點。這一比例高居各類污染投訴的首位。由于環(huán)境噪聲污染影響范圍較大，近年來因噪聲擾民引起的糾紛不斷出現(xiàn)，其中以反映商業(yè)、飲食服務(wù)業(yè)和建筑施工場所噪聲擾民居多。5.1.2噪聲源識別與控制方法現(xiàn)狀噪聲控制的目的就是根據(jù)實際需要和可能性，用最經(jīng)濟的辦法把噪聲限制在某種合適的范圍內(nèi)。為了有效地控制噪聲，一般根據(jù)噪聲傳播的具體情況，分別在噪聲源部位，噪聲傳播途徑或噪聲接受部位采取措施。目前經(jīng)常采用的措施是對噪聲傳播途徑的控制。即在噪聲傳播的過程中利用吸聲材料或吸收結(jié)構(gòu)來吸收聲能，用隔聲設(shè)旌隔離噪聲，或者用消聲設(shè)備降低空氣中聲的傳播。實際上我們對噪聲傳播途徑的控制也是基于對主次聲源識別的基礎(chǔ)上的。識別了主要噪聲源，我們才能有針對性地治理噪聲。在目前對噪聲源識別的研究中，主要有以下幾種方法：5.1.2.1消去法所謂消去法，就是首先把附帶全部裝備的試驗對象（如汽車、發(fā)動機等）在一定條件下測定其總和的工作噪聲，然后去除其中的一部分裝備或控制這部分噪聲傳出，再按同樣試驗條件測定去除部分裝備后的試驗對象的工作噪聲，則去除這部分裝備前后試驗對象噪聲交化值，即被視作被拆卸或控制傳出裝備噪聲的方法。這一方法試驗難度大，耗時耗力，由于對部分設(shè)備要求停機而不能進行在線測量，因而不能獲得實驗對象噪聲分布的整體概念.5.1.2.2聲強測量法聲強測量法是利用雙點聲壓梯度的積分來近似空氣質(zhì)點的振動速度，并利用FFT來實現(xiàn)聲強的實時測量。采用聲強測量法分析噪聲可獲得豐富的噪聲信息，不但可測出噪聲的大小，而且還能測出聲能的流向。該方法對環(huán)境要求低?，F(xiàn)場識別能力強.但是很強的背景聲會對聲強測試產(chǎn)生不良的影響，產(chǎn)生較大的誤差，并且其測試精度較難控制。5.1.2.3相干函數(shù)法相十函數(shù)法是根據(jù)相關(guān)理論，用互相關(guān)函數(shù)描述被測部件的振動信號與外界噪聲信號之間的關(guān)系.根據(jù)結(jié)果判斷噪聲是否由該部件振動引起，從而判斷出各個聲源.在實際噪聲測量中，常常是復(fù)雜的多輸入單輸出系統(tǒng)。而且輸入噪聲源又可能是相互依賴的，因此這種方法也不適合判斷多輸入單輸出的系統(tǒng)5.1.2.4頻譜分析法頻譜分析法，由于聲源噪聲的形成機理不同。使每個聲源的嗓聲特點有差異，能量分布的頻譜范圍有的位于全頻帶，有的位于中、低頻或高頻。因而，在了解各組成聲源噪聲頻譜特點的情況下，測量出總噪聲頻譜。再取得各部分聲源的頻譜，并將其與總聲源頻譜比較，就可以分析出各組成部分噪聲貢獻幅度。從而找出需控制的主要聲源，以有效降低高能頻率段的噪聲。頻譜分析法中，除應(yīng)用幅值頻譜外，還常應(yīng)用功率譜。但是由于采用了傳統(tǒng)的功率譜估計，所以存在頻率分辨率低，旁瓣泄漏嚴(yán)重的缺點。由于現(xiàn)代譜估計方法為信號建立了一個準(zhǔn)確或近似的模型，從根本上摒棄了對數(shù)據(jù)序列加窗的隱含假設(shè)，因而能夠更好地估計信號功率譜.利用現(xiàn)代譜方法估計功率譜時，不要求噪聲源中只有一個信號為高斯分布，這與盲信號分離有很大區(qū)別（盲信號分離要求源信號中只能有一個信號為高斯分布）。只要利用能量隨頻率的分布（即功率譜）來求取噪聲源的比重?，F(xiàn)代譜估計的上述優(yōu)點為噪聲源識別提供了更為精確的方法。5.2研究內(nèi)容由于缺少實際的噪聲源，在此采用仿真的噪聲信號來模擬噪聲。我們知道AR模型的參數(shù)與其階次有關(guān)。階選得太低，功率譜受到的平滑太厲害，平滑后的譜已經(jīng)分辨不出真實譜中的兩個峰了。階選的太高，固然會提高譜估計的分辨率，但同時會產(chǎn)生虛假的譜峰或譜的細節(jié)。所以必須選擇合適的AR模型的階次。本文采用選用Akaike(alC)準(zhǔn)則作為確定模型階的依據(jù)。將信號x1,x2按不同的比例混合在一起組成混合噪聲以模仿生活中常見的噪聲信號。5.2.1仿真信號的選取N=512;fs=1;a1=5;a2=3;a3=2;w=2火pi/fs;x1=a1*sin(0.5火w*(0:N-1))+a2火cos(w*(0:N-1)).八2+randn(1,N);x2=a2火sin(w*(0:N-1))+a2火cos(0.8火w*(0:N-1)).八2+randn(1,N);y=1:512;plot(y,x1);figure(2);plot(y,x2);save