周 云,李世平,羅 鵬,周秋平
(1.第二炮兵工程學院,陜西 西安 710025;2.重慶城市職業學院,重慶 402160)
人們為了方便分析,將信號簡化為只有一個恒定頻率的平穩信號,然而在實際生活中遇到的大多數信號都是頻率隨時間變化的非平穩信號,因此瞬時頻率的概念就顯得尤為重要。
瞬時頻率的定義最早是由Car son和Fry在研究調頻信號時分別提出的,在Gabor提出了解析信號的概念之后,Ville將二者結合起來,提出了現在普遍接受的實信號的瞬時頻率的定義,即:實信號的瞬時頻率就是該信號所對應的解析信號的相位關于時間的導數[1]。該定義只對單分量信號有意義,針對單分量信號求解瞬時頻率目前有多種方法:相位法、譜峰檢測法[1]、過零點法、求根估計算法[2]以及Teager能量算子法[3]等。按照上述定義,文獻[4]指出,要使多分量信號的瞬時頻率有意義,對信號的要求十分苛刻,只有將多分量信號分解為單分量信號后瞬時頻率才有意義。因此,基于HHT、局部均值等方法求解瞬時頻率的最終落腳點都在求解單分量信號的瞬時頻率。該文的研究對象就是非平穩的單分量信號。
Hilbert變換求解瞬時頻率是相位法求解瞬時頻率的一種方法[5]。實信號經過Hilbert變換后變換為解析信號,解析信號是一種復信號,將解析信號的虛部除以實部就得到信號的相位,將信號的相位對時間求導即可得到解析信號的瞬時頻率。解析信號的瞬時頻率即是信號的瞬時頻率[1]。
對某一實信號s(t)經Hilbert變換得到解析信號為
其中:A(t)=[s(t)2+H(s(t))2]1/2
瞬時頻率為
通過Hilbert變換求解瞬時頻率產生誤差的主要原因在于Hilbert變換是通過傅里葉變換得到,在傅里葉變換過程中會產生頻譜泄漏和柵欄效應從而使信號失真。并且只有幅值A(t)和相位φ(t)在頻域完全分開而沒有重疊時將信號變換為解析信號才能表達信號的真實物理意義,從而用式(3)準確求出瞬時頻率。
W-V分布是譜峰檢測法中的一種。1932年,Wigner提出了Wigner分布,最初應用于量子力學研究。1948年,Ville將其引入信號分析領域。1980年,Claasen和Mecklenbraker聯合發表的論文詳細闡述了Wigner-Ville分布的概念、定義、性質以及數值計算等問題。
Wigner-Ville分布的定義為
Wigner-Ville分布精確地定位了信號s(t)的瞬時頻率:
Wigner-Ville分布滿足許多優良的時頻分布數學特性,如邊緣特性,實值性,時、頻移不變性,一致性等。但Wigner-Ville分布的不足在于不能保證非負性,尤其是對多分量信號會產生嚴重的交叉項干擾,使得2個單分量信號在時頻平面上相距很遠。不包含交叉干擾且具有Wigner-Ville分布聚集性的時頻分布是不存在的,如何降低交叉項干擾產生的誤差進而提高Wigner-Ville分布估計瞬時頻率的準確度是研究的熱點問題,可以通過對時間和頻率分別加窗的方法來減少交叉干擾[1,6]。同時,為了消除更多的噪聲干擾,可對檢測信號進行多次峰值檢測,即對第1次得到的估計信號進行重構,然后對重構后的信號重復進行峰值檢測。經仿真驗證可知,此迭代方法能夠提高瞬時頻率的估計準確度[7]。
大量仿真實驗證明,AR模型對平穩信號的測頻準確度高于Hilbert變換和W-V變換。可以設想:將AR模型用于非平穩信號瞬時頻率的測量能否有較高的測頻準確度?但AR模型不能處理非平穩信號,而其改進算法時變AR模型卻能求解非平穩信號的瞬時頻率。
時變AR模型求解瞬時頻率的一種方法是通過求根估計算法求解瞬時頻率。首先,通過時變AR模型求解信號的參數模型;其次,求解參數模型的根;最后,尋找與單位圓最近的根,求解這個根的頻率值。距離單位圓最近根的頻率值即是此時刻的頻率。通過時變AR模型求解瞬時頻率的另一方法是求出功率譜峰值所對應的頻率和時間t,這樣即可得到信號的頻率-時間變化圖。本文采用前一種方法。
設一非平穩過程的觀測值 s(t),t=1,2,…,N 滿足時變AR模型的條件,則:
式中:p——模型階數;
ai(t)——t時刻模型參數;
ν(t)——服從N(0,1)且與s(t-i)相互獨立的高斯白噪聲殘差序列。
Grenier提出的時間基擴展方法是對ai(t)參數辨識最為常用的方法,在式(6)的基礎上進一步假設ai(t)可分為時間基函數的線性組合:
式中:gj(t)——基函數,可以為多項式基、切比雪夫基、傅里葉基、離散余弦基函數[8-10]。
這樣就可以將時變的模型轉化為時不變模型。
利用前后向無約束最小二乘法[2,8,11],前向預測值為
將式(6)與式(8)相減得前向預測誤差為
從而,前向預測誤差能量為
僅當p+1≤t≤N時,式(10)才包括所有時變AR(p)模型參數。相應地,后向預測誤差為
后向預測誤差能量為
僅當 1≤t≤N-p時,式(12)才包括所有時變AR(p)模型參數。總的預測誤差能量為
將式(7)、式(9)~式(12)帶入式(13)并求解方程組即可求出參數,具體求解方程組過程可參見文獻[12]。
前面所述的時變AR模型求解是建立在階數p已知的情況下。然而,實際問題中,階數p往往是未知的或者根本就不存在確定的階數p,所以要考慮階數p的確定。
目前,確定AR模型階數p的方法包括AIC準則、貝葉斯方法、最大似然估計方法等[2,8]。
求解多項式
式中:z=ejw,得到根為 z1,z2,…,zp。
由fi(t)=angle(zi(t))·Fs/(2π)即可得到瞬時頻率,其中 zi為 z1,z2,…,zp中距離單位圓最近的值,Fs為采樣頻率。
設信號s(t)由2個分段線性調頻信號組成
則理論上它的瞬時頻率f(t)為
圖1為信號的理想瞬時頻率,在0<t<0.256s時信號的頻率單調遞增,在0.256s<t<0.512s時信號單調遞減,在0.256s處的頻率為299.8Hz。圖2為使用Hilbert變換求解并歸一化后的信號瞬時頻率,從圖2可以看出在信號兩端及轉折點處都出現了較大的失真,而在中間部分也有震蕩情況發生,這與Hilbert變換本身存在的邊界效應有關。圖3為使用W-V分布求解并歸一化后的信號瞬時頻率,從圖中可以得知可以基本分離出信號瞬時頻率的變化趨勢,具有較好的時頻聚集度,但存在較嚴重的交叉項干擾。文獻[6]證明了不含交叉項干擾且具有W-V分布聚集性的時頻分布是不存在的,并且不含交叉項干擾而聚集性充分接近W-V分布聚集性的時頻分布一般也不存在。圖4為使用時變AR模型求解的信號瞬時頻率,采樣頻率Fs=1kHz,模型階數p=4,采用gj(t)=tj的基時間函數進行分析。從圖4中可以看出,除了在0.256s附近出現了一定偏差外,其他地方效果都不錯,與圖2和圖3相比則有較大的改善。
圖1 信號理想頻率
圖2 Hilbert變換求解瞬時頻率
圖3 W-V分布求解瞬時頻率
圖4 時變AR模型求解瞬時頻率
本文分析了瞬時頻率測量的重要意義、發展現狀和目前主要研究方法,通過Hilbert變換、W-V分布和時變AR模型3種方法實現瞬時頻率的測量,并比較了3種方法的測量準確度和優缺點。
通過理論推導和實驗仿真可以得出:對于單分量信號,時變AR模型測量瞬時頻率具有最高的準確度;Hilbert變換具有最明確的物理意義,是最直觀、最易讓人理解的求解瞬時頻率方法,但在兩端存在較大的誤差;W-V分布具有最高的時頻聚集度,但是存在不可避免的交叉相干擾。時變AR模型測量瞬時頻率具有的優勢將吸引更多的學者進行研究,并解決其存在的階數難以確定等問題,在以后的應用中有較好的發展前景。
[1]張賢達,保錚.非平穩信號分析與處理[M].北京:國防工業出版社,1998.
[2]王宏禹.非平穩隨機信號分析與處理[M].北京:國防工業出版社,1999.
[3]樊新海,何嘉武,王戰軍,等.基于Teager能量算子的解調方法[J].裝甲兵工程學院學報,2009,23(4):60-64.
[4]鄒巖崑,馬孝江,朱泓,等.論局域波法中瞬時頻率分析[J].大連理工大學學報,2003,43(6):775-779.
[5]劉小丹,孫曉奇,沈濱.基于EMD的信號瞬時頻率估計[J].遼寧師范大學學報:自然科學版,2009,32(3):50-57.
[6]鄒紅星,戴瓊海,李衍達,等.不含交叉項干擾且具有WVD聚集時頻分布之不存在性[J].中國科學(E輯),2001,31(4):348-355.
[7]趙宏偉.基于時頻分布的跳頻信號參數檢測方法研究[M].西安:西北工業大學,2006.
[8]史杰,鐘雄虎,宋叔飚,等.基于前后向線性預測的時變參數估計方法研究[J].計算機測量與控制,2006,14(11):1515-1518.
[9]劉影,謝弛,楊府.激光跟蹤儀系統測量不確定度預測方法研究[J].中國測試,2011,36(4):25-27.
[10]鄧衛強,王躍鋼,楊穎濤,等.時變自回歸模型參數拉格朗日乘子估計方法[J].探測與控制學報,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]王永民,何幼樺,忻莉莉,等.時變自回歸模型系數的估計及預測[J].應用數學與計算數學學報,2007,21(2):35-41.