基于時(shí)變AR模型的瞬時(shí)頻率測量

周 云，李世平，羅 鵬，周秋平

（1.第二炮兵工程學(xué)院，陜西 西安 710025；2.重慶城市職業(yè)學(xué)院，重慶 402160）

0 引 言

人們?yōu)榱朔奖惴治?，將信號簡化為只有一個(gè)恒定頻率的平穩(wěn)信號，然而在實(shí)際生活中遇到的大多數(shù)信號都是頻率隨時(shí)間變化的非平穩(wěn)信號，因此瞬時(shí)頻率的概念就顯得尤為重要。

瞬時(shí)頻率的定義最早是由Car son和Fry在研究調(diào)頻信號時(shí)分別提出的，在Gabor提出了解析信號的概念之后，Ville將二者結(jié)合起來，提出了現(xiàn)在普遍接受的實(shí)信號的瞬時(shí)頻率的定義，即：實(shí)信號的瞬時(shí)頻率就是該信號所對應(yīng)的解析信號的相位關(guān)于時(shí)間的導(dǎo)數(shù)[1]。該定義只對單分量信號有意義，針對單分量信號求解瞬時(shí)頻率目前有多種方法：相位法、譜峰檢測法[1]、過零點(diǎn)法、求根估計(jì)算法[2]以及Teager能量算子法[3]等。按照上述定義，文獻(xiàn)[4]指出，要使多分量信號的瞬時(shí)頻率有意義，對信號的要求十分苛刻，只有將多分量信號分解為單分量信號后瞬時(shí)頻率才有意義。因此，基于HHT、局部均值等方法求解瞬時(shí)頻率的最終落腳點(diǎn)都在求解單分量信號的瞬時(shí)頻率。該文的研究對象就是非平穩(wěn)的單分量信號。

1 Hilbert變換求瞬時(shí)頻率

Hilbert變換求解瞬時(shí)頻率是相位法求解瞬時(shí)頻率的一種方法[5]。實(shí)信號經(jīng)過Hilbert變換后變換為解析信號，解析信號是一種復(fù)信號，將解析信號的虛部除以實(shí)部就得到信號的相位，將信號的相位對時(shí)間求導(dǎo)即可得到解析信號的瞬時(shí)頻率。解析信號的瞬時(shí)頻率即是信號的瞬時(shí)頻率[1]。

對某一實(shí)信號s（t）經(jīng)Hilbert變換得到解析信號為

其中：A（t）=[s（t）2+H（s（t））2]1/2

瞬時(shí)頻率為

通過Hilbert變換求解瞬時(shí)頻率產(chǎn)生誤差的主要原因在于Hilbert變換是通過傅里葉變換得到，在傅里葉變換過程中會產(chǎn)生頻譜泄漏和柵欄效應(yīng)從而使信號失真。并且只有幅值A(chǔ)（t）和相位φ（t）在頻域完全分開而沒有重疊時(shí)將信號變換為解析信號才能表達(dá)信號的真實(shí)物理意義，從而用式（3）準(zhǔn)確求出瞬時(shí)頻率。

2 W-V分布

W-V分布是譜峰檢測法中的一種。1932年，Wigner提出了Wigner分布，最初應(yīng)用于量子力學(xué)研究。1948年，Ville將其引入信號分析領(lǐng)域。1980年，Claasen和Mecklenbraker聯(lián)合發(fā)表的論文詳細(xì)闡述了Wigner-Ville分布的概念、定義、性質(zhì)以及數(shù)值計(jì)算等問題。

Wigner-Ville分布的定義為

Wigner-Ville分布精確地定位了信號s（t）的瞬時(shí)頻率：

Wigner-Ville分布滿足許多優(yōu)良的時(shí)頻分布數(shù)學(xué)特性，如邊緣特性，實(shí)值性，時(shí)、頻移不變性，一致性等。但Wigner-Ville分布的不足在于不能保證非負(fù)性，尤其是對多分量信號會產(chǎn)生嚴(yán)重的交叉項(xiàng)干擾，使得2個(gè)單分量信號在時(shí)頻平面上相距很遠(yuǎn)。不包含交叉干擾且具有Wigner-Ville分布聚集性的時(shí)頻分布是不存在的，如何降低交叉項(xiàng)干擾產(chǎn)生的誤差進(jìn)而提高Wigner-Ville分布估計(jì)瞬時(shí)頻率的準(zhǔn)確度是研究的熱點(diǎn)問題，可以通過對時(shí)間和頻率分別加窗的方法來減少交叉干擾[1，6]。同時(shí)，為了消除更多的噪聲干擾，可對檢測信號進(jìn)行多次峰值檢測，即對第1次得到的估計(jì)信號進(jìn)行重構(gòu)，然后對重構(gòu)后的信號重復(fù)進(jìn)行峰值檢測。經(jīng)仿真驗(yàn)證可知，此迭代方法能夠提高瞬時(shí)頻率的估計(jì)準(zhǔn)確度[7]。

3 時(shí)變AR模型求解瞬時(shí)頻率

大量仿真實(shí)驗(yàn)證明，AR模型對平穩(wěn)信號的測頻準(zhǔn)確度高于Hilbert變換和W-V變換?？梢栽O(shè)想：將AR模型用于非平穩(wěn)信號瞬時(shí)頻率的測量能否有較高的測頻準(zhǔn)確度？但AR模型不能處理非平穩(wěn)信號，而其改進(jìn)算法時(shí)變AR模型卻能求解非平穩(wěn)信號的瞬時(shí)頻率。

時(shí)變AR模型求解瞬時(shí)頻率的一種方法是通過求根估計(jì)算法求解瞬時(shí)頻率。首先，通過時(shí)變AR模型求解信號的參數(shù)模型；其次，求解參數(shù)模型的根；最后，尋找與單位圓最近的根，求解這個(gè)根的頻率值。距離單位圓最近根的頻率值即是此時(shí)刻的頻率。通過時(shí)變AR模型求解瞬時(shí)頻率的另一方法是求出功率譜峰值所對應(yīng)的頻率和時(shí)間t，這樣即可得到信號的頻率-時(shí)間變化圖。本文采用前一種方法。

3.1 時(shí)變AR模型建模

設(shè)一非平穩(wěn)過程的觀測值 s（t），t=1，2，…，N 滿足時(shí)變AR模型的條件，則：

式中：p——模型階數(shù)；

ai（t）——t時(shí)刻模型參數(shù)；

ν（t）——服從N（0，1）且與s（t-i）相互獨(dú)立的高斯白噪聲殘差序列。

Grenier提出的時(shí)間基擴(kuò)展方法是對ai（t）參數(shù)辨識最為常用的方法，在式（6）的基礎(chǔ)上進(jìn)一步假設(shè)ai（t）可分為時(shí)間基函數(shù)的線性組合:

式中：gj（t）——基函數(shù)，可以為多項(xiàng)式基、切比雪夫基、傅里葉基、離散余弦基函數(shù)[8-10]。

這樣就可以將時(shí)變的模型轉(zhuǎn)化為時(shí)不變模型。

3.2 時(shí)變AR模型求解

利用前后向無約束最小二乘法[2，8，11]，前向預(yù)測值為

將式（6）與式（8）相減得前向預(yù)測誤差為

從而，前向預(yù)測誤差能量為

僅當(dāng)p+1≤t≤N時(shí)，式（10）才包括所有時(shí)變AR（p）模型參數(shù)。相應(yīng)地，后向預(yù)測誤差為

后向預(yù)測誤差能量為

僅當(dāng) 1≤t≤N-p時(shí)，式（12）才包括所有時(shí)變AR（p）模型參數(shù)?？偟念A(yù)測誤差能量為

將式（7）、式（9）～式（12）帶入式（13）并求解方程組即可求出參數(shù)，具體求解方程組過程可參見文獻(xiàn)[12]。

3.3 時(shí)變AR模型定階

前面所述的時(shí)變AR模型求解是建立在階數(shù)p已知的情況下。然而，實(shí)際問題中，階數(shù)p往往是未知的或者根本就不存在確定的階數(shù)p，所以要考慮階數(shù)p的確定。

目前，確定AR模型階數(shù)p的方法包括AIC準(zhǔn)則、貝葉斯方法、最大似然估計(jì)方法等[2，8]。

3.4 由時(shí)變AR模型求解瞬時(shí)頻率

求解多項(xiàng)式

式中：z=ejw，得到根為 z1，z2，…，zp。

由fi（t）=angle（zi（t））·Fs/（2π）即可得到瞬時(shí)頻率，其中 zi為 z1，z2，…，zp中距離單位圓最近的值，F(xiàn)s為采樣頻率。

4 計(jì)算機(jī)仿真例

設(shè)信號s（t）由2個(gè)分段線性調(diào)頻信號組成

則理論上它的瞬時(shí)頻率f（t）為

圖1為信號的理想瞬時(shí)頻率，在0＜t＜0.256s時(shí)信號的頻率單調(diào)遞增，在0.256s＜t＜0.512s時(shí)信號單調(diào)遞減，在0.256s處的頻率為299.8Hz。圖2為使用Hilbert變換求解并歸一化后的信號瞬時(shí)頻率，從圖2可以看出在信號兩端及轉(zhuǎn)折點(diǎn)處都出現(xiàn)了較大的失真，而在中間部分也有震蕩情況發(fā)生，這與Hilbert變換本身存在的邊界效應(yīng)有關(guān)。圖3為使用W-V分布求解并歸一化后的信號瞬時(shí)頻率，從圖中可以得知可以基本分離出信號瞬時(shí)頻率的變化趨勢，具有較好的時(shí)頻聚集度，但存在較嚴(yán)重的交叉項(xiàng)干擾。文獻(xiàn)[6]證明了不含交叉項(xiàng)干擾且具有W-V分布聚集性的時(shí)頻分布是不存在的，并且不含交叉項(xiàng)干擾而聚集性充分接近W-V分布聚集性的時(shí)頻分布一般也不存在。圖4為使用時(shí)變AR模型求解的信號瞬時(shí)頻率，采樣頻率Fs=1kHz，模型階數(shù)p=4，采用gj（t）=tj的基時(shí)間函數(shù)進(jìn)行分析。從圖4中可以看出，除了在0.256s附近出現(xiàn)了一定偏差外，其他地方效果都不錯(cuò)，與圖2和圖3相比則有較大的改善。

圖1 信號理想頻率

圖2 Hilbert變換求解瞬時(shí)頻率

圖3 W-V分布求解瞬時(shí)頻率

圖4 時(shí)變AR模型求解瞬時(shí)頻率

5 結(jié)束語

本文分析了瞬時(shí)頻率測量的重要意義、發(fā)展現(xiàn)狀和目前主要研究方法，通過Hilbert變換、W-V分布和時(shí)變AR模型3種方法實(shí)現(xiàn)瞬時(shí)頻率的測量，并比較了3種方法的測量準(zhǔn)確度和優(yōu)缺點(diǎn)。

通過理論推導(dǎo)和實(shí)驗(yàn)仿真可以得出：對于單分量信號，時(shí)變AR模型測量瞬時(shí)頻率具有最高的準(zhǔn)確度；Hilbert變換具有最明確的物理意義，是最直觀、最易讓人理解的求解瞬時(shí)頻率方法，但在兩端存在較大的誤差；W-V分布具有最高的時(shí)頻聚集度，但是存在不可避免的交叉相干擾。時(shí)變AR模型測量瞬時(shí)頻率具有的優(yōu)勢將吸引更多的學(xué)者進(jìn)行研究，并解決其存在的階數(shù)難以確定等問題，在以后的應(yīng)用中有較好的發(fā)展前景。

[1]張賢達(dá)，保錚.非平穩(wěn)信號分析與處理[M].北京：國防工業(yè)出版社，1998.

[2]王宏禹.非平穩(wěn)隨機(jī)信號分析與處理[M].北京：國防工業(yè)出版社，1999.

[3]樊新海，何嘉武，王戰(zhàn)軍，等.基于Teager能量算子的解調(diào)方法[J].裝甲兵工程學(xué)院學(xué)報(bào)，2009，23（4）：60-64.

[4]鄒巖崑，馬孝江，朱泓，等.論局域波法中瞬時(shí)頻率分析[J].大連理工大學(xué)學(xué)報(bào)，2003，43（6）：775-779.

[5]劉小丹，孫曉奇，沈?yàn)I.基于EMD的信號瞬時(shí)頻率估計(jì)[J].遼寧師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，32（3）：50-57.

[6]鄒紅星，戴瓊海，李衍達(dá)，等.不含交叉項(xiàng)干擾且具有WVD聚集時(shí)頻分布之不存在性[J].中國科學(xué)（E輯），2001，31（4）：348-355.

[7]趙宏偉.基于時(shí)頻分布的跳頻信號參數(shù)檢測方法研究[M].西安：西北工業(yè)大學(xué)，2006.

[8]史杰，鐘雄虎，宋叔飚，等.基于前后向線性預(yù)測的時(shí)變參數(shù)估計(jì)方法研究[J].計(jì)算機(jī)測量與控制，2006，14（11）：1515-1518.

[9]劉影，謝弛，楊府.激光跟蹤儀系統(tǒng)測量不確定度預(yù)測方法研究[J].中國測試，2011，36（4）：25-27.

[10]鄧衛(wèi)強(qiáng)，王躍鋼，楊穎濤，等.時(shí)變自回歸模型參數(shù)拉格朗日乘子估計(jì)方法[J].探測與控制學(xué)報(bào)，2010，32（5）：59-64.

[11]Reyes B A，Villalobos S C，Camarena R G，et al.Time-frequency representations for second heart sound analysis[C]∥30th Annual International IEEE EMBS Conference：Vancouver，2008（8）：3616-3620.

[12]王永民，何幼樺，忻莉莉，等.時(shí)變自回歸模型系數(shù)的估計(jì)及預(yù)測[J].應(yīng)用數(shù)學(xué)與計(jì)算數(shù)學(xué)學(xué)報(bào)，2007，21（2）：35-41.