基于廣義解調和廣義S變換的時頻域去噪方法

柏 晨,肖澤龍,許建中

(南京理工大學電子工程與光電技術學院,江蘇南京 210094)

0 引言

利用電磁波對膛內彈丸運動參數測試是目前比較主要的測量方式。然而,在利用該方法進行測試時,由于受到信號選擇、反射面確定、電離干擾及外界環境的影響,得到的干涉信號中存在大量的噪聲,信號信噪比差。因而,要得到高精度的測試結果,就對后續的數據處理提出了較高的要求[1]。

為了改善信噪比,可以引入小波分解的方法對干涉信號進行預處理以抑制部分噪聲,隨后利用短時傅里葉變換在時頻域對其進行分析[1-2]。而短時傅里葉變換的窗函數固定,對快速變化的信號,其時頻分辨率難以得到理想的結果。為了解決這一問題,文獻[3]提出利用希爾伯特變換的方法提取信號的相位信息,從相位信息中得到彈丸的運動表達式。但此方法仍是基于信號信噪比較好的情況下得出的。并且與EMD分解相比,小波分解缺乏自適應性,分解后的各個分量仍有可能是多分量信號。文獻[4]在對低信噪比信號研究分析的基礎上,提出了一種基于希爾伯特黃變換的分析處理方法。但是抑噪效果并不是非常理想,處理后信號的時頻脊線圖中仍包含有大量奇異點,必然使得計算結果誤差增大。

由于傳統的時域或頻域去噪方法經常受到使用條件的限制,本文提出了基于廣義解調和廣義S變換的時頻去噪方法。

1 理論基礎

1.1 廣義解調

信號x(t)的廣義傅里葉變換(Generalized Fourier Transform,GFT)定義為[5-6]:

式(1)中,s0(t)是隨時間變化的實值函數。實際上是對x(t)e-j2πs0(t)作標準的傅里葉變換。同樣,可對XG(f)進行逆GFT變換,得到x(t):

如果 X G(f)=δ(f-f 0),則有 x(t)=ej2π[f0 t+s0(t)],也就是說一個瞬時頻率為 f(t)=f0+s′0(t)的信號經過合適的廣義傅里葉變換可以將能量集中于頻率 f=f0上,其中s′0(t)=d s0(t)/d t。因此只需找到一個近似于s0(t)的相位函數s(t),對信號x(t)進行廣義傅里葉變換,得到的解調函數x(t)e-j2πs(t)的時頻譜將是一條近似于平行時間坐標軸的,由f=f0確定的直線。這種信號方法稱為廣義解調[7-8]。

1.2 廣義S變換

Stockwell等學者于1996年首次提出了S變換(S-Transform)[9-10]:

式(4)中,w0(t-τ)為高斯窗口,τ為控制高斯窗口在時間t軸位置的參數,f為頻率。并且逆S變換可以表示為:

由于S變換中的基本小波形態固定,使得S變換在實際應用中受到一定限制。對S變換的高斯窗進行改造,便得到了廣義 S變換(Generalized S-transform)[11]:

式(6)中,λ、p為調節因子。λ＞0,選定 p后,當λ＞1時,則時窗寬度隨信號頻率呈反比變換的速度加快,λ＜1則減慢[12],使其具有更好的局部特性。當λ=p=1時,廣義S變換即簡化為基本的S變換,說明廣義S變換與S變換的計算相似,不會增加額外計算量。

2 基于廣義解調和廣義S變換的信號時頻濾波去噪

由前述可知,對于信號 x(t)=A ej2π[f0t+s0(t)],其廣義解調信號為 x(t)e-j2πs(t)=A ej2π[f0t+s0t]e-j2πs(t)≈A ej2πf0t,s(t)≈s0(t)為信號 x(t)相位函數的估計值。由式(1)和式(3)共同得到廣義解調信號的廣義S變換為:

式(7)中,W(j f)為改進的窗函數w(t)的傅里葉變換。式(7)表明,經廣義解調后的信號,其廣義S變換的時頻表示近似為一條平行于時間軸的由f=f0所確定的直線,其分辨率由窗函數本身所決定。

因此,首先對信號進行廣義解調和廣義S變換,將信號的有用部分在時頻域內近似表示為平行于時間軸的一條直線,然后通過時頻濾波濾除噪聲部分,結合式(7),時頻濾波器H(τ,f)的表達式為:式(8)中,f 0是由式(7)所確定的信號時頻表示的中心頻率,B為由時頻分辨率所決定的濾波器的帶寬。從式(8)可以看出,濾波器實際上為一平行于時間軸的二維矩形窗“帶通”濾波器。考慮到對信號突然截斷容易產生邊界效應等負面影響,因此對濾波器H(τ,f)進行改進,構造二維梯形窗濾器:

式(9)中,ω1為濾波器衰減起始頻率,ω2為截止頻率,可以根據信號具體的時頻分辨率選取合適的ω1,ω2。由此,可以重構降噪后的信號,結合式(5)和式(9),其重構過程為:

此時,y(t)仍為廣義解調信號,需對y(t)進行逆廣義解調,即:

則可得到經降噪后重構的原始信號y(t)。最后由y(t)時頻分布提取出相對準確的彈丸運動曲線。

時頻域去噪,關鍵是在時頻域設計出高效的二維濾波器。而對于瞬時頻率為非線性函數的信號,時頻濾波器需要兼顧時間和頻率兩個方面,往往比較復雜。由式(8)和式(9)可以看出,雖然需要在時頻域構造二維濾波器,但是通過廣義解調,此濾波器的參數可以完全和時間無關,有效地降低了濾波器的構造難度。并且,與小波變換和短時傅里葉變換相比,逆S變換實際上是廣義S變換的時頻分布對時間積分后,對其結果進行的傅里葉反變換。它所具有的能量無損特性,使其能夠在降噪后的計算過程中對數據精度不產生影響。此外,廣義S變換它所具有的高時頻分辨率以及其基本小波不必滿足容許性條件等獨特的優點,使得在克服常規時頻濾波缺陷的同時,可以盡可能完整地保留有用信息的時頻成分,對即使是信噪比較差的信號,也能得到較為精確的時頻曲線。

3 仿真及分析

3.1 試驗信號

高精度毫米波干涉運動參數測試系統,通常包括兩種微波源:低頻微波源和高頻微波源。低頻微波通常在1～12 GHz頻段內,而高頻微波通常選在35 GHz,55 GHz和95 GHz頻點上,此時火炮身管被視為自由場,彈丸被視為位置可變的反射體。由于采用空間自由場耦合方式,在波束范圍內的任何運動都會反映在干涉信號中[1]。因此干涉信號中的高頻干擾為噪聲的主要部分,所以對這部分噪聲的結果直接影響到最終的信噪比和精確度。

在此基礎上首先考察一測試信號,并對本文方法的有效性與基于RLS自適應濾波和EMD分解分析的方法進行對比。測試信號為:x(t)=sin(502.66t3-753.98t2+568.48 t+20.13)(加入2 dB高斯白噪聲以及20%的高頻干擾),信號采樣頻率為f s=1 000 Hz,采樣點數N=1 024,其時域波形如圖 1(a)所示 ,時頻譜如圖1(b)所示。由RLS自適應濾波和EMD分解的方法得到的時頻脊線圖噪聲抑制效果不理想,含有大量奇異點。在此基礎上提取的脊線擬合曲線并不十分準確。

采用本文的方法對試驗信號進行處理,首先估計信號的相位函數s(t)。從圖2(a)可以看到,當t=0 s時,f=85 Hz;t=0.5 s時,f=30 Hz;t=1 s時,f=85 Hz。假設瞬時頻率曲線為f(t)=s′(t)=at2+bt+c,根據曲線擬合可以確定相位函數s(t)=80 t2+120 t。這里并沒有利用已知信號參數獲得參數a、b,而是從圖中得到參數的近似值。隨后依照前述步驟進行處理。由圖2(b)看出信號的主要能量都集中在一條平行于時間軸的直線上,從而可以確定一個只與頻率有關的時頻域二維“帶通”濾波器。經過濾波和重構后的信號廣義S變換分布如圖2(d)所示,信號的噪聲幾乎被完全抑制。

為了說明此方法的通用性,在原始模擬信號的基礎上,再次加入瑞利分布噪聲,其廣義S分布如圖2(g)所示,去噪效果如2(h)所示。可以看出,對于其他類型的噪聲,此方法仍然能夠達到預期效果。

圖1 測試信號及基于RLS濾波和EMD分解的處理結果Fig.1 Test signal and the output through the method based on RLS filtering and EMD decomposition

圖2 廣義解調和廣義S變換的去噪結果Fig.2 Denoising output of generalized demodulation and generalized S-transform

圖3 對兩種方法得出的擬合曲線和理想瞬時頻率曲線進行了對比,在相同數量級下,文獻[4]方法所得曲線與理想曲線方差的有效數字為8.168 9;而本文方法的方差有效數字為1.785 8。可以看出:采用本文方法,信噪比在明顯提高的同時,精確度也得到大幅提高。在此基礎上提取的時頻脊線和曲線擬合表達準確,擬合程度高。

圖3 幾種擬合曲線和理想瞬時頻率曲線對比圖Fig.3 Comparison of fitting curves and the ideal instantaneous frequency curve

本文的方法對相位函數的選擇并不十分敏感。當相位函數估計不準確時,廣義解調后信號的時頻分布曲線將不會是平行于時間軸的直線[7,13]。然而,只要誤差不超過一定的范圍,保證在廣義解調后信號的時頻分布曲線包含在時頻空間的某一個矩形時頻帶中,由此相應的調節時頻濾波器參數,也能達到比較理想的處理結果。

3.2 實際數據分析

為了模擬內膛彈丸運動,采用自行研制的94 GHz干涉儀對某自制土火炮發射過程進行測試,采樣頻率2 MHz,獲得實測數據時域波形圖如圖4所示(限于篇幅限制,本文已截取有用信號段)。由于信號長度過長,將信號分兩組進行運算:第一組對應采樣時間0～0.12 s,第二組對應采樣時間0.12～0.21 s。如圖4所示,可以看出,采用本文所述方法進行處理,顯著提高了信噪比,并在此基礎上提取出幾乎不含奇異點的時頻脊線。

圖4 實測數據處理流程Fig.4 Test data processing

將時頻脊線進行組合,如圖5(a)所示。在此基礎上進行二次項曲線擬合,并根據多普勒頻率公式推導出彈丸在堂內運動的速度曲線及加速度曲線。由圖5(c)可知,彈丸出膛后速度在0.205 s左右達到最大,為28.16 m/s。相應彈丸運動距離,即反射板到炮口距離以及炮筒長度共2.310 2 m。與實際相符,驗證了結果的準確性。

圖5 彈丸運動曲線Fig.5 Curves of bullet motion

4 結論

本文提出一種基于廣義解調和廣義S變換的方法,使得廣義S變換的時頻分布為近似平行于時間軸的一條直線;然后在時頻域構造與時間因子無關的二維時頻帶通濾波器去除噪聲,降低了濾波器的構造難度;最后重構信號并由此提取彈丸運動曲線。仿真結果表明:與 RLS自適應濾波和希爾伯特黃(HHT)變換相比,本文方法顯著提高了信號的信噪比,增加了所得結果的準確度。此外,廣義S正反變換具有的能量無損特性,與奇異值分解(SVD)等降噪方法相比,不需要假設噪聲與信號的相關性差,并且計算量少,運算速度快。本文提出的方法能夠準確地表述彈丸運動模型,是一種有效的處理方法。但是,如何有效地選取相位函數以及設計更高效的時頻濾波器,還需要進一步的研究和完善。

[1]王黎明,趙英亮,韓焱.高速運動目標測試數據處理精度的方法研究[J].華北工學院學報,2005,26(1):53-57.WANG Liming,ZHAO Yingliang,HAN Yan.The methods of increasing the data processing precision of highspeed moving object[J].Journal of North China Institute of Technologe,2005,26(1):53-57.

[2]王黎明,趙昕,劉洪濤,等.微波干涉儀測試數據處理方法[J].火炮發射與控制學報,2004(4):63-66.WANG Liming,ZHAO Xin,LIU Hongtao,et al.Method of test data processing on microwave interferometer[J].Gun Launch&Control Journal,2004(4):63-66.

[3]肖劍,柳斌,郭亞龍,等.基于希爾伯特變換的干涉儀信號處理[J].探測與控制學報,2010,32(1):80-83.XIAO Jian,LIU Bin,GUO Yalong,et al.Signal processing of microwave interferometer based on hilbert transform[J].Journal of Detection&Control,2010,32(1):80-83.

[4]陶金泉.毫米波干涉儀信號處理研究[D].南京:南京理工大學電光學院,2009.

[5]Olhede S,Walden A T.A generalized demodulation approach to time-f requency projections formulticom-ponent signals[J].Proceeding s of the Royal Society A,2005,461(2059):2 159-2 179.

[6]Zhi Pengfang,Fulei Chu,Ming J Zuo.Time-frequency analysis of time-varying modulated signals based on improved energy separation by iterative generalized[J].Journal of Sound and Vibration,2011,33(6):1 225-1 243.

[7]楊宇,程軍圣,于德介.廣義解調時頻分析方法中的若干問題探討[J].振動與沖擊,2008,27(2):19-24.YANG Yu,CHENG Junsheng,YU Dejie.Study on some problems in the generalized demodulation time-f requency analysis method[J].Jouranl of Vibration and Shock,2008,27(2):19-24.

[8]羅向龍,高靜懷.基于廣義解調和奇異值分解的時頻表示增強[J].數據采集與處理,2010,25(4):419-424.LUO Xianglong,GAO Jinghuai.Enhancing time-frequency representation based on generalized demodulation and SVD[J].Journal of Data Acquisition&Processing,2010,25(4):419-424.

[9]Stockwell R G,et al.Localization of the complex spectrum:The S-transform[J].IEEE Transactions on Signal Processing,1996,44(4):998-1 001.

[10]Stockwell RG.A basis for efficient representation of the S-transform[J].Digital Signal Processing,2007,17(1):371-393.

[11]高靜懷,陳文超,李幼銘,等.廣義S變換與薄互層地震響應分析[J].地球物理學報,2003,46(4):526-532.GAO Jinghuai,CHEN Wenchao,LI Youming,et al.Generalized S transform and seismic responseanalysis of thin interbeds[J].Chinese Journal of Geophysics,2003,46(4):526-532.

[12]趙淑紅,朱光明.S變換時頻濾波去噪方法[J].石油地球物理勘探,2007,42(4):402-408.ZHAO Shuhong,ZHU Guangming.Time-frequency filtering to denoise by S transform[J].OGP,2007,42(4):402-408.

[13]Junsheng Cheng,Yu Yang,Dejie Yu.The envelope order spectrum based on generalized demodulation time-frequency analysis and its application to gear fault diagnosis[J].Mechanical Systems and Signal Processing,2010,24(2):508-521.