毛飛龍,焦義文,馬 宏,李 超,高澤夫,王育欣

(航天工程大學(xué) 電子與光學(xué)工程系,北京 101416)

0 引 言

隨著航天技術(shù)與深空探測(cè)技術(shù)的迅猛發(fā)展,各國(guó)對(duì)空間的探測(cè)向著更深更遠(yuǎn)的范圍拓展,探測(cè)距離將從目前月球探測(cè)的40萬(wàn)km拓展到火星探測(cè)的4億km和木星探測(cè)的10億km[1]。這將導(dǎo)致嚴(yán)重的傳輸損耗[2],同時(shí)也對(duì)深空探測(cè)航天器和地面設(shè)備提出了兩點(diǎn)要求:一是遠(yuǎn)距離通信能力;二是高速率傳輸能力[3]。由于深空探測(cè)航天器的等效全向發(fā)射功率(Effective Isotropic Radiated Power,EIRP)有限,為了實(shí)現(xiàn)對(duì)低信噪比信號(hào)的高質(zhì)量接收,就必須增大地面站的天線口徑。根據(jù)NASA的任務(wù)參數(shù)計(jì)算可知,火星探測(cè)任務(wù)中,傳輸1 Mb/s的信號(hào)需要直徑80 m的天線。而目前NASA深空網(wǎng)(Deep Space Network,DSN)的直徑70 m天線已接近工程極限,伺服機(jī)動(dòng)十分困難,靈活性較差。同時(shí),天線的口徑越大,其信號(hào)的波束越窄,越難以捕獲探測(cè)器的下行信號(hào)。因此,單天線的性能提升已經(jīng)到了緩慢發(fā)展接近停滯的狀態(tài),通過多天線合成的方法獲得更高的深空通信增益是值得大力發(fā)展的方法[4]。

隨著高速率通信系統(tǒng)的發(fā)展,信號(hào)的帶寬、傳輸速率在近年來也有了大幅度提升。信號(hào)帶寬可達(dá)幾十MHz甚至幾百M(fèi)Hz,這將為地面接收設(shè)備帶來巨大挑戰(zhàn)。寬帶信號(hào)多天線合成系統(tǒng)是一個(gè)數(shù)據(jù)密集、計(jì)算密集、邏輯復(fù)雜的計(jì)算系統(tǒng),使用單純的中央處理器(Central Processing Unit,CPU)計(jì)算已不能滿足運(yùn)算要求,需要選用新的計(jì)算架構(gòu)來提高系統(tǒng)實(shí)時(shí)處理的性能。近些年來,圖形處理單元(Graphics Processing Unit,GPU)在計(jì)算領(lǐng)域得到了充分應(yīng)用,已發(fā)展成為一種高度并行化、多線程、多核的通用計(jì)算設(shè)備,具有杰出的計(jì)算能力和極高的存儲(chǔ)器帶寬,并被越來越多地應(yīng)用于圖形處理之外的計(jì)算領(lǐng)域。基于GPU的信號(hào)處理技術(shù)成為眾多領(lǐng)域的熱點(diǎn),如射電天文、雷達(dá)、無(wú)線通信、人工智能等。因此,基于GPU的寬帶信號(hào)多天線合成具有一定的工程研究意義。多天線信號(hào)合成技術(shù)的本質(zhì)就是在消除不同天線信號(hào)間的時(shí)延差和相位差后進(jìn)行相干相加,因此天線間時(shí)延差和相位差的估計(jì)、時(shí)延差與相位差的補(bǔ)償為多天線合成中的兩項(xiàng)關(guān)鍵技術(shù)。本文旨在現(xiàn)有寬帶信號(hào)時(shí)延差與相位差估計(jì)方法的基礎(chǔ)上,設(shè)計(jì)基于GPU的寬帶信號(hào)時(shí)延差與相位差估計(jì)方法,通過算法的并行化映射尋求一種精度高且實(shí)時(shí)性好的估計(jì)方法。

1 寬帶信號(hào)估計(jì)方法分析

1.1 多天線信號(hào)合成方法

多天線信號(hào)合成的概念最早是由美國(guó)噴氣推進(jìn)實(shí)驗(yàn)室于1965年提出的[25]。該技術(shù)通過將不同天線接收到同一信源的信號(hào)進(jìn)行合成,從而提高接收信號(hào)的信噪比(Signal-to-Noise Ratio,SNR)。圖1為多天線信號(hào)合成原理示意圖。以雙天線為例,兩天線接收到來自同一信源的信號(hào),通過AD采集為數(shù)字信號(hào),然而由于各天線信號(hào)傳播路徑的不同,導(dǎo)致信號(hào)間存在各種各樣的差異,這些差異將損失信號(hào)間的相干性。此時(shí)兩天線信號(hào)存在時(shí)延差與相位差。將其中值較高的天線當(dāng)作參考天線,另一天線為待補(bǔ)償天線,對(duì)其進(jìn)行時(shí)延差與相位差估計(jì)并補(bǔ)償,將兩天線的時(shí)延與相位對(duì)齊。之后對(duì)兩天線信號(hào)相加,就得到了信號(hào)幅度翻倍、信噪比提升的合成信號(hào)。多天線信號(hào)合成技術(shù)的本質(zhì)就是在消除不同天線信號(hào)間的時(shí)延差和相位差后進(jìn)行相干相加。由于各天線間的噪聲是隨機(jī)非相關(guān)的,理論上使用N個(gè)天線合成信號(hào)的信噪比是單天線接收信號(hào)的N倍。

圖1 多天線信號(hào)合成原理示意

多天線信號(hào)合成技術(shù)經(jīng)過幾十年的發(fā)展,逐漸形成了5種合成方案[6]:全頻譜合成、基帶合成(Baseband Combining,BC)、符號(hào)流合成(Symbol-Stream Combining,SSC)、載波組陣合成(Carrier Arraying,CA)、復(fù)符號(hào)合成(Complex-Symbol Combining,CSC)[7]。除全頻譜合成方法之外,其余四種方法都屬于基于載波跟蹤技術(shù)的合成方法,都需要陣元載波鎖定后合成,不適用于當(dāng)前大規(guī)模的多天線信號(hào)合成。本文在考慮分析了各種方法的利弊后,選用全頻譜合成方法對(duì)寬帶信號(hào)進(jìn)行合成。

1.2 現(xiàn)有寬帶信號(hào)時(shí)延差與相位差估計(jì)方法

寬帶信號(hào)合成的關(guān)鍵是對(duì)寬帶信號(hào)間的殘余時(shí)延和相位進(jìn)行高精度估計(jì)和補(bǔ)償,目前可以分為時(shí)域合成方法和頻域合成方法。

1.2.1 時(shí)域合成方法

NASA深空網(wǎng)在其全頻譜處理陣(Full Spectrum Processing Array,FSPA)設(shè)備中,使用了一種基于上下邊帶的殘余時(shí)延和相位差聯(lián)合估計(jì)方法,原理如圖2所示。該方法首先通過提取信號(hào)濾波器(中心頻率可調(diào)帶通濾波器)提取各天線信號(hào)的上邊帶和下邊帶,然后利用相位估計(jì)算法分別估計(jì)待校準(zhǔn)天線與參考天線上下邊帶的相位差,最終兩相位點(diǎn)構(gòu)成的直線的斜率即為殘余時(shí)延,截距即為載波相位差。

圖2 基于上下邊帶的殘余時(shí)延和相位估計(jì)方法的原理

在圖2中,L1和L2表示兩天線信號(hào)經(jīng)過提取信號(hào)濾波器(中心頻率可調(diào)帶通濾波器)得到的下邊帶,其中心頻率為fL;U1和U2表示兩天線信號(hào)的上邊帶,其中心頻率為fU;相位差φL和φU為相位差估計(jì)算法估計(jì)得到的下邊帶間和上邊帶間的相位差。則兩相位點(diǎn)構(gòu)成直線的斜率即為殘余時(shí)延,截距即為載波處的殘余相位。設(shè)殘余時(shí)延為Δτ,載波處的殘余相位為Δφc,則有

(1)

(2)

該方法的優(yōu)勢(shì)在于運(yùn)算量小,而且將時(shí)延估計(jì)問題轉(zhuǎn)化為相位估計(jì)問題,在低信噪比條件下,各邊帶可以采用高性能相位估計(jì)算法來提高估計(jì)精度。然而,該方法存在兩個(gè)問題,下面進(jìn)行詳細(xì)分析。

1)相位卷繞問題

由于相位估計(jì)算法只能提取初相差,當(dāng)上下邊帶頻率差較大或者殘余時(shí)延較大時(shí)出現(xiàn)相位卷繞問題。而該方法只取了上下邊帶兩個(gè)相位差值進(jìn)行估計(jì),在沒有時(shí)延差與相位差先驗(yàn)知識(shí)的情況下無(wú)法正確地進(jìn)行解卷繞處理,進(jìn)而導(dǎo)致殘余時(shí)延和殘余相位與其真實(shí)值差距較大。

本文研究寬帶信號(hào)的合成方法,以典型的Ka頻段信號(hào)為例,其信號(hào)帶寬可達(dá)500 MHz,若采用基于上下邊帶的殘余時(shí)延和相位差聯(lián)合估計(jì)方法,在接收信號(hào)的上下邊帶分別進(jìn)行互相關(guān),當(dāng)上下邊帶頻率差較大或者殘余時(shí)延較大時(shí),得到上下邊帶的相位結(jié)果為其卷繞的值。因此,為了得到正確的殘余時(shí)延與殘余相位,就要保證不發(fā)生相位卷繞(φU-φL<2π),此時(shí)

(3)

當(dāng)Δf為250 MHz時(shí),Δτ為4×10-9s,即為4 ns,而Ka頻段的寬帶信號(hào)通常殘余時(shí)延為幾十ns甚至幾百ns,因此引入可消除相位卷繞現(xiàn)象的殘余時(shí)延與殘余相位估計(jì)方法是寬帶信號(hào)合成中的重點(diǎn)之一。

2)頻帶內(nèi)相位滑動(dòng)問題

整數(shù)時(shí)延補(bǔ)償都采用整數(shù)時(shí)鐘移位的補(bǔ)償方法,其最高精度為時(shí)鐘的周期Ts,整數(shù)時(shí)延補(bǔ)償之后必然存在Δτ

Δφ=2π·fU·Δτ-2π·fL·Δτ=2π·Δf·Δτ。

(4)

式中:Δφ為相位滑動(dòng);Δf為頻率高處與頻率低處的差值。相位滑動(dòng)現(xiàn)象將影響合成信號(hào)的質(zhì)量,當(dāng)Δφ=0時(shí),信號(hào)的帶寬較窄,可近似看作同一頻率,各處的相位補(bǔ)償值相同,信號(hào)合成的性能較好,如圖3所示;當(dāng)Δφ=π時(shí),相位補(bǔ)償為真實(shí)值的負(fù)值,將導(dǎo)致合成信號(hào)出現(xiàn)反向相消的現(xiàn)象,如圖4所示。

圖3 無(wú)相位滑動(dòng)現(xiàn)象示意圖

圖4 有相位滑動(dòng)現(xiàn)象示意圖

為了抑制反向相消現(xiàn)象,需保證Δφ<<π,即2π·Δf·Δτ<<π,若Δf為250 MHz,需滿足Δτ<<2 ns,同樣地,需要滿足采樣時(shí)鐘Ts<<2 ns,即采樣率需遠(yuǎn)大于500 MHz。考慮到信號(hào)帶寬為500 MHz時(shí),根據(jù)奈奎斯特采樣定理,只需采樣率大于1 000 MHz即可。而遠(yuǎn)大于500 MHz的采樣率將帶來采樣效率極大地浪費(fèi),同時(shí)也增加了系統(tǒng)實(shí)時(shí)處理的難度。因此,引入可消除相位滑動(dòng)現(xiàn)象的殘余時(shí)延與殘余相位估計(jì)方法也是寬帶信號(hào)合成中的重點(diǎn)之一。

由上述分析可知,時(shí)域合成方法都存在補(bǔ)償精度受采樣率限制的問題。即使假設(shè)時(shí)延估計(jì)誤差為零,補(bǔ)償后的小數(shù)時(shí)延(也稱作分?jǐn)?shù)時(shí)鐘時(shí)延,最大為半個(gè)采樣時(shí)鐘周期)也將導(dǎo)致頻帶內(nèi)的相位滑動(dòng),從而惡化合成性能。尤其當(dāng)帶寬較大時(shí),上述方法不再適用。因此,引入更合理的小數(shù)時(shí)延補(bǔ)償方法是寬帶信號(hào)合成的重點(diǎn)和難點(diǎn)問題之一。

1.2.2 頻域合成方法

頻域合成方法的處理流程為首先利用分析濾波器組對(duì)寬帶信號(hào)進(jìn)行信道化,得到低采樣速率的子帶信號(hào),再利用相位差估計(jì)技術(shù)得到不同天線子帶間的相位差,進(jìn)而估計(jì)寬帶信號(hào)間的殘余時(shí)延和相位并對(duì)其進(jìn)行補(bǔ)償,然后對(duì)各天線子帶進(jìn)行合成,最終利用綜合濾波器組對(duì)各天線合成后的子帶進(jìn)行重構(gòu),得到原始采樣率的寬帶合成信號(hào)。在對(duì)殘余時(shí)延進(jìn)行補(bǔ)償時(shí),首先采用數(shù)字延遲線進(jìn)行整數(shù)時(shí)延補(bǔ)償,殘余的小數(shù)時(shí)延以及相位的補(bǔ)償則在頻域進(jìn)行,通過對(duì)各個(gè)子帶進(jìn)行線性相位偏移來實(shí)現(xiàn)。相比較時(shí)域合成方法,頻域合成方法具有如下優(yōu)勢(shì):一是通過對(duì)子帶構(gòu)成的相頻曲線進(jìn)行解卷繞,可以有效解決時(shí)延差引起的2π模糊問題;二是通過將各子帶中心頻率獨(dú)立進(jìn)行相位對(duì)齊,可將由小數(shù)時(shí)延引起的整帶內(nèi)相位滑動(dòng)限制在子帶范圍內(nèi),大大減小了小數(shù)時(shí)延引起的合成損失。

由上述分析可知,頻域合成方法可以有效解決時(shí)域合成方法存在的問題,已得到廣泛關(guān)注和研究,是目前寬帶信號(hào)全頻譜合成的最佳實(shí)現(xiàn)方法。可以預(yù)見,該方法將是未來大規(guī)模寬帶天線組陣信號(hào)合成的主流方法。

1.3 本文采用的寬帶信號(hào)時(shí)延差與相位差估計(jì)方法

由1.2節(jié)分析可知,現(xiàn)有的基于上下邊帶的殘余時(shí)延和相位差聯(lián)合估計(jì)方法在處理本文中帶寬500 MHz的信號(hào)時(shí)存在兩方面問題,分別是相位卷繞問題和頻帶內(nèi)相位滑動(dòng)問題。對(duì)相位卷繞問題分析可知,其主要出現(xiàn)在上下邊帶頻率差較大或者殘余時(shí)延較大時(shí),而信號(hào)的殘余時(shí)延差Δτ無(wú)法改變,只能通過減小頻率差Δf來改進(jìn)。對(duì)頻帶內(nèi)相位滑動(dòng)問題分析可知,Δφ<<π時(shí),該現(xiàn)象可消除。而Δφ=2π·Δf·Δτ,同樣地,由于Δτ無(wú)法改變,只能通過減小頻率差Δf來減小相位滑動(dòng)問題帶來的影響。綜合上述因素,減小頻率差Δf是一種改進(jìn)思路,結(jié)合軟件無(wú)線電中信道化的思想,將各天線的信號(hào)拆分為多個(gè)獨(dú)立的子帶,之后對(duì)不同天線對(duì)應(yīng)的子帶信號(hào)進(jìn)行處理。下面進(jìn)行詳細(xì)分析:

將寬帶信號(hào)拆分為K個(gè)子帶,對(duì)相應(yīng)的子帶進(jìn)行相位差估計(jì),得到K個(gè)子帶對(duì)應(yīng)的相位差信息。此時(shí),只要保證相位差出現(xiàn)卷繞現(xiàn)象是由于某一子帶的相位差值在[-π,π]之外引起的,而并非相鄰兩個(gè)相位差的差值大于2π引起的,換句話說,只要保證相鄰兩個(gè)相位差不存在卷繞現(xiàn)象就可正確進(jìn)行解卷繞。在各子帶中,頻率間隔在數(shù)值上為fs/K,要使相鄰兩個(gè)相位差不存在卷繞現(xiàn)象,即

(5)

式中:τd為子帶內(nèi)的最大時(shí)延。式(5)可化簡(jiǎn)為

(6)

由式(3)可知,K需滿足

K>2fs·τd。

(7)

因此,當(dāng)時(shí)延差τd的精度確定時(shí),子帶的最小數(shù)量即確定。

對(duì)于頻帶內(nèi)相位滑動(dòng)問題,將寬帶信號(hào)拆為子帶信號(hào)進(jìn)行處理,可將相位滑動(dòng)限制在子帶范圍內(nèi),如圖5所示。可以看出,雖然沒有徹底消除相位滑動(dòng)問題,但相位滑動(dòng)被限制在了很低的范圍內(nèi),可大幅度降低相位滑動(dòng)問題對(duì)合成性能的影響,使合成信號(hào)的功率譜保持在可接收的范圍內(nèi),如圖6所示。

圖5 相位滑動(dòng)現(xiàn)象被限制在子帶內(nèi)示意

圖6 減弱相位滑動(dòng)前后合成信號(hào)功率譜示意

寬帶信號(hào)合成需考慮相位卷繞問題和頻帶內(nèi)相位滑動(dòng)問題,根據(jù)上述分析可知,將寬帶信號(hào)拆為子帶進(jìn)行處理可有效地降低這兩種問題帶來的影響。同時(shí),寬帶信號(hào)合成方案必須考慮相頻特性差異造成的相位非相干性以及采樣率高帶來的數(shù)據(jù)吞吐率大的問題。因此,本文選用寬帶信號(hào)頻域合成方法,如圖7所示。

圖7 寬帶信號(hào)頻域全頻譜合成方法

基本思路為首先根據(jù)頻率范圍對(duì)寬帶信號(hào)進(jìn)行子帶拆分,得到低采樣速率的子帶信號(hào)。利用相位差估計(jì)技術(shù)得到不同天線子帶信號(hào)間的相位差,進(jìn)而估計(jì)并補(bǔ)償寬帶信號(hào)間的殘余時(shí)延和相位,然后對(duì)各天線子帶進(jìn)行合成。最終通過對(duì)各天線合成后的子帶進(jìn)行重構(gòu),得到原始采樣率的寬帶合成信號(hào)。

寬帶信號(hào)時(shí)延差與相位差估計(jì)示意圖如圖8所示。以雙天線為例,天線信號(hào)經(jīng)AD采集與子帶拆分之后,每個(gè)天線信號(hào)得到了頻率范圍不同的K個(gè)子帶。對(duì)兩個(gè)天線對(duì)應(yīng)頻率的子帶信號(hào)進(jìn)行相位差估計(jì),得到K個(gè)子帶的相位差結(jié)果。由于寬帶信號(hào)的相位信息φ(ω)是頻率ω的線性函數(shù),且滿足

圖8 寬帶信號(hào)時(shí)延差與相位差估計(jì)示意

φ(ω)=ω·Δτk+θk。

(8)

則寬帶信號(hào)的殘余延遲差Δτk多項(xiàng)式為一次項(xiàng)參數(shù),即斜率,殘余相位差θk為零次項(xiàng)參數(shù),即截距。因此,在估計(jì)得到兩個(gè)天線子帶間的相位差信息后,通過線性最小二乘擬合就可以得到兩天線寬帶信號(hào)的殘余時(shí)延與殘余相位。由于子帶相位差估計(jì)方法最后需要進(jìn)行反正切運(yùn)算,因此,當(dāng)殘余時(shí)延對(duì)應(yīng)整周期內(nèi)的角度大于2π時(shí),就會(huì)出現(xiàn)相位卷繞現(xiàn)象。因此,需要在擬合之前進(jìn)行解卷繞得到正確的相位值。

圖9為基于GPU的寬帶信號(hào)時(shí)延差與相位差估計(jì)流程,具體步驟如下:

Step1 參考天線與待補(bǔ)償天線首先經(jīng)過FFT均勻信道化器將寬帶信號(hào)拆為多個(gè)子帶信號(hào)。

Step2 待補(bǔ)償天線的子帶信號(hào)均與參考天線對(duì)應(yīng)的子帶信號(hào)進(jìn)行共軛相乘運(yùn)算,完成互相關(guān)。

Step3 將多段互相關(guān)的結(jié)果進(jìn)行累加平均并求反正切,得到對(duì)應(yīng)子帶相位差的估計(jì)值。

Step4 待補(bǔ)償天線各子帶的相位值進(jìn)行解卷繞運(yùn)算得到正確的相位值。

Step5 對(duì)待補(bǔ)償天線所有子帶的相位差結(jié)果進(jìn)行最小二乘線性擬合,得到寬帶信號(hào)的時(shí)延差與相位差。

其中Step 1～3即為子帶相位差估計(jì)算法中典型的Simple算法。

2 GPU并行優(yōu)化

2.1 線程并行優(yōu)化

由1.3節(jié)可知,子帶相位差估計(jì)算法中Simple算法主要包括FFT、共軛相乘、累加平均三個(gè)模塊,下面對(duì)各模塊的GPU實(shí)現(xiàn)方法進(jìn)行研究。

在CPU平臺(tái)下,N點(diǎn)FFT的運(yùn)算量為N/2×lbN次的復(fù)數(shù)乘法和N×lbN次的復(fù)數(shù)加法,N表示做傅里葉變換的點(diǎn)數(shù)。對(duì)于如此大量的高速并行數(shù)據(jù),在CPU平臺(tái)中無(wú)法做到實(shí)時(shí)處理。但是由于FFT計(jì)算具有高度并行性這一顯著特點(diǎn),本文考慮將快速傅里葉變換部分放在GPU中進(jìn)行并行加速處理,利用CUDA中自帶的CUFFT庫(kù),可實(shí)現(xiàn)大量高速數(shù)據(jù)的實(shí)時(shí)FFT計(jì)算。調(diào)用CUFFT函數(shù)庫(kù)內(nèi)的cufftPlanD函數(shù)對(duì)數(shù)據(jù)進(jìn)行快速傅里葉變換。在GPU上進(jìn)行快速傅里葉變換需要完成以下幾步工作:

Step1 使用 malloc()函數(shù)為輸入的數(shù)據(jù)分配主機(jī)端(CPU)內(nèi)存,并存儲(chǔ)輸入數(shù)據(jù)。

Step2 使用 cudaMalloc()函數(shù)為輸入數(shù)據(jù)分配設(shè)備端(GPU)顯存。

Step3 使用 cudaMemcpy()函數(shù)將主機(jī)(CPU)內(nèi)存數(shù)據(jù)傳輸?shù)皆O(shè)備端(GPU)顯存。

Step4 創(chuàng)建一個(gè)plan,調(diào)用 cufftPlanD函數(shù)創(chuàng)建快速傅里葉變換。

Step5 執(zhí)行plan,使用 cufftExecC2C()函數(shù)完成plan的計(jì)算任務(wù)。

Step6 用cudaMemcpy()將GPU端顯存數(shù)據(jù)傳輸?shù)街鳈C(jī)內(nèi)存。

上述工作執(zhí)行完成以后如果這個(gè)plan不再需要,則調(diào)用cufftDestroy()函數(shù)銷毀該plan及分配給它的計(jì)算資源。以上各步的具體流程和所調(diào)用的CUDA函數(shù)如圖10所示。在CUFFT庫(kù)中,有3種不同的庫(kù)函數(shù):cufftExecC2C(復(fù)數(shù)到復(fù)數(shù))、cufftExecC2R(復(fù)數(shù)到實(shí)數(shù))、cufftExecR2C(實(shí)數(shù)到復(fù)數(shù)),它們的計(jì)算耗時(shí)由于本身運(yùn)算量的差異而有所不同。因此,為最大限度地縮短運(yùn)算時(shí)間,必須針對(duì)本提相過程的實(shí)際情況,選擇最合適的庫(kù)函數(shù)類型。在本系統(tǒng)輸入信號(hào)的FFT計(jì)算中,由于輸入是復(fù)信號(hào),而輸出是復(fù)數(shù)信號(hào),因此最宜采用cufftExecC2C類型的FFT。

圖11是GPU中FFT與共軛相乘模塊的實(shí)現(xiàn)框圖。調(diào)用cufft庫(kù)進(jìn)行FFT、多個(gè)線程同時(shí)進(jìn)行共軛相乘運(yùn)算,代碼如下:

圖11 FFT與共軛相乘GPU實(shí)現(xiàn)框圖

cufftPlan1d(&plan,NFFT,CUFFT_C2C,nBatch);

cufftExecC2C(plan,d_signal,d_signal_FFT,CUFFT_FORWARD);

cufftDestroy(plan);

Conj<<>>(d_signal_1_FFT,d_signal_4_FFT,d_Conj,nsize);

__global__ void Conj(float2 *src,float2 *send,float2 *dst,int length)

{ int tid = threadIdx.x + blockIdx.x * blockDim.x;

if(tid >= length)

return;

dst[tid].x =(src[tid].x * send[tid].x + src[tid].y * send[tid].y);

dst[tid].y =(src[tid].y * send[tid].x - src[tid].x * send[tid].y);}

圖12為累加平均并求相位核函數(shù)實(shí)現(xiàn)框圖。將共軛相乘后數(shù)據(jù)傳輸至累加平均并求相位核函數(shù)中進(jìn)行能量累積和提相等處理。將每個(gè)天線的N個(gè)數(shù)據(jù)等分成x段進(jìn)行相加平均,每段M/x點(diǎn),累加平均后使用atan 2函數(shù)求出M/x個(gè)復(fù)數(shù)的相角,所得相角即為子帶間的相位差。

圖12 累加平均并求相位GPU實(shí)現(xiàn)框圖

本文設(shè)計(jì)的累加平均并求相位核函數(shù)部分的CUDA代碼如下所示:

FAdd<<>>(d_Conj,d_phase,d_amp,fft_time,NFFT);

__global__ void FAdd(float2 *src,float *dst,float *dst2,int jft,int length)

{ int tid = threadIdx.x + blockIdx.x * blockDim.x;

if(tid >= length)

return;

int curpos;

double sum_x = 0.0;

double sum_y = 0.0;

for(int i = 0;i

{ curpos = tid + i * length;

sum_x = sum_x +((double)src[curpos].x);

sum_y = sum_y +((double)src[curpos].y);}

sum_x = sum_x / jft;

sum_y = sum_y / jft;

dst[tid]= atan2(sum_y,sum_x);

dst2[tid]= 10 * log10(sum_x * sum_x + sum_y * sum_y);}

2.2 異步流并發(fā)優(yōu)化

如圖9所示,基于GPU的寬帶信號(hào)時(shí)延差與相位差估計(jì)流程主要包括在GPU平臺(tái)上進(jìn)行FFT、子帶相位差估計(jì)、解卷繞、擬合等。子帶相位差估計(jì)具體包括共軛相乘、累加平均等運(yùn)算,FFT過程可用CUDA的庫(kù)函數(shù)CUFFT進(jìn)行計(jì)算,而CUFFT庫(kù)提供的二維FFT計(jì)算模型為本文的多路FFT高效并行實(shí)現(xiàn)提供了可能。

流(stream)是CUDA環(huán)境中的一個(gè)并行概念,在GPU中可使用流來處理并發(fā)操作[8]。流本質(zhì)上是一個(gè)操作隊(duì)列,該隊(duì)列可由主機(jī)上不同的線程發(fā)出并以異步方式在GPU上按順序執(zhí)行。不同的流可以并發(fā)執(zhí)行,多個(gè)流并發(fā)執(zhí)行不同路數(shù)據(jù)的FFT、子帶相位差估計(jì)、解卷繞、擬合等運(yùn)算,可大大提高運(yùn)算的并行性?；贕PU的寬帶信號(hào)時(shí)延差與相位差估計(jì),各路數(shù)據(jù)完全獨(dú)立,每一路數(shù)據(jù)的處理包含F(xiàn)FT、共軛相乘、累加平均、解卷繞、擬合等多個(gè)操作,因此,可進(jìn)一步利用異步流并發(fā)技術(shù)對(duì)估計(jì)方法進(jìn)行優(yōu)化。

設(shè)天線數(shù)為N,若將每一路天線數(shù)據(jù)的處理過程用一個(gè)CUDA流來處理,則基于CUDA流式架構(gòu)的相位差估計(jì)算法實(shí)現(xiàn)框圖如圖13所示。由圖中虛線框所示部分可知,在該算法中,同一時(shí)刻最多有N個(gè)Kernel在同時(shí)執(zhí)行。

圖13 基于CUDA流式架構(gòu)的相位差估計(jì)算法實(shí)現(xiàn)框圖

圖14為N=3時(shí)流式架構(gòu)相位差估計(jì)算法流程。從圖中可知,基于線程的并行方法中,數(shù)據(jù)傳輸和核函數(shù)順序執(zhí)行,即一次操作完全完成后才進(jìn)行下一項(xiàng)操作,在某次操作執(zhí)行期間其他所有操作都處于等待狀態(tài)。而在流式結(jié)構(gòu)并發(fā)算法中,數(shù)據(jù)傳輸和核函數(shù)的執(zhí)行被拆分為多個(gè)并行的操作序列,隨著CUDA流的并發(fā)啟動(dòng),數(shù)據(jù)傳輸過程被核函數(shù)執(zhí)行過程較好地隱藏。另外,從圖14的算法結(jié)構(gòu)對(duì)比結(jié)果可看出,3路并行流式結(jié)構(gòu)并行算法在執(zhí)行時(shí),同一時(shí)刻駐留在GPU上的核函數(shù)最多可達(dá)3個(gè),其執(zhí)行效率明顯高于線程并行算法。

圖14 基于CUDA流式架構(gòu)的4天線相位差估計(jì)算法實(shí)現(xiàn)框圖

2.3 小結(jié)

首先利用線程并行的方式對(duì)時(shí)延差與相位差估計(jì)方法中的共軛相乘、累加平均、解卷繞等模塊的實(shí)現(xiàn)算法進(jìn)行了相應(yīng)的改進(jìn)和優(yōu)化,實(shí)現(xiàn)了多個(gè)線程同時(shí)處理多個(gè)數(shù)據(jù)的運(yùn)算,提高了每個(gè)核函數(shù)本身的運(yùn)算效率。之后利用異步流并發(fā)的方式對(duì)估計(jì)方法進(jìn)一步優(yōu)化,每路信號(hào)采用獨(dú)立的流進(jìn)行運(yùn)算,從結(jié)構(gòu)層面有效提高了數(shù)據(jù)的并行處理效率。

3 仿真驗(yàn)證

本文采用Simple算法對(duì)寬帶信號(hào)進(jìn)行合成,利用相位干涉儀中的4個(gè)天線進(jìn)行實(shí)測(cè)數(shù)據(jù)信號(hào)合成。系統(tǒng)由4個(gè)陣元構(gòu)成,陣元擺放形式如圖15所示。系統(tǒng)硬件平臺(tái)選用NVIDIA Tesla V100顯卡,接收信號(hào)頻段采用S頻段,利用相位校正(Phase CALibration,PCAL)信號(hào)進(jìn)行性能分析,PCAL信號(hào)的頻率特性具有梳狀頻譜特性和線性相位特性。本系統(tǒng)的核心設(shè)計(jì)指標(biāo)如表1所示。

表1 仿真實(shí)驗(yàn)參數(shù)

圖15 相位干涉儀布陣示意

3.1 正確性驗(yàn)證

首先驗(yàn)證GPU線程并行與GPU流并行估計(jì)方法的正確性,設(shè)置天線4為參考天線,天線1與天線4的時(shí)延差為2個(gè)采樣點(diǎn)、相位差為-π/2。利用CPU串行Simple算法、GPU線程并行Simple算法與GPU流并行Simple算法估計(jì)天線1與參考天線的時(shí)延差與相位差。每次估計(jì)出時(shí)延差與相位差后對(duì)天線1信號(hào)進(jìn)行修正,按式(9)更新時(shí)延和相位補(bǔ)償值:

(9)

式中:delaynow和phasenow為本次迭代計(jì)算出的時(shí)延差與相位差;delayi-1和phasei-1為上一次迭代的時(shí)延和相位補(bǔ)償值;delayi和phasei為本次更新后的時(shí)延補(bǔ)償值。步長(zhǎng)因子為0.5,記錄迭代30次的時(shí)延和相位補(bǔ)償值。

圖16為CPU串行Simple算法、GPU線程并行Simple算法與GPU流并行Simple算法的時(shí)延和相位補(bǔ)償值隨迭代次數(shù)變化的曲線。20次迭代后,三種算法的時(shí)延和相位補(bǔ)償值都收斂于真值,都可實(shí)現(xiàn)高精度時(shí)延和相位估計(jì),驗(yàn)證了GPU線程并行與GPU流并行估計(jì)方法的正確性。

利用本文中的并行Simple算法與傳統(tǒng)串行Simple算法在不同單天線信噪比的條件下分別進(jìn)行合成,進(jìn)行蒙特卡洛仿真并得出合成損失與合成效率的平均值,得到不同信噪比條件下兩種方法的合成損失與合成效率。圖17給出了兩種方法在不同單天線信噪比條件下合成損失與合成效率的仿真結(jié)果。

(a)合成損失

從圖17可以看出,本文中的并行Simple算法與傳統(tǒng)串行Simple算法的合成損失與合成效率曲線重合,兩者的合成性能一致,與分析得出的結(jié)論一致。同時(shí),兩種算法隨著單天線信噪比的升高,合成損失降低,合成效率升高。仿真結(jié)果表明,本文中的并行Simple算法與傳統(tǒng)串行Simple算法在信噪比高于-15 dB的條件下有較好的合成性能。

3.2 加速效果驗(yàn)證

為了詳細(xì)分析核函數(shù)運(yùn)行過程,利用NVIDIA Nsight Systems對(duì)流式架構(gòu)Simple算法進(jìn)行分析,結(jié)果如圖18和圖19所示。

圖18 CUDA流式架構(gòu)Simple算法時(shí)序圖

圖19 CUDA流式架構(gòu)Sumple算法時(shí)序圖

圖18中上方為未采用異步流并發(fā)方式優(yōu)化的Simple算法時(shí)序圖,可以看到同一時(shí)刻最多只有一個(gè)核函數(shù)駐留,不同天線的相同運(yùn)算模塊順序運(yùn)行。圖18中下方為采用異步流并發(fā)方式優(yōu)化的Simple算法時(shí)序圖。由于Simple算法將其中1個(gè)天線設(shè)置為參考天線,因此,在天線數(shù)為4的情況下只需進(jìn)行3路信號(hào)流并發(fā)即可。從圖18中可以看到同一時(shí)刻最多有3個(gè)核函數(shù)駐留,不同天線的相同運(yùn)算模塊分別在3個(gè)不同流上執(zhí)行。

圖19中上方為未采用異步流并發(fā)方式優(yōu)化的Sumple算法時(shí)序圖,可以看到同一時(shí)刻最多只有一個(gè)核函數(shù)駐留,不同天線的相同運(yùn)算模塊順序運(yùn)行。圖19中下方為采用異步流并發(fā)方式優(yōu)化的Sumple算法時(shí)序圖。由于Sumple算法將所有天線的加權(quán)值設(shè)置為參考天線,每路天線信號(hào)都要與虛擬的參考天線進(jìn)行時(shí)延差與相位差估計(jì),因此,在天線數(shù)為4的情況下需要進(jìn)行4路信號(hào)的流并發(fā)。從圖19中可以看到同一時(shí)刻最多有4個(gè)核函數(shù)駐留,不同天線的相同運(yùn)算模塊分別在4個(gè)不同流上執(zhí)行。從時(shí)序圖中可明顯地看出異步流并發(fā)適合于時(shí)延差與相位差估計(jì)方法的優(yōu)化加速。

下面對(duì)GPU線程并行與GPU流并行估計(jì)方法相比CPU串行方法的加速比進(jìn)行測(cè)試,分別對(duì)30次迭代中CPU串行、GPU線程并行與GPU流并行估計(jì)方法的耗時(shí)進(jìn)行記錄。共測(cè)試4組,每組的區(qū)別只有迭代的數(shù)據(jù)量,分別為512,5 120,51 200,512 000。表2加速比測(cè)試結(jié)果。

表2 加速比測(cè)試結(jié)果

圖20為GPU線程并行與GPU流并行估計(jì)方法相比CPU串行方法的加速比測(cè)試結(jié)果圖,可以看出,當(dāng)數(shù)據(jù)量為512時(shí),GPU線程并行方法相比CPU串行方法的加速比為0.98,表明此時(shí)的GPU線程并行方法的平均耗時(shí)大于CPU串行方法的平均耗時(shí)。因?yàn)楫?dāng)數(shù)據(jù)量較小時(shí),每個(gè)核函數(shù)調(diào)用的線程數(shù)較少,不能充分地發(fā)揮GPU多線程的優(yōu)勢(shì),因此沒有加速效果。GPU流并行方法相比CPU串行方法的加速比為2.01,雖然數(shù)據(jù)量較小時(shí)不能發(fā)揮GPU多線程的優(yōu)勢(shì),但是流并行方法在結(jié)構(gòu)上對(duì)運(yùn)算順序做了優(yōu)化可以利用多流并發(fā)的方式同時(shí)對(duì)多路信號(hào)進(jìn)行運(yùn)算,因此加速效果明顯。隨著數(shù)據(jù)量的升高,GPU線程并行與GPU流并行估計(jì)方法相比CPU串行方法的加速比也升高,表明GPU線程并行與GPU流并行估計(jì)方法適合大規(guī)模數(shù)據(jù)的實(shí)時(shí)處理。從表2可以看出,CPU串行、GPU線程并行、GPU流并行估計(jì)方法在數(shù)據(jù)量為512 000的估計(jì)耗時(shí)約為17.43 ms,0.33 ms與0.14 ms。本文中系統(tǒng)的采樣率為56 000 000,處理512 000數(shù)據(jù)量耗時(shí)若小于10 ms即可實(shí)現(xiàn)實(shí)時(shí)估計(jì)。因此,CPU串行估計(jì)方法無(wú)法實(shí)現(xiàn)實(shí)時(shí)估計(jì),本文中的GPU線程并行、GPU流并行估計(jì)方法可實(shí)現(xiàn)實(shí)時(shí)估計(jì)。此外,本文中的數(shù)據(jù)為仿真數(shù)據(jù),提前讀入GPU,而實(shí)際系統(tǒng)中,受限于系統(tǒng)吞吐率與GPU并行性能的限制,實(shí)際系統(tǒng)的加速比不能無(wú)限增大。

圖20 加速比測(cè)試結(jié)果

3 結(jié) 論

針對(duì)大規(guī)模多天線信號(hào)合成系統(tǒng)的天線數(shù)量越來越多,接收信號(hào)的帶寬越來越大,單獨(dú)使用CPU已不能滿足運(yùn)算需求這一問題,本文研究了寬帶信號(hào)時(shí)延差與相位差估計(jì)方法,分析了算法并行的可行性并設(shè)計(jì)了基于GPU的寬帶信號(hào)時(shí)延差與相位差估計(jì)方法,通過仿真對(duì)比實(shí)驗(yàn)驗(yàn)證了本文提出的并行算法的優(yōu)越性。

本文的并行算法可以運(yùn)用到未來的大規(guī)模寬帶信號(hào)多天線信號(hào)合成系統(tǒng)中,實(shí)現(xiàn)高速、準(zhǔn)確、實(shí)時(shí)的信號(hào)合成,有力支持未來的一系列深空探測(cè)任務(wù)。