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

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

當(dāng)N—s時(shí)，Bias—0。由于我們假定各段周期圖是相互獨(dú)立的，所以可按式(3-8)得到下式:Vark(Vark(rn)11Var\lxxLM(①)](2-11)由此可見，隨著L的增加Vai^P(w)k下降的，當(dāng)L—s時(shí)，Var^P(w)］—0。因xxxx此Bartlett估計(jì)是一致估計(jì)。比較式(2-10)的eP(w)］與式(2-1)的血(°)］可見在二種情況的估計(jì)量的期望值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愈大說明譜分辨率愈差。一個(gè)固定的記錄長度N，周期圖分段的數(shù)目L愈大將使方差愈小，但M也愈小，因而使Bias今式(2-10)中w(eg)的主瓣寬度大于后式中的主瓣寬度，因而Bias。(W)］，而主瓣愈寬顯然使分辨率愈差。因此Bias可用來說明譜的N偏倚△xxBias由此可見Bartlett法使譜估計(jì)的方差減小是用增加Bias以及降低譜分辨率的代價(jià)換來的。實(shí)際上，當(dāng)N一定時(shí)，Bias與Var的互換性是譜估計(jì)的一個(gè)固有特性。為了說明經(jīng)平均后的周期圖作為功率譜估計(jì)的實(shí)際效果，設(shè)有一零均值高斯分布的隨機(jī)過程，其功率譜密度為■?”)■?”)(1-az-1)(1-az)a=1—*'0.1這一功率譜密度是由一零值高斯分布單位方差的噪聲序列通過一個(gè)其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的曲線比較可見在這種估計(jì)中誤差的大部分起因于Bias而不是方差(因?yàn)槎N情況均為M=8，Bias相仿，誤差也相仿)。M=16,L=2及8的周期圖表示在圖2-3。圖2-3中L不同造成的影響也是明

顯的。但是這二個(gè)曲線的起伏都很大，因此有理由認(rèn)為誤差主要起因于方差。比較M=8與M=16的周期圖可見，M愈大（即§2.3中的N愈大），將使周期圖的起伏愈增快的結(jié)論，在這里也同樣成立。比較圖2-2與圖2-3發(fā)現(xiàn)在這個(gè)例子中最好的選擇是應(yīng)用L=16,M=8的估計(jì)而不是L=8、M=16的估計(jì)，即寧可減小方差，犧牲Bias。在實(shí)際中，當(dāng)然功率譜密度的真值是不知道的。但是譜的窗口函數(shù)以及關(guān)于功率譜密度的某些信息往往是預(yù)函先知道的。通過改變M和L以及利用預(yù)先知道的情況，通?？梢院芎玫剡M(jìn)行選擇。平均周期圖的方法特別適合于應(yīng)用FFT算法。因此在FFT出現(xiàn)以后這個(gè)方法比下面將要討論的利用窗口函數(shù)的處理法用得更多。而在FFT出現(xiàn)以前主要是用窗口函數(shù)處理法平滑周期圖。圖2-1P(o圖2-1P(o)與E[18(0)]的特性圖2-2平滑后的周期圖（每段取8個(gè)數(shù)據(jù)）圖2-3平均后的周期圖（每段取16個(gè)數(shù)據(jù)）2.4Welch法Welch提出了對Bartlett法的修正使更適合于用FFT進(jìn)行計(jì)算。他主要提出二方面的修正，其一是選擇適當(dāng)?shù)拇昂瘮?shù)o⑴，并在周期圖計(jì)算前直接加進(jìn)法，這樣得到的每一段的周期圖為電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)點(diǎn)是無論什么樣的窗函數(shù)均可使譜估計(jì)非負(fù)。其二是在分段時(shí)，可使各段之間有重迭，這樣將會使方差減小(當(dāng)N與M一定時(shí))。他建議可以重迭50%。3經(jīng)典譜估計(jì)方法編程與分析3.1直接法直接法又稱周期圖法，它是把隨機(jī)序列x(n)的N個(gè)觀測數(shù)據(jù)視為一能量有限的序列，直接計(jì)算x(n)的離散傅立葉變換，得X(k),然后再取其幅值的平方，并除以N，作為序列x(n)真實(shí)功率譜的估計(jì)。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功率譜隨著采樣點(diǎn)數(shù)的增加，該股集是漸進(jìn)無騙的。從圖中可以看出，采用周期突發(fā)估計(jì)得出的功率譜很不平滑，相應(yīng)的估計(jì)協(xié)方差比較大。而且采用增加采樣點(diǎn)的辦法也不能吃周期圖變得更加平滑，這是周期圖法的缺點(diǎn)。周期圖法得出的估計(jì)譜方差特性不好：當(dāng)數(shù)據(jù)長度N太大時(shí)，撲線的起伏加??；N太小時(shí)譜的分辨率又不好。對其改進(jìn)的主要方法有二種，即平均和平滑，平均就是將截取的數(shù)據(jù)段七(n)再分成L個(gè)小段，分別計(jì)算功率譜后取功率譜的平均，這種方法使估計(jì)的方差減少，但偏差加大，分辨率下降。平滑是用一個(gè)適當(dāng)?shù)拇昂瘮?shù)W(em)與算出的功率譜5{ejw)進(jìn)行卷積，使譜線平滑。這種方法得出的譜估計(jì)是無偏的，方差也小，但X分辨率下降。3.2間接法間接法先由序列x(n)估計(jì)出自相關(guān)函數(shù)R(n)，然后對R(n)進(jìn)行傅立葉變換，便得到x(n)的功率譜估計(jì)。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');%計(jì)算序列的自相關(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時(shí),BT法與周期圖法估計(jì)出的功率譜是一樣的;當(dāng)M<N-1時(shí),BT法的偏差大于周期圖法，在窗函數(shù)滿足一定條件時(shí)是漸進(jìn)無偏估計(jì);方差小于周期圖的方差;分辨率比周期圖法低,與窗函數(shù)的選擇有關(guān)。BT法的缺點(diǎn)在于當(dāng)M一N時(shí)，R(m)的方差很大,使譜估計(jì)質(zhì)量下降;由R(m)得到的$(明不一定為正值，從而可能失去功率譜的物理意義。3.3改進(jìn)的直接法：對于直接法的功率譜估計(jì)，當(dāng)數(shù)據(jù)長度N太大時(shí)，譜曲線起伏加劇，若N太小，譜的分辨率又不好，因此需要改進(jìn)。3.3.1Bartlett法對周期圖法改進(jìn)的思想是將信號分段進(jìn)行估計(jì)，然后再將這些估計(jì)結(jié)果進(jìn)行平均，從而減小估計(jì)的協(xié)方差，是功率譜圖變得比直接法更平滑。增加分段數(shù)可以進(jìn)一步減低估計(jì)的協(xié)方差，然而若每段中的數(shù)據(jù)點(diǎn)數(shù)太少，就會使估計(jì)的頻率分辨率下降很多。從樣本信號序列總點(diǎn)數(shù)一定的條件下，可以采用使分段相互重疊的方法來增加分段數(shù)，從而保持每段信號點(diǎn)數(shù)不變，這樣就在保證頻率分辨率的前提下進(jìn)一步降低估計(jì)協(xié)方差。而一般的Bartlett法是通過降低分辨率來降低其方差的。Bartlett平均周期圖的方法是將N點(diǎ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法估計(jì)功率譜圖現(xiàn)在比較常用的改進(jìn)方法是Welch法，又叫加權(quán)交疊平均法，簡記為WOSA法，這種方法以加窗（加權(quán)）求取平滑，以分段重疊求得平均，因此集平均與平滑的優(yōu)點(diǎn)于一體，同時(shí)也不可避免帶有兩者的缺點(diǎn)，因此歸根到底是一種折中。其主要步驟是：（1）將N長的數(shù)據(jù)段分成L個(gè)小段，每小段M點(diǎn)，相鄰小段間交疊M/2點(diǎn)（即2:1分段）。因?yàn)長（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（?）這里，財(cái)口代表窗函數(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^進(jìn)亍了兩方面的修正，一是選擇適當(dāng)?shù)拇昂瘮?shù)w(n),并再周期圖計(jì)算前直接加進(jìn)去，加窗的優(yōu)點(diǎn)是無論什么樣的窗函數(shù)均可使譜估計(jì)非負(fù)。二是在分段時(shí)，可使各段之間有重疊，這樣會使方差減小。不同窗函數(shù)的Welch譜估計(jì)在選擇窗函數(shù)時(shí)，一般有如下要求：窗口寬度M要遠(yuǎn)小于樣本序列長度N，以排除不可靠的自相關(guān)值；當(dāng)平穩(wěn)信號為實(shí)過程時(shí)，為保證平滑周期圖和真是功率譜也是實(shí)偶函數(shù)，平滑窗函數(shù)必須是實(shí)偶對稱的；平滑窗函數(shù)應(yīng)當(dāng)在m=0出游峰值，并且m隨絕對值增加而單調(diào)下降，使可靠的自相關(guān)值有較大的權(quán)值；功率譜是頻率的非負(fù)函數(shù)且周期圖是非負(fù)的，因而要求窗函數(shù)的Fourier變換是非負(fù)的。仍以上述平穩(wěn)隨機(jī)信號為例，采樣頻率、采樣點(diǎn)數(shù)、FFT點(diǎn)數(shù)、窗長度及重疊數(shù)據(jù)不變，窗函數(shù)采用矩形窗、Blackman窗、Hamming窗，仿真結(jié)果如圖所示。由仿真結(jié)果可以看出：使用不同的窗函數(shù)譜估計(jì)的質(zhì)量是不一樣的，矩形窗的主瓣較窄，分辨率較好，但方差較大，噪聲水平較高；而Blackman窗和Hamming窗的主瓣較寬，分辨率較低，但方差較小，噪聲水平較低。因此，在進(jìn)行譜分析時(shí)選擇何種窗函數(shù)，要視具體情況而定。如果強(qiáng)調(diào)高分辨率，能精確讀出主瓣頻率，而不關(guān)心幅度的精度，例如測量震動(dòng)物體的自震頻率，可以選用主瓣寬度比較窄的矩形窗；對受到強(qiáng)干擾的窄帶信號若干擾靠近信號，則可選用旁瓣幅度較小的窗函數(shù)，若離開通帶較遠(yuǎn)，則可選用漸近線衰減速度比較快的窗函數(shù)。總之，要針對不同的信號和不同的處理目的來選擇合適的窗函數(shù)，這樣才能得到良好的效果。3.4經(jīng)典譜估計(jì)方法的比較總結(jié)周期圖法和BT法的特點(diǎn)是離散性大，曲線粗糙，方差較大，但是分辨率較高；Bartlett法和Welch法的收斂性較好，曲線平滑，方差較小，但是功率譜主瓣較寬，分辨率低，這是由于對隨機(jī)序列加窗截?cái)嗨鸬腉ibbs效應(yīng)造成的；(3)與Bartlett法相比，Welch法的估計(jì)曲線比較粗糙，但是分辨率較好，原因是Welch法中對數(shù)據(jù)進(jìn)行截?cái)鄷r(shí)加的是Hanning窗，而在Bartlett法中使用的是矩形窗，相對于矩形窗，Hanning窗的主瓣包含更多的能量，因而使功率譜的主瓣較窄，分辨率較高。由上述理論分析及仿真實(shí)驗(yàn)可知，Welch法采用加窗交疊求功率譜，可以有效減小方差和偏差，一般情況下能接近一致估計(jì)的要求，因而得到廣泛應(yīng)用。同時(shí)還可以發(fā)現(xiàn)，對信號加不同的窗函數(shù)，譜估計(jì)的質(zhì)量是不同的。在經(jīng)典譜估計(jì)中，無論是周期圖法還是其改進(jìn)的方法，都存在著頻率分辨率低、方差性能不好的問題，原因是譜估計(jì)時(shí)需要對數(shù)據(jù)加窗截?cái)?，用有限個(gè)數(shù)據(jù)或其自相關(guān)函數(shù)來估計(jì)無限個(gè)數(shù)據(jù)的功率譜，這其實(shí)是假定了窗以外的數(shù)據(jù)或自相關(guān)函數(shù)全為零，這種假定是不符合實(shí)際的，正是由于這些不符合實(shí)際的假設(shè)造成了經(jīng)典譜估計(jì)分辨率較差。另外，經(jīng)典譜估計(jì)的功率譜定義中既無求均值運(yùn)算又無求極限運(yùn)算，因而使得譜估計(jì)的方差性能較差，當(dāng)數(shù)據(jù)很短時(shí)，這個(gè)問題更為突出。如何選取最佳窗函數(shù)、提高頻譜分辨率，如何在短數(shù)據(jù)情況下提高信號譜估計(jì)質(zhì)量，還需要進(jìn)一步研究。4.現(xiàn)代譜估計(jì)現(xiàn)代功率譜估計(jì)即參數(shù)譜估計(jì)方法是通過觀測數(shù)據(jù)估計(jì)參數(shù)模型再按照求參數(shù)模型輸出功率的方法估計(jì)信號功率譜，主要是針對經(jīng)典譜估計(jì)的分辨率低和方差性能不好等問題提出的。常用模型有ARMA模型、AR模型、MA模型,AR模型的正則方程是一組線性方程，而MA和ARMA模型是非線性方程。由于AR模型具有一系列良好的性能，因此被研究最多也得到最廣泛的應(yīng)用。本節(jié)將較為詳細(xì)的討論AR模型，并對MA和ARMA模型譜估計(jì)方法做簡要的討論。4.1AR模型AR模型又稱為自回歸模型，是一個(gè)全極點(diǎn)的模型,可用如下差分方程來表示：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在功率譜估計(jì)中，若觀測的數(shù)據(jù)x(n)是平穩(wěn)隨機(jī)過程，則該系統(tǒng)的輸入w(n)也可認(rèn)為是平穩(wěn)的，因而根據(jù)線性系統(tǒng)對平穩(wěn)隨機(jī)信號的響應(yīng)理論可得觀測數(shù)據(jù)的功率譜為(4-3)P(w)=b2ih(ejw)|2=ZWxw11+2^ae~jkwk(4-2)(4-3)k=1由式可知，利用ar模型進(jìn)行功率譜估計(jì)的實(shí)質(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模型進(jìn)行譜估計(jì)只需較少的觀測數(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)進(jìn)而按式(3)求得信號功率譜的估上式就是著名的Yule-Walker(Y—W)方程.它表明，只要已知觀測數(shù)據(jù)的自相關(guān)數(shù)，就能求出AR模型參數(shù)｛ak｝進(jìn)而按式(3)求得信號功率譜的估該模型的現(xiàn)在輸出值是它本身過去

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