信號參數(shù)估計

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

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

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

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

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

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

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

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

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

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

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

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

13、于經(jīng)典譜估計中的自相關法，若單純?yōu)榱藢σ欢斡邢揲L數(shù)據(jù)做譜估計，就沒必要求解MA模型了。但在系統(tǒng)分析于辨識以及ARMA譜估計中都要用到MA模型，因此仍有必要討論MA系數(shù)的求解方法。目前提出的方法主要有：1，譜分解法；2用高階的AR模型來近似MA模型；3，最大似然法。以第二種方法為例討論MA模型參數(shù)的提取。 有N點數(shù)據(jù)x（n）建立一個p階的AR模型，pq，可用AR模型參數(shù)的計算方法求出p階的AR系數(shù)，k=1,2p； 利用，k=1,2p進行線性預測，等效于一個q階的AR模型，再一次利用AR系數(shù)的求解方法得到b（k），k=1,2，q。 進而通過兩次求AR系數(shù)可得MA模型的系數(shù)。2.3ARMA模型理論上

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

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

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

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

18、R模型的階次。本文采用選用Akaike(alC)準則作為確定模型階的依據(jù)。將信號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 x2;首先判斷AR

19、模型的階數(shù)。這里應用的就是A1C準則。我們用仿真信號x1為例進行了階數(shù)測定，與模型階的關系曲線如圖所示。從圖中可以發(fā)現(xiàn)模型的AIC存在最小值，此時信號模型對應的階數(shù)為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,2:N).*e

20、2(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); a(k)=(1

21、-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(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.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 ab