基于廣義解調(diào)和廣義S變換的時(shí)頻域去噪方法

柏 晨,肖澤龍,許建中

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

0 引言

利用電磁波對(duì)膛內(nèi)彈丸運(yùn)動(dòng)參數(shù)測(cè)試是目前比較主要的測(cè)量方式。然而,在利用該方法進(jìn)行測(cè)試時(shí),由于受到信號(hào)選擇、反射面確定、電離干擾及外界環(huán)境的影響,得到的干涉信號(hào)中存在大量的噪聲,信號(hào)信噪比差。因而,要得到高精度的測(cè)試結(jié)果,就對(duì)后續(xù)的數(shù)據(jù)處理提出了較高的要求[1]。

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

由于傳統(tǒng)的時(shí)域或頻域去噪方法經(jīng)常受到使用條件的限制,本文提出了基于廣義解調(diào)和廣義S變換的時(shí)頻去噪方法。

1 理論基礎(chǔ)

1.1 廣義解調(diào)

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

式(1)中,s0(t)是隨時(shí)間變化的實(shí)值函數(shù)。實(shí)際上是對(duì)x(t)e-j2πs0(t)作標(biāo)準(zhǔn)的傅里葉變換。同樣,可對(duì)XG(f)進(jìn)行逆GFT變換,得到x(t):

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

1.2 廣義S變換

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

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

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

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

2 基于廣義解調(diào)和廣義S變換的信號(hào)時(shí)頻濾波去噪

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

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

因此,首先對(duì)信號(hào)進(jìn)行廣義解調(diào)和廣義S變換,將信號(hào)的有用部分在時(shí)頻域內(nèi)近似表示為平行于時(shí)間軸的一條直線,然后通過(guò)時(shí)頻濾波濾除噪聲部分,結(jié)合式(7),時(shí)頻濾波器H(τ,f)的表達(dá)式為:式(8)中,f 0是由式(7)所確定的信號(hào)時(shí)頻表示的中心頻率,B為由時(shí)頻分辨率所決定的濾波器的帶寬。從式(8)可以看出,濾波器實(shí)際上為一平行于時(shí)間軸的二維矩形窗“帶通”濾波器?？紤]到對(duì)信號(hào)突然截?cái)嗳菀桩a(chǎn)生邊界效應(yīng)等負(fù)面影響,因此對(duì)濾波器H(τ,f)進(jìn)行改進(jìn),構(gòu)造二維梯形窗濾器:

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

此時(shí),y(t)仍為廣義解調(diào)信號(hào),需對(duì)y(t)進(jìn)行逆廣義解調(diào),即:

則可得到經(jīng)降噪后重構(gòu)的原始信號(hào)y(t)。最后由y(t)時(shí)頻分布提取出相對(duì)準(zhǔn)確的彈丸運(yùn)動(dòng)曲線。

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

3 仿真及分析

3.1 試驗(yàn)信號(hào)

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

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

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

為了說(shuō)明此方法的通用性,在原始模擬信號(hào)的基礎(chǔ)上,再次加入瑞利分布噪聲,其廣義S分布如圖2(g)所示,去噪效果如2(h)所示?？梢钥闯?對(duì)于其他類型的噪聲,此方法仍然能夠達(dá)到預(yù)期效果。

圖1 測(cè)試信號(hào)及基于RLS濾波和EMD分解的處理結(jié)果Fig.1 Test signal and the output through the method based on RLS filtering and EMD decomposition

圖2 廣義解調(diào)和廣義S變換的去噪結(jié)果Fig.2 Denoising output of generalized demodulation and generalized S-transform

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

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

本文的方法對(duì)相位函數(shù)的選擇并不十分敏感。當(dāng)相位函數(shù)估計(jì)不準(zhǔn)確時(shí),廣義解調(diào)后信號(hào)的時(shí)頻分布曲線將不會(huì)是平行于時(shí)間軸的直線[7,13]。然而,只要誤差不超過(guò)一定的范圍,保證在廣義解調(diào)后信號(hào)的時(shí)頻分布曲線包含在時(shí)頻空間的某一個(gè)矩形時(shí)頻帶中,由此相應(yīng)的調(diào)節(jié)時(shí)頻濾波器參數(shù),也能達(dá)到比較理想的處理結(jié)果。

3.2 實(shí)際數(shù)據(jù)分析

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

圖4 實(shí)測(cè)數(shù)據(jù)處理流程Fig.4 Test data processing

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

圖5 彈丸運(yùn)動(dòng)曲線Fig.5 Curves of bullet motion

4 結(jié)論

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

[1]王黎明,趙英亮,韓焱.高速運(yùn)動(dòng)目標(biāo)測(cè)試數(shù)據(jù)處理精度的方法研究[J].華北工學(xué)院學(xué)報(bào),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]王黎明,趙昕,劉洪濤,等.微波干涉儀測(cè)試數(shù)據(jù)處理方法[J].火炮發(fā)射與控制學(xué)報(bào),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]肖劍,柳斌,郭亞龍,等.基于希爾伯特變換的干涉儀信號(hào)處理[J].探測(cè)與控制學(xué)報(bào),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]陶金泉.毫米波干涉儀信號(hào)處理研究[D].南京:南京理工大學(xué)電光學(xué)院,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]楊宇,程軍圣,于德介.廣義解調(diào)時(shí)頻分析方法中的若干問(wèn)題探討[J].振動(dòng)與沖擊,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]羅向龍,高靜懷.基于廣義解調(diào)和奇異值分解的時(shí)頻表示增強(qiáng)[J].數(shù)據(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變換與薄互層地震響應(yīng)分析[J].地球物理學(xué)報(bào),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變換時(shí)頻濾波去噪方法[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.