王宏偉 趙國慶
(西安電子科技大學 電子對抗研究所 陜西 西安710071)
傳統(tǒng)的時頻分析方法包括短時傅立葉變換(STFT)[1],維格納變換(WVD)[2],小波變換(WT)[3],滑動 DFT[4-6],滑窗 FFT[7-9]等 。STFT 受測不準原理制約,時、頻分辨率較差;WVD存在交叉項干擾;WT計算量大;單點滑動DFT運算速度快,但數(shù)據(jù)不能在時域加窗以減小頻譜泄露[6];滑窗FFT由于使用了技術(shù)成熟的FFT器件,運算速度較DFT快,但只能毫無選擇地計算出全部離散頻率點的頻譜。
遞歸算法可以有選擇地計算局部頻段的連續(xù)頻譜,實時地進行時頻分析,時域參數(shù)測量和頻域參數(shù)測量。通過工作參數(shù)的設(shè)置與調(diào)整,表明遞歸算法不僅時間分辨率和頻率分辨率可以調(diào)整,而且所分析的頻率或頻段可靈活設(shè)置。該算法具有運算速度快、數(shù)據(jù)存儲量少、資源占用量與工作參數(shù)的變化無關(guān)等特點。
圖1為遞歸算法的基本遞歸單元(單一通路)結(jié)構(gòu)圖,時域數(shù)據(jù)流經(jīng)數(shù)據(jù)整理,獲取以n時刻為起點,長度為 N的一段時域數(shù)據(jù)幀,加窗濾波減少頻譜泄露(一般選擇漢寧窗,窗口寬度為 N),進入遞歸運算,N次循環(huán)迭代后,乘以權(quán)系數(shù)e便得到了n時刻起長度為N的時域數(shù)據(jù)在頻率點f處的離散傅立葉變換系數(shù)y(n,f),其中 fs為采樣頻率。

圖1 基本遞歸單元的結(jié)構(gòu)圖
遞歸算法滿足循環(huán)迭代公式

式中:r為迭代次數(shù);y(n,f,r)為第r次迭代結(jié)果。令初始狀態(tài)y(n,f,0)=0,則


那么

在表達式y(tǒng)(n,f,N)中:n為時域數(shù)據(jù)幀的起始時刻;N為循環(huán)迭代次數(shù),其值等于參與循環(huán)迭代運算的數(shù)據(jù)幀長度和窗口寬度;f為所分析頻率,在頻率軸上可以選擇任意實數(shù)值。當·k時,,此處 k=0,1,…,N-1,式(3)可改寫為

式(4)就是傳統(tǒng)的滑窗FFT運算,它得到了離散頻率點k處的傅立葉變換系數(shù)。表明式(4)只是式(3)的特例,式(3)得到的是任意實數(shù)(包括 f s/N整數(shù)倍)頻率點的離散時間(DTFT)傅立葉系數(shù),而式(4)像滑窗FFT算法一樣只能得到 f s/N整數(shù)倍頻率點的離散傅立葉(DFT)系數(shù)。
遞歸算法是建立在有限長數(shù)據(jù)的離散傅立葉變換基礎(chǔ)上的,因此,具有和離散時間傅立葉變換相似的性質(zhì),包括下列性質(zhì):
1)線性性質(zhì)。遞歸算法屬于線性變換,不會產(chǎn)生新的頻率分量,不會受交叉項干擾。
從濾波器角度看,這種帶反饋結(jié)構(gòu)的基本遞歸單元可以被看成用一階 IIR濾波器實現(xiàn)的 N階FIR濾波器。
雖然通過對工作參數(shù)n,f,N的設(shè)置和調(diào)整,整個二維時頻譜的計算完全可以由單一通路完成,但是為了提高時頻分析效率,需要將多個通路集成為不同結(jié)構(gòu)的IIR濾波器組,如圖2,文中給出了兩種結(jié)構(gòu)。濾波器組內(nèi)部各通路相互獨立,并行工作,通路數(shù)目M根據(jù)需要可獨立增減和使用。
如果將IIR濾波器組內(nèi)各通路的參數(shù) N和n設(shè)置為相同規(guī)律變化,而各通路的參數(shù) f m,m=1,…,M不相同,結(jié)構(gòu)如圖2(a)。設(shè)計圖2(a)結(jié)構(gòu)的濾波器組時,可將數(shù)據(jù)整理工作移到濾波器組外進行。
根據(jù)實際需要,利用圖2(a)結(jié)構(gòu),可以有選擇地計算任何時刻起點n的一幀N長時域數(shù)據(jù)對應(yīng)的任意局部頻段的譜線,比只能計算全部離散頻譜的FFT算法更加經(jīng)濟有效,更加靈活。
利用各通路的獨立性,如果將所有通路都設(shè)置為對某特定頻率點 f 0(一般選信號中心頻率或載頻)的遞歸運算,N仍采用相同規(guī)律變化,而各個通路所利用的數(shù)據(jù)幀起點nm不同,m=1,…,M,結(jié)構(gòu)如圖2(b)。各通路初始的數(shù)據(jù)幀起始時刻 的設(shè)置可參考

循環(huán)流水作業(yè)時,各通路數(shù)據(jù)幀的起始時刻nm,j可設(shè)置為


式中:M 為通路數(shù)目;j=1,…,∞;nm,j表示第m路第j級流水作業(yè)時數(shù)據(jù)幀的起始時刻;n0為研究開始時刻點;d為相鄰?fù)窋?shù)據(jù)幀起始時刻的時間間隔,即相鄰?fù)窋?shù)據(jù)幀間滑移量,d∈1,…,N.時域數(shù)據(jù)流經(jīng)圖2(b)結(jié)構(gòu)處理后得到某特定頻率f0信號不同時刻點的頻域信息,滿足時間分辨率的前提下,用門限檢測的方法就可獲取信號的起始時刻,終止時刻,持續(xù)時間等時域信息。
由于圖2(b)結(jié)構(gòu)的各通路以時空轉(zhuǎn)化方式,將長時間跨度的時域數(shù)據(jù)轉(zhuǎn)化到并行的通路空間上,并行循環(huán)流水作業(yè)或依次循環(huán)流水作業(yè),并且相鄰?fù)窋?shù)據(jù)幀間滑移量d可以調(diào)整,因此,運算的實時性、時域參數(shù)測量精度與整體運算量之間的矛盾問題容易協(xié)調(diào)解決。
在圖2中,每個通路的計算量僅為N+1次復(fù)數(shù)乘(加)。由于各通路的并行工作,因此,多路遞歸算法的計算量也為 N+1次復(fù)數(shù)乘(加),相比于滑窗FFT算法[9]的 N log2N次復(fù)數(shù)乘(加),計算量小。
通路數(shù)目M較多時,從宏觀上看,濾波器組規(guī)模龐大,但每一路的遞歸運算只需要一個加法器、兩個乘法器和一個鎖存器,通路內(nèi)部器件簡單,各通路結(jié)構(gòu)相同,易于多個通路集成,數(shù)據(jù)存儲量少,資源占用量與工作參數(shù)n,f,N的變化無關(guān)。
通過對工作參數(shù)n,f,N的設(shè)置與調(diào)整,結(jié)合具體實例來詳細討論遞歸算法的其它特點。
設(shè)有同時到達的兩個單載頻脈沖信號,一個信號的載頻 f 1=30 MHz,占空比40%;另一個信號的載頻 f 2=30.8 MHz,占空比60%;脈沖重復(fù)頻率均為 fp=100 k Hz,采樣頻率為 fs=100 MHz;信噪比SNR=15 d B;圖2的兩種IIR濾波器組內(nèi)單一通路的集成數(shù)目M均為100。
信號檢測的目的是檢測在有效全局頻段0~50 MHz內(nèi)有無信號,存在信號的數(shù)目及其所在的局部頻段。若指定相鄰譜線間隔Δfgrid=0.1 MHz,那么實際需要的通路數(shù)目M=50/0.1=500>100路。解決方案1:并行方式。將五塊圖2(a)結(jié)構(gòu)的IIR濾波器組并行使用,各通路頻率分別設(shè)置為 f m=0.1×m(MHz),通路編號m=0,1,…,499。解決方案2:時分方式。僅用一塊圖2(a)結(jié)構(gòu)的IIR濾波器組,對于輸入的每一幀數(shù)據(jù),數(shù)據(jù)不變(d=0)僅改變各通路被測頻率,運算五次。其中第l次將各通路頻率分別設(shè)置為 f m=0.1·(l-1)+m/2(MHz),l=1,…,5,m=0,1,…,99。當其它參數(shù):n=0,N分別取100和350時,得到全頻段信號檢測的幅頻譜和相頻圖如圖3(a)。
由圖3(a)可知,在30 MHz左右的局部頻段存在信號,并且只有當循環(huán)迭代次數(shù)N比較大時才能可靠地發(fā)現(xiàn)此局部頻段存在兩個信號。
為了進一步測量信號的頻域參數(shù),調(diào)整各通路頻率,對有信號(或感興趣)的局部頻段,進行頻譜細化分析。局部頻段分析時,各通路頻率設(shè)置可參考下式

通過對被測頻段中心的調(diào)整并配合頻譜細化技術(shù),可以使相鄰譜線間隔Δf grid=α·f s/N(調(diào)整頻譜細化系數(shù)α)變得很小(滿足頻率分辨率的前提下,譜線間隔Δfgrid越小,頻域參數(shù)測量精度越高),得到一幀可移動中心位置(調(diào)整)的可變局部頻段的近似連續(xù)頻譜,從而可以分析任意感興趣的局部頻段的頻譜細節(jié)。若不考慮噪聲、頻譜混疊等因素的影響,利用圖3(b),理論上頻域參數(shù)的測量誤差可以做到趨于零。

雖然利用遞歸算法可以完成全頻段信號檢測任務(wù),但需要消耗大量資源或以犧牲部分實時性為代價。而在局部頻段分析時,無須增加資源或同幀數(shù)據(jù)多次利用便可保證信號實時檢測、頻譜細化分析和頻域參數(shù)精確測量。可見,遞歸算法比較適合局部頻段的時頻分析。在實際應(yīng)用中,也可配合其它算法先進行頻段粗引導(dǎo),再用遞歸算法進行局部頻段精細分析。
前面用一幀數(shù)據(jù)已獲得某時刻起N長時域數(shù)據(jù)對應(yīng)的頻譜分布,隨著參數(shù)n的連續(xù)變化(由數(shù)據(jù)幀的滑移量d體現(xiàn),d∈{1,2,…,N}),計算不同起始時刻數(shù)據(jù)幀(N保持不變)對應(yīng)的頻譜分布,聯(lián)合起來便可以得到二維時頻譜圖。
仍以2.1節(jié)中的兩脈沖信號為例。當工作參數(shù):f m=(30.4+m·0.15·f s/N)(MHz),m=-50,-49,…,49;N=350;d=10時得到的二維時頻圖,如圖4。
在二維時頻圖中,時頻分析系統(tǒng)可以根據(jù)各路輸出結(jié)果,調(diào)整工作參數(shù),在顯示器上觀察不同參數(shù)條件下,信號實時動態(tài)變化的特點。系統(tǒng)也可以在脫機條件下,針對信號特點,尋找最佳工作參數(shù)并精確測量時頻參數(shù)。
圖2(a)結(jié)構(gòu)雖然采用并行工作模式,每幀數(shù)據(jù)對單一通路而言,計算量不大,但滑移量d取得較小時(滿足時間分辨率的前提下,d越小,時頻譜網(wǎng)格的時域間隔Δtgrid=d/fs越小,則時域參數(shù)測量精度就越高),總體運算量變得很大,難以保證二維時頻分析的實時性,尤其是時域參數(shù)實時測量。利用圖2(b)結(jié)構(gòu)來解決單載頻信號時域參數(shù)實時測量問題。
當工作參數(shù)f0=30.8 MHz;N=350;而各通路數(shù)據(jù)幀起始時刻按照式(5)和式(6)設(shè)置,其中d=1(數(shù)據(jù)幀逐點滑移),n0=0時,得到載頻 f 0=30.8 MHz的脈沖信號時域參數(shù)測量圖,如圖5(a),圖5(b)為該脈沖信號理論占空比示意圖。
當窗寬N小于脈內(nèi)(或時限信號持續(xù)時間內(nèi))有效采樣點數(shù)時,圖5(a)形狀為等腰梯形的脈沖串,通常取脈頂頻譜幅度平均值的一半作為檢測門限,統(tǒng)計某單個等腰梯形脈沖首次高于和首次低于檢測門限的時刻n1,n2。由于遞歸算法公式中采用左對齊方式加窗,因此,由式(8)計算該脈沖的起始時刻n r,終止時刻n f和脈寬 n w。

因為信號在頻域具有能量集中的特性[9],遞歸算法采用頻域數(shù)據(jù)來測量時域參數(shù),相比用時域數(shù)據(jù)測量時域參數(shù),抗噪性能加強了,提高了檢測靈敏度。

雖然通過對頻譜細化系數(shù)α和數(shù)據(jù)幀的滑移量d的調(diào)整,可以使時頻網(wǎng)格的頻域間隔Δf grid和時域間隔Δtgrid達到很小,但是并不意味著頻率分辨率Δf(指能夠區(qū)分頻率軸上靠得很近的兩信號或兩頻率分量的最小頻率間隔)或時域分辨率Δt(指能夠區(qū)分時間軸上靠得很近的兩信號或時限信號的起始時刻、終止時刻的最小時間間隔)的改善。
遞歸運算的頻率分辨率Δf和時間分辨率 Δt受采樣頻率f s、窗口寬度N(等于循環(huán)迭代次數(shù))、窗函數(shù)形狀等因素的影響[7]。Δf正比于f s/N,Δt反比于fs/N。N變大,頻域分辨率提高,而時域分辨率下降,二者仍然受制于“測不準原理”。但在一定的范圍內(nèi),對循環(huán)迭代次數(shù) N的調(diào)整,使得遞歸運算具有了可調(diào)可控的時間分辨率、頻率分辨率。
假設(shè)有一個分段正弦信號

其中:f 1=52 Hz,f 2=54 Hz,f 3=56 Hz,f 4=58 Hz,f s=200 Hz。當工作參數(shù) f m=(50+0.01m)(Hz),m=0,1,…,99;d=10;N 分別取 100、200和400時,利用遞歸算法得到的結(jié)果,如圖6。通過對比可以看到,N較小時,具有較好的時間分辨率,而相應(yīng)的頻率分辨率則不高,隨著N的增加,頻率分辨率變得越高,但此時的時間分辨率則會相應(yīng)地下降。
另外,由圖6還可以看到,參數(shù) N除了影響時、頻分辨率外,還影響著譜線的幅度。N越大,譜線的幅度越大,則系統(tǒng)的抗噪性能越好。
遞歸算法具有很大的靈活性,非常適合對局部頻段實時地進行頻域參數(shù)測量、時域參數(shù)測量和二維時頻譜分析,三項任務(wù)可獨立完成,也可相互配合完成。時頻分析系統(tǒng)需要根據(jù)某種規(guī)則做出工作參數(shù)的設(shè)置與調(diào)整,使其在保證實時性,時間分辨率,頻率分辨率,時、頻參數(shù)測量精度,抗噪性能等性能指標時統(tǒng)籌調(diào)度通路資源,協(xié)調(diào)工作。
[1] 李岳霖,龐偉正,等.多普勒雷達測量輕武器彈丸轉(zhuǎn)速方法研究[J].電波科學學報,2007,22(3):502-507.
LI Yuelin,PANG Weizheng,et al.Measuring method of rotation velocity for small arms projectiles by Doppler radar[J].Chinese Journal of Radio Science,2007,22(3):502-507.(in Chinese)
[2] 袁偉明,王 敏,吳順君.一種新的 LPI信號的截獲方法[J].電波科學學報,2005,20(1):73-76.
YUAN Weiming,WANG Min,WU Shunjun.A novel interception of LPI signals[J].Chinese Journal of Radio Science,2005,20(1):73-76.(in Chinese)
[3] 徐玉清,郎 銳,高 攀.基于小波的信息保持型雷場圖像融合算法[J].電波科學學報,2009,24(2):233-237.
XU Yuqing,LANG Rui,GAO Pan.Wavelet-based information restorefusion method of landmine field images[J].Chinese Journal of Radio Science,2009,24(2):233-237.(in Chinese)
[4] 王小龍.DFT和交替DFT調(diào)制濾波器組設(shè)計算法研究[D].陜西西安:西安電子科技大學,2008.
[5] JACOBSEN E,LYONS R.The sliding DFT[J].IEEE Signal Processing Magazine,2003,20(3):74-80.
[6] 黃寒華.滑動DFT算法研究[D].江蘇南京:東南大學,2006.
[7] DRESSLER K.Sinusoidal extraction using an efficient implementation of a multi-resolution FFT[C]∥Proc.of the 9th Int.Conference on Digital Audio Effects,Montreal,Canada,September 18-20,2006.
[8] FARHANG-BOROUJENY B,GAZORS.Generalized sliding FFT and its application to implementation of block LMS adaptive filters[J].IEEE Transactions on Signal Processing,1994,42(3):532-538.
[9] 劉 平,靳成英,陳曾平.一種基于短時FFT的寬帶數(shù)字接收機設(shè)計[J].信號處理,2008,24(6):988-991.
LIU Ping,JIN Chengying,CHEN Zengping.A shorttime FFT based design for wideband digital reconnaissance receiver[J].Signal Processing,2008,24(6):988-991.(in Chinese)