王宏偉 趙國(guó)慶

(西安電子科技大學(xué) 電子對(duì)抗研究所 陜西 西安710071)

1.引 言

傳統(tǒng)的時(shí)頻分析方法包括短時(shí)傅立葉變換(STFT)[1],維格納變換(WVD)[2],小波變換(WT)[3],滑動(dòng) DFT[4-6],滑窗 FFT[7-9]等 。STFT 受測(cè)不準(zhǔn)原理制約,時(shí)、頻分辨率較差;WVD存在交叉項(xiàng)干擾;WT計(jì)算量大;單點(diǎn)滑動(dòng)DFT運(yùn)算速度快,但數(shù)據(jù)不能在時(shí)域加窗以減小頻譜泄露[6];滑窗FFT由于使用了技術(shù)成熟的FFT器件,運(yùn)算速度較DFT快,但只能毫無(wú)選擇地計(jì)算出全部離散頻率點(diǎn)的頻譜。

遞歸算法可以有選擇地計(jì)算局部頻段的連續(xù)頻譜,實(shí)時(shí)地進(jìn)行時(shí)頻分析,時(shí)域參數(shù)測(cè)量和頻域參數(shù)測(cè)量。通過(guò)工作參數(shù)的設(shè)置與調(diào)整,表明遞歸算法不僅時(shí)間分辨率和頻率分辨率可以調(diào)整,而且所分析的頻率或頻段可靈活設(shè)置。該算法具有運(yùn)算速度快、數(shù)據(jù)存儲(chǔ)量少、資源占用量與工作參數(shù)的變化無(wú)關(guān)等特點(diǎn)。

2.遞歸算法的介紹

2.1 基本遞歸單元的結(jié)構(gòu)

圖1為遞歸算法的基本遞歸單元(單一通路)結(jié)構(gòu)圖,時(shí)域數(shù)據(jù)流經(jīng)數(shù)據(jù)整理,獲取以n時(shí)刻為起點(diǎn),長(zhǎng)度為 N的一段時(shí)域數(shù)據(jù)幀,加窗濾波減少頻譜泄露(一般選擇漢寧窗,窗口寬度為 N),進(jìn)入遞歸運(yùn)算,N次循環(huán)迭代后,乘以權(quán)系數(shù)e便得到了n時(shí)刻起長(zhǎng)度為N的時(shí)域數(shù)據(jù)在頻率點(diǎn)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,則

那么

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

式(4)就是傳統(tǒng)的滑窗FFT運(yùn)算,它得到了離散頻率點(diǎn)k處的傅立葉變換系數(shù)。表明式(4)只是式(3)的特例,式(3)得到的是任意實(shí)數(shù)(包括 f s/N整數(shù)倍)頻率點(diǎn)的離散時(shí)間(DTFT)傅立葉系數(shù),而式(4)像滑窗FFT算法一樣只能得到 f s/N整數(shù)倍頻率點(diǎn)的離散傅立葉(DFT)系數(shù)。

遞歸算法是建立在有限長(zhǎng)數(shù)據(jù)的離散傅立葉變換基礎(chǔ)上的,因此,具有和離散時(shí)間傅立葉變換相似的性質(zhì),包括下列性質(zhì):

1)線性性質(zhì)。遞歸算法屬于線性變換,不會(huì)產(chǎn)生新的頻率分量,不會(huì)受交叉項(xiàng)干擾。

從濾波器角度看,這種帶反饋結(jié)構(gòu)的基本遞歸單元可以被看成用一階 IIR濾波器實(shí)現(xiàn)的 N階FIR濾波器。

2.2 多路遞歸算法的結(jié)構(gòu)

雖然通過(guò)對(duì)工作參數(shù)n,f,N的設(shè)置和調(diào)整,整個(gè)二維時(shí)頻譜的計(jì)算完全可以由單一通路完成,但是為了提高時(shí)頻分析效率,需要將多個(gè)通路集成為不同結(jié)構(gòu)的IIR濾波器組,如圖2,文中給出了兩種結(jié)構(gòu)。濾波器組內(nèi)部各通路相互獨(dú)立,并行工作,通路數(shù)目M根據(jù)需要可獨(dú)立增減和使用。

如果將IIR濾波器組內(nèi)各通路的參數(shù) N和n設(shè)置為相同規(guī)律變化,而各通路的參數(shù) f m,m=1,…,M不相同,結(jié)構(gòu)如圖2(a)。設(shè)計(jì)圖2(a)結(jié)構(gòu)的濾波器組時(shí),可將數(shù)據(jù)整理工作移到濾波器組外進(jìn)行。

根據(jù)實(shí)際需要,利用圖2(a)結(jié)構(gòu),可以有選擇地計(jì)算任何時(shí)刻起點(diǎn)n的一幀N長(zhǎng)時(shí)域數(shù)據(jù)對(duì)應(yīng)的任意局部頻段的譜線,比只能計(jì)算全部離散頻譜的FFT算法更加經(jīng)濟(jì)有效,更加靈活。

利用各通路的獨(dú)立性,如果將所有通路都設(shè)置為對(duì)某特定頻率點(diǎn) f 0(一般選信號(hào)中心頻率或載頻)的遞歸運(yùn)算,N仍采用相同規(guī)律變化,而各個(gè)通路所利用的數(shù)據(jù)幀起點(diǎn)nm不同,m=1,…,M,結(jié)構(gòu)如圖2(b)。各通路初始的數(shù)據(jù)幀起始時(shí)刻 的設(shè)置可參考

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

式中:M 為通路數(shù)目;j=1,…,∞;nm,j表示第m路第j級(jí)流水作業(yè)時(shí)數(shù)據(jù)幀的起始時(shí)刻;n0為研究開(kāi)始時(shí)刻點(diǎn);d為相鄰?fù)窋?shù)據(jù)幀起始時(shí)刻的時(shí)間間隔,即相鄰?fù)窋?shù)據(jù)幀間滑移量,d∈1,…,N.時(shí)域數(shù)據(jù)流經(jīng)圖2(b)結(jié)構(gòu)處理后得到某特定頻率f0信號(hào)不同時(shí)刻點(diǎn)的頻域信息,滿足時(shí)間分辨率的前提下,用門限檢測(cè)的方法就可獲取信號(hào)的起始時(shí)刻,終止時(shí)刻,持續(xù)時(shí)間等時(shí)域信息。

由于圖2(b)結(jié)構(gòu)的各通路以時(shí)空轉(zhuǎn)化方式,將長(zhǎng)時(shí)間跨度的時(shí)域數(shù)據(jù)轉(zhuǎn)化到并行的通路空間上,并行循環(huán)流水作業(yè)或依次循環(huán)流水作業(yè),并且相鄰?fù)窋?shù)據(jù)幀間滑移量d可以調(diào)整,因此,運(yùn)算的實(shí)時(shí)性、時(shí)域參數(shù)測(cè)量精度與整體運(yùn)算量之間的矛盾問(wèn)題容易協(xié)調(diào)解決。

2.3 多路遞歸算法的特點(diǎn)

在圖2中,每個(gè)通路的計(jì)算量?jī)H為N+1次復(fù)數(shù)乘(加)。由于各通路的并行工作,因此,多路遞歸算法的計(jì)算量也為 N+1次復(fù)數(shù)乘(加),相比于滑窗FFT算法[9]的 N log2N次復(fù)數(shù)乘(加),計(jì)算量小。

通路數(shù)目M較多時(shí),從宏觀上看,濾波器組規(guī)模龐大,但每一路的遞歸運(yùn)算只需要一個(gè)加法器、兩個(gè)乘法器和一個(gè)鎖存器,通路內(nèi)部器件簡(jiǎn)單,各通路結(jié)構(gòu)相同,易于多個(gè)通路集成,數(shù)據(jù)存儲(chǔ)量少,資源占用量與工作參數(shù)n,f,N的變化無(wú)關(guān)。

通過(guò)對(duì)工作參數(shù)n,f,N的設(shè)置與調(diào)整,結(jié)合具體實(shí)例來(lái)詳細(xì)討論遞歸算法的其它特點(diǎn)。

3.參數(shù)設(shè)置與調(diào)整

3.1 參數(shù) f的設(shè)置與調(diào)整——實(shí)現(xiàn)信號(hào)檢測(cè)、頻譜細(xì)化分析、頻域參數(shù)測(cè)量

設(shè)有同時(shí)到達(dá)的兩個(gè)單載頻脈沖信號(hào),一個(gè)信號(hào)的載頻 f 1=30 MHz,占空比40%;另一個(gè)信號(hào)的載頻 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。

信號(hào)檢測(cè)的目的是檢測(cè)在有效全局頻段0～50 MHz內(nèi)有無(wú)信號(hào),存在信號(hào)的數(shù)目及其所在的局部頻段。若指定相鄰譜線間隔Δfgrid=0.1 MHz,那么實(shí)際需要的通路數(shù)目M=50/0.1=500＞100路。解決方案1:并行方式。將五塊圖2(a)結(jié)構(gòu)的IIR濾波器組并行使用,各通路頻率分別設(shè)置為 f m=0.1×m(MHz),通路編號(hào)m=0,1,…,499。解決方案2:時(shí)分方式。僅用一塊圖2(a)結(jié)構(gòu)的IIR濾波器組,對(duì)于輸入的每一幀數(shù)據(jù),數(shù)據(jù)不變(d=0)僅改變各通路被測(cè)頻率,運(yùn)算五次。其中第l次將各通路頻率分別設(shè)置為 f m=0.1·(l-1)+m/2(MHz),l=1,…,5,m=0,1,…,99。當(dāng)其它參數(shù):n=0,N分別取100和350時(shí),得到全頻段信號(hào)檢測(cè)的幅頻譜和相頻圖如圖3(a)。

由圖3(a)可知,在30 MHz左右的局部頻段存在信號(hào),并且只有當(dāng)循環(huán)迭代次數(shù)N比較大時(shí)才能可靠地發(fā)現(xiàn)此局部頻段存在兩個(gè)信號(hào)。

為了進(jìn)一步測(cè)量信號(hào)的頻域參數(shù),調(diào)整各通路頻率,對(duì)有信號(hào)(或感興趣)的局部頻段,進(jìn)行頻譜細(xì)化分析。局部頻段分析時(shí),各通路頻率設(shè)置可參考下式

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

雖然利用遞歸算法可以完成全頻段信號(hào)檢測(cè)任務(wù),但需要消耗大量資源或以犧牲部分實(shí)時(shí)性為代價(jià)。而在局部頻段分析時(shí),無(wú)須增加資源或同幀數(shù)據(jù)多次利用便可保證信號(hào)實(shí)時(shí)檢測(cè)、頻譜細(xì)化分析和頻域參數(shù)精確測(cè)量?？梢?jiàn),遞歸算法比較適合局部頻段的時(shí)頻分析。在實(shí)際應(yīng)用中,也可配合其它算法先進(jìn)行頻段粗引導(dǎo),再用遞歸算法進(jìn)行局部頻段精細(xì)分析。

3.2 參數(shù)n的設(shè)置與調(diào)整——實(shí)現(xiàn)二維時(shí)頻譜分析和時(shí)域參數(shù)的實(shí)時(shí)測(cè)量

前面用一幀數(shù)據(jù)已獲得某時(shí)刻起N長(zhǎng)時(shí)域數(shù)據(jù)對(duì)應(yīng)的頻譜分布,隨著參數(shù)n的連續(xù)變化(由數(shù)據(jù)幀的滑移量d體現(xiàn),d∈{1,2,…,N}),計(jì)算不同起始時(shí)刻數(shù)據(jù)幀(N保持不變)對(duì)應(yīng)的頻譜分布,聯(lián)合起來(lái)便可以得到二維時(shí)頻譜圖。

仍以2.1節(jié)中的兩脈沖信號(hào)為例。當(dāng)工作參數(shù):f m=(30.4+m·0.15·f s/N)(MHz),m=-50,-49,…,49;N=350;d=10時(shí)得到的二維時(shí)頻圖,如圖4。

在二維時(shí)頻圖中,時(shí)頻分析系統(tǒng)可以根據(jù)各路輸出結(jié)果,調(diào)整工作參數(shù),在顯示器上觀察不同參數(shù)條件下,信號(hào)實(shí)時(shí)動(dòng)態(tài)變化的特點(diǎn)。系統(tǒng)也可以在脫機(jī)條件下,針對(duì)信號(hào)特點(diǎn),尋找最佳工作參數(shù)并精確測(cè)量時(shí)頻參數(shù)。

圖2(a)結(jié)構(gòu)雖然采用并行工作模式,每幀數(shù)據(jù)對(duì)單一通路而言,計(jì)算量不大,但滑移量d取得較小時(shí)(滿足時(shí)間分辨率的前提下,d越小,時(shí)頻譜網(wǎng)格的時(shí)域間隔Δtgrid=d/fs越小,則時(shí)域參數(shù)測(cè)量精度就越高),總體運(yùn)算量變得很大,難以保證二維時(shí)頻分析的實(shí)時(shí)性,尤其是時(shí)域參數(shù)實(shí)時(shí)測(cè)量。利用圖2(b)結(jié)構(gòu)來(lái)解決單載頻信號(hào)時(shí)域參數(shù)實(shí)時(shí)測(cè)量問(wèn)題。

當(dāng)工作參數(shù)f0=30.8 MHz;N=350;而各通路數(shù)據(jù)幀起始時(shí)刻按照式(5)和式(6)設(shè)置,其中d=1(數(shù)據(jù)幀逐點(diǎn)滑移),n0=0時(shí),得到載頻 f 0=30.8 MHz的脈沖信號(hào)時(shí)域參數(shù)測(cè)量圖,如圖5(a),圖5(b)為該脈沖信號(hào)理論占空比示意圖。

當(dāng)窗寬N小于脈內(nèi)(或時(shí)限信號(hào)持續(xù)時(shí)間內(nèi))有效采樣點(diǎn)數(shù)時(shí),圖5(a)形狀為等腰梯形的脈沖串,通常取脈頂頻譜幅度平均值的一半作為檢測(cè)門限,統(tǒng)計(jì)某單個(gè)等腰梯形脈沖首次高于和首次低于檢測(cè)門限的時(shí)刻n1,n2。由于遞歸算法公式中采用左對(duì)齊方式加窗,因此,由式(8)計(jì)算該脈沖的起始時(shí)刻n r,終止時(shí)刻n f和脈寬 n w。

因?yàn)樾盘?hào)在頻域具有能量集中的特性[9],遞歸算法采用頻域數(shù)據(jù)來(lái)測(cè)量時(shí)域參數(shù),相比用時(shí)域數(shù)據(jù)測(cè)量時(shí)域參數(shù),抗噪性能加強(qiáng)了,提高了檢測(cè)靈敏度。

3.3 參數(shù)N的設(shè)置與調(diào)整——影響時(shí)間分辨率、頻率分辨率、抗噪性能

雖然通過(guò)對(duì)頻譜細(xì)化系數(shù)α和數(shù)據(jù)幀的滑移量d的調(diào)整,可以使時(shí)頻網(wǎng)格的頻域間隔Δf grid和時(shí)域間隔Δtgrid達(dá)到很小,但是并不意味著頻率分辨率Δf(指能夠區(qū)分頻率軸上靠得很近的兩信號(hào)或兩頻率分量的最小頻率間隔)或時(shí)域分辨率Δt(指能夠區(qū)分時(shí)間軸上靠得很近的兩信號(hào)或時(shí)限信號(hào)的起始時(shí)刻、終止時(shí)刻的最小時(shí)間間隔)的改善。

遞歸運(yùn)算的頻率分辨率Δf和時(shí)間分辨率 Δt受采樣頻率f s、窗口寬度N(等于循環(huán)迭代次數(shù))、窗函數(shù)形狀等因素的影響[7]。Δf正比于f s/N,Δt反比于fs/N。N變大,頻域分辨率提高,而時(shí)域分辨率下降,二者仍然受制于“測(cè)不準(zhǔn)原理”。但在一定的范圍內(nèi),對(duì)循環(huán)迭代次數(shù) N的調(diào)整,使得遞歸運(yùn)算具有了可調(diào)可控的時(shí)間分辨率、頻率分辨率。

假設(shè)有一個(gè)分段正弦信號(hào)

其中:f 1=52 Hz,f 2=54 Hz,f 3=56 Hz,f 4=58 Hz,f s=200 Hz。當(dāng)工作參數(shù) f m=(50+0.01m)(Hz),m=0,1,…,99;d=10;N 分別取 100、200和400時(shí),利用遞歸算法得到的結(jié)果,如圖6。通過(guò)對(duì)比可以看到,N較小時(shí),具有較好的時(shí)間分辨率,而相應(yīng)的頻率分辨率則不高,隨著N的增加,頻率分辨率變得越高,但此時(shí)的時(shí)間分辨率則會(huì)相應(yīng)地下降。

另外,由圖6還可以看到,參數(shù) N除了影響時(shí)、頻分辨率外,還影響著譜線的幅度。N越大,譜線的幅度越大,則系統(tǒng)的抗噪性能越好。

4.結(jié) 論

遞歸算法具有很大的靈活性,非常適合對(duì)局部頻段實(shí)時(shí)地進(jìn)行頻域參數(shù)測(cè)量、時(shí)域參數(shù)測(cè)量和二維時(shí)頻譜分析,三項(xiàng)任務(wù)可獨(dú)立完成,也可相互配合完成。時(shí)頻分析系統(tǒng)需要根據(jù)某種規(guī)則做出工作參數(shù)的設(shè)置與調(diào)整,使其在保證實(shí)時(shí)性,時(shí)間分辨率,頻率分辨率,時(shí)、頻參數(shù)測(cè)量精度,抗噪性能等性能指標(biāo)時(shí)統(tǒng)籌調(diào)度通路資源,協(xié)調(diào)工作。

[1] 李岳霖,龐偉正,等.多普勒雷達(dá)測(cè)量輕武器彈丸轉(zhuǎn)速方法研究[J].電波科學(xué)學(xué)報(bào),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信號(hào)的截獲方法[J].電波科學(xué)學(xué)報(bào),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] 徐玉清,郎 銳,高 攀.基于小波的信息保持型雷場(chǎng)圖像融合算法[J].電波科學(xué)學(xué)報(bào),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è)計(jì)算法研究[D].陜西西安:西安電子科技大學(xué),2008.

[5] JACOBSEN E,LYONS R.The sliding DFT[J].IEEE Signal Processing Magazine,2003,20(3):74-80.

[6] 黃寒華.滑動(dòng)DFT算法研究[D].江蘇南京:東南大學(xué),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] 劉 平,靳成英,陳曾平.一種基于短時(shí)FFT的寬帶數(shù)字接收機(jī)設(shè)計(jì)[J].信號(hào)處理,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)