1、 摘要：信號參數估計是現代信號處理的重要研究內容之一，在頻域中進行傅里葉變換研究信號，是研究確定性信號最簡單且有效的手段，但在現代信號分析中，對于常見的隨機信號，不可能用清楚的數學關系式來描述，其傅里葉變換更不存在，轉而可以利用給定的若干個樣本數據估計來估計信號的參數。本學期在導師的指導下我學習了這門課程，了解到相關的知識，深刻體會了信號參數估計的理論基礎。本文主要介紹我對信號參數估計中的現代譜估計的理解和有關體會。關鍵字：參數估計；隨機信號；譜估計引言： 功率譜估計是隨機信號處理的重要內容，其技術淵源很長，而且在過去的40余年中獲得了飛速的發展。涉及到信號與系統、隨機信號分析、概率統計、矩陣

2、代數等一系列的基礎學科，廣泛應用于人民的日常生活及軍事、工業、農業活動中，是一個具有強大生命力的研究領域。現代譜估計的方法又大致可分為參數模型譜估計和非參數模型譜估計，前者有AR模型、MA模型、ARMA模型、PRONY模型等，后者有最小方差方法、多分量的MUSIC方法等。一 現代譜估計方法的發展1.1功率譜研究的發展過程功率譜估計是數字信號處理的主要內容之一，主要研究信號在頻域中的各種特征，目的是根據有限數據在頻域內提取被淹沒在噪聲中的有用信號。現代譜估計的提出主要是針對經典譜估計（周期圖和自相關法）的分辨率和方差性能不好的問題。1967 年，Burg 提出的最大嫡譜估計，即是朝著高分辨率譜估

3、計所作的最有意義的努力。雖然，Bartlett 在 1948年，Parzem 于 1957 年都曾經建議用自回歸模型做譜估計，但在 Burg 的論文發表之前，都沒有引起注意。現代譜估計的內容極其豐富，涉及的學科及應用領域也相當廣泛，至今，每年都有大量的論文出現。非參數模型譜估計的特點是其模型不是用有限參數來描述，而直接由相關函數序列得到，這種方法能提高低信噪比時的譜分辨率。參數模型譜估計是先根據過程的先驗信息或者一些假定，建立一個數學模型來表示所給定采樣數據的過程，或者選擇一個較好的近似實際模型，而后利用采樣數據序列或者自相關序列，估計該模型的參數，最后把參數代入到該模型對應的理論功率譜表達式

4、，得到所需要的譜估計。1.2功率譜估計應用及用途功率譜估計有著極其廣泛的應用，不僅在認識一個隨機信號時，需要估計它的功率譜。它還被廣泛地應用于各種信號處理中。在信號處理的許多場所，要求預先知道信號的功率譜密度(或自相關函數)。功率譜估計就是通過信號的相關性估計出接受到信號的功率隨頻率的變化關系，實際用途有濾波，信號識別（分析出信號的頻率），信號分離，系統辨識等。譜估計技術是現代信號處理的一個重要部分，還包括空間譜估計，高階譜估計等。維納濾波、卡爾曼濾波，可用于自適應濾波，信號波形預測等（火控系統中的飛機航跡預判）。二現代譜估計現代功率譜估計即參數譜估計方法是通過觀測數據估計參數模型再按照求參數

5、模型輸出功率的方法估計信號功率譜,主要是針對經典譜估計的分辨率低和方差性能不好等問題提出的。常用模型有 ARMA 模型、 AR 模型、 MA 模型,AR模型的正則方程是一組線性方程，而MA和ARMA模型是非線性方程。由于AR模型具有一系列良好的性能，因此被研究最多也得到最廣泛的應用。2.1 AR模型AR 模型又稱為自回歸模型,是一個全極點的模型,可用如下差分方程來表示: （2-1）其中，p是 A R模型的階數， =l ， 2 ， ， p是p階 A R模型的參數將該模型記為AR( p) ，它的系統轉移函數為 (2-2) 在功率譜估計中，若觀測的數據 x(n) 是平穩隨機過程，則該系統的輸入w(

6、n ) 也可認為是平穩的，因 而根據線性系統對平穩隨機信號的響應理論可得觀測數據的功率譜為 (2-3)由式可知，利用A R模型進行功率譜估計的實質是求解模型系數 和的問 題將式 ( 1 )兩端乘以x ( n-m ) 求平均 ( 數學期望) ，可以求得觀測數據的A R( p) 模型參數與自相關函數的關系式為 (2-4)可見， 階 AR模型輸出的相關函數具有遞推的性質， 因而選用 AR模型進行譜估計只需較少的觀測數據將式 ( 4 )寫成矩陣形式得 (2-5) 上式就是著名的Yule-Walker( YW)方程它表明，只要已知觀測數據的自相關數，就能求出AR模型參數 和，進而按式 ( 3 )求得信號

7、功率譜的估值。 另外，從 AR模型的差分方程式可知，該模型的現在輸出值是它本身過去值的回歸，這與預測器存在著一定的相似性，它們之間有著非常密切的關系，即它們的系統函數互為倒數，也就是說預測誤差濾波器 就是 AR信號模型H(z)的逆濾波器因此通過預測誤差濾波器優化設計使預測均方誤差最小就可求得A R信號模型的最優參數，即P階線性預測器的預測系數 等于p階A R模型的系數，其最小均方預測誤差等于自噪聲方差。 因此， 根據上述的Y-W方程以及 AR模型與預測誤差濾波器之間的關系，就可提取 A R模型參數目前主要有三種：Levinson-Durbin算法、Burg算法和Marple算法。這三種方法都可

8、以用時間平均代替集合平均的最小平方準則推導得到。在本文中，我們主要采用TBurg法來估計信號的模型參數，因此主要介紹一下Burg法。Levinson-Durbin算法 L -D 遞推算法是在滿足前向預測均方誤差最小的前提下，先求得觀測數據的自相關函數，然后利用Y - W方程的遞推性質求得模型參數，進而根據式 ( 3 )求得功率譜的估值它是模型階次逐次加大的一種算法，即先計算階次m=l 時的預測系數 = 和，再計算m= 2 時的，和 按此依次計算到 階次m=p時的， 及當滿足精度要求時即可停止遞推遞推公式為： (2-6) k=1,2,m-1 (2-7)其中 （2-8）Burg算法 Burg算法的

9、基本思想是直接從觀測的數據利用線性預測器的前向和后向預測的總均方誤差之和為最小的準則來估計反射系數，進而通過 L -D算法的遞推公式求出 AR模型優化參數設觀察到的N個數據為x( 0 ),x(1)， ， x ( N1 ) ，則具體算法如下： (1)取m=1,初始化：= =x(n), n=0,1N-1 (2)計算反射系數： （3）計算濾波器系數及預測誤差功率： k=1,2,m-1 (4)遞推高一階前、后向預測誤差： 把m更新為m+1 ，重復直到滿足要求。Marple算法 Marple算法又稱為不受約束的最小二乘法，它的主要思路是為了擺脫因采用遞推運算對確定預測系數的約束，讓每一預測系數 ( 模型

10、參數)的確定直接與前、后向預測的總的平方誤差最小 ( 最小二乘法)聯系起來即令總的平方誤差 最小由上式可見，總的平方誤差是系數的函數若把對各預測系數而非單一地對=求導數，并令其為零，就可以得到一組線性方程 解此方程組所得的就是在最小平方誤差準則下的最優預測系數但由于方程組系數矩陣不是Toeplitz型，所以不能利用 L D算法求解為了減少運算量，Ma r p l e 提出一種格型結構的高效遞推算法。2.1.1 AR模型的穩定性及階的確定AR(p)模型穩定的充分必要條件是H(z)的極點(即A(z)的根)都在單位圓內。穩定的AR(p)模型將具有以下的性質：(1) H(z)的全部極點或A(z)的所有

11、根都在單位圓內。(2) 自相關矩陣是正定的。(3) 激勵信號的方差(能量)隨階次的增加而遞。(4) 反射系數的模恒小于1。但是在實際應用中，levinson算法的已知數據(自相關值)是由 來估計的，有限字長效應有可能造成大的誤差，致使估計出來的AR(p)參數所構成的A(z)的根跑到單圓上或單位圓外，從而使模型失去穩定。在遞推計算的過程中如果出現這種情況，將致，或|即停止遞推計算。通常事先并不知道AR模型的階。階選得太低，功率譜受到的平滑太厲害，平滑后的譜已經分辨不出真實譜中的兩個峰了。階選的太高，固然會提高譜估計的分辨率，但同時會產生虛假的譜峰或譜的細節。一種簡單而直觀的確定AR模型的階的方法

12、，是不斷增加模型的階，同時觀察預測誤差功率，當其下降到最小時，對應的階便可選定為模型的階。但是預測誤差功率(或AR模型激勵源的方差)是隨著階次的增加而單調下降的，因此，很難確定降到什么程度才合適。另一方面，應該注意到，隨著階次增加，模型參數的數目亦增多了，譜估計的方差會變大(表現在虛假譜峰的出現)。因此，不能簡單地依靠觀察預測誤差功率的下降來確定模型的階。與此相應的另一種簡單方法是觀察各階模型預測誤差序列的周期圖，當它很接近于平坦(白色譜)時即對應最佳的階。2.2 MA模型及功率譜估計 MA模型是一個全零點模型，其功率譜不易體現信號中的峰值，即分辨率較低。從譜估計的角度來看，MA模型譜估計等效

13、于經典譜估計中的自相關法，若單純為了對一段有限長數據做譜估計，就沒必要求解MA模型了。但在系統分析于辨識以及ARMA譜估計中都要用到MA模型，因此仍有必要討論MA系數的求解方法。目前提出的方法主要有：1，譜分解法；2用高階的AR模型來近似MA模型；3，最大似然法。以第二種方法為例討論MA模型參數的提取。 有N點數據x（n）建立一個p階的AR模型，p>>q，可用AR模型參數的計算方法求出p階的AR系數，k=1,2p； 利用，k=1,2p進行線性預測，等效于一個q階的AR模型，再一次利用AR系數的求解方法得到b（k），k=1,2，q。 進而通過兩次求AR系數可得MA模型的系數。2.3A

14、RMA模型理論上可以證明，在一定條件下，A R MA和 M A模型都可以用一個階次足夠大的A R模型來近似，所以在實際使用中我們用一個階次比較高的A R模型代替A R MA模型來進行功率譜估計，避免了求解非線性方程組帶來的困難。3基于AR模型的噪聲源識別3.1噪聲源識別的方法：3.1.1 消去法所謂消去法，就是首先把附帶全部裝備的試驗對象(如汽車、發動機等)在一定條件下測定其總和的工作噪聲，然后去除其中的一部分裝備或控制這部分噪聲傳出，再按同樣試驗條件測定去除部分裝備后的試驗對象的工作噪聲，則去除這部分裝備前后試驗對象噪聲交化值，即被視作被拆卸或控制傳出裝備噪聲的方法。這一方法試驗難度大，耗時

15、耗力，由于對部分設備要求停機而不能進行在線測量，因而不能獲得實驗對象噪聲分布的整體概念.3.1.2 聲強測量法聲強測量法是利用雙點聲壓梯度的積分來近似空氣質點的振動速度，并利用FFT來實現聲強的實時測量。采用聲強測量法分析噪聲可獲得豐富的噪聲信息，不但可測出噪聲的大小，而且還能測出聲能的流向。該方法對環境要求低。現場識別能力強但是很強的背景聲會對聲強測試產生不良的影響，產生較大的誤差，并且其測試精度較難控制。3.1.3 相干函數法相干函數法是根據相關理論，用互相關函數描述被測部件的振動信號與外界噪聲信號之間的關系根據結果判斷噪聲是否由該部件振動引起，從而判斷出各個聲源在實際噪聲測量中，常常是復

16、雜的多輸入單輸出系統。而且輸入噪聲源又可能是相互依賴的，因此這種方法也不適合判斷多輸入單輸出的系統3.1.4 頻譜分析法頻譜分析法，由于聲源噪聲的形成機理不同。使每個聲源的嗓聲特點有差異，能量分布的頻譜范圍有的位于全頻帶，有的位于中、低頻或高頻。因而，在了解各組成聲源噪聲頻譜特點的情況下，測量出總噪聲頻譜。再取得各部分聲源的頻譜，并將其與總聲源頻譜比較，就可以分析出各組成部分噪聲貢獻幅度。從而找出需控制的主要聲源，以有效降低高能頻率段的噪聲。頻譜分析法中，除應用幅值頻譜外，還常應用功率譜。但是由于采用了傳統的功率譜估計，所以存在頻率分辨率低，旁瓣泄漏嚴重的缺點。由于現代譜估計方法為信號建立了一

17、個準確或近似的模型，從根本上摒棄了對數據序列加窗的隱含假設，因而能夠更好地估計信號功率譜利用現代譜方法估計功率譜時，不要求噪聲源中只有一個信號為高斯分布，這與盲信號分離有很大區別(盲信號分離要求源信號中只能有一個信號為高斯分布)。只要利用能量隨頻率的分布(即功率譜)來求取噪聲源的比重。現代譜估計的上述優點為噪聲源識別提供了更為精確的方法。3.2研究內容由于缺少實際的噪聲源，在此采用仿真的噪聲信號來模擬噪聲。我們知道AR模型的參數與其階次有關。階選得太低，功率譜受到的平滑太厲害，平滑后的譜已經分辨不出真實譜中的兩個峰了。階選的太高，固然會提高譜估計的分辨率，但同時會產生虛假的譜峰或譜的細節。所以

18、必須選擇合適的AR模型的階次。本文采用選用Akaike(alC)準則作為確定模型階的依據。將信號x1,x2按不同的比例混合在一起組成混合噪聲以模仿生活中常見的噪聲信號。3.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 x1 x1;save x2 x

19、2;首先判斷AR模型的階數。這里應用的就是A1C準則。我們用仿真信號x1為例進行了階數測定，與模型階的關系曲線如圖所示。從圖中可以發現模型的AIC存在最小值，此時信號模型對應的階數為12。所以選取12階的模型。clear;load x1.mat;x=x1;N=length(x);M=20;e1=zeros(M,N);e2=zeros(M,N);e1(1,1:N)=x;e2(1,1:N)=x;r=zeros(1,M);a=zeros(1,M);q1=zeros(1,M);q2=zeros(1,M);AIC=zeros(1,M);a(1)=(1/N)*sum(x.2);q1(1)=sum(e1(1

20、,2:N).*e2(1,1:N-1);q2(1)=sum(e1(1,2:N).2)+(e2(1,1:N-1).2);r(1)=2*q1(1)/q2(1);AIC(1)=log(a(1);for k=2:M for n=k:N e1(k,n)=e1(k-1,n)+(-1)*r(k-1)*e2(k-1,n-1); e2(k,n)=(-1)*r(k-1)*e1(k-1,n)+e2(k-1,n-1); end q1(k)=sum(e1(k,k+1:N).*e2(k,k:N-1); q2(k)=sum(e1(k,k+1:N).2)+(e2(k,k:N-1).2); r(k)=2*q1(k)/q2(k);

21、 a(k)=(1-r(k-1)2)*a(k-1); AIC(k)=log(a(k)+k*log(N)/N;endX,Y=min(AIC);k=1:M;plot(k,AIC);X,Y3.2.3功率譜估計及噪聲源識別y1分別取10x1+x2，x1+10x2，x1+x2;我們嘗試著去識別主噪聲源在進行信號混合時利用功率譜估計出哪個信號是主噪聲源，求功率譜的方法如下：function argonglv(y1,y2,y3)clcclearload x1.matload x2.matload x3.maty1=10*x1+x2;k1=sqrt(sum(y1.2);y1=y1./k1p1= arburg(y

22、1,12);la=length(p1);ln=128-la;ts=1/5000;ab=zeros(1,ln);bb=p1 abdw=1/128;w=dw*0:127f=fft(bb)f1=abs(f);f2=f1.2;f3=1./f2;f4=10*log10(f3);plot(w(1:128),f4(1:128),'-kX')hold onfigure (1);y1=x1;k1=sqrt(sum(y1.2);y1=y1./k1p1= arburg(y1,12);la=length(p1);ln=128-la;ts=1/5000;ab=zeros(1,ln);bb=p1 abdw=1/128;w=dw*0:127f=fft(bb)f1=abs(f);f2=f1.