陳永強(qiáng)，馬 宏，黨宏杰，焦義文，劉燕都

(1.航天工程大學(xué) 電子與光學(xué)工程系，北京 101416；2.北京通信與跟蹤技術(shù)研究所，北京 100094)

0 引言

多相濾波器組(Polyphase Filter Bank, PFB)是數(shù)字信號(hào)濾波抽取的一種高效實(shí)現(xiàn)結(jié)構(gòu)，利用該技術(shù)可將串行的信道化過(guò)程分解為并行的多路處理流程，從而提高數(shù)字信號(hào)處理效率[1]。1991年,SETI首次將PFB引入射電系統(tǒng)并使用該技術(shù)設(shè)計(jì)了頻譜分析儀[2]，此后，PFB被越來(lái)越多地用于信號(hào)處理數(shù)字后端，目前國(guó)際上多個(gè)重要系統(tǒng)均使用了PFB技術(shù)進(jìn)行觀測(cè)頻譜的信道化[3]。

20世紀(jì)60年代,美國(guó)Geraid.Estrin[4]率先提出可重構(gòu)計(jì)算(Reconfigurable Computing)概念。該技術(shù)核心思想在于，在通用平臺(tái)上通過(guò)系統(tǒng)軟硬件結(jié)構(gòu)的靈活重構(gòu)實(shí)現(xiàn)不同的功能[5-6]。圖形處理器件(Graphic Processing Unit, GPU)由于具有眾多的運(yùn)算核心，特別適合于大量數(shù)據(jù)的并行處理，2007年NVIDIA提出的計(jì)算統(tǒng)一設(shè)備架構(gòu)(Compute Unified Device Architecture, CUDA)從軟件和硬件層面大大簡(jiǎn)化了基于GPU的系統(tǒng)開(kāi)發(fā)流程，使得GPU在通用計(jì)算領(lǐng)域得到更為廣泛應(yīng)用。當(dāng)前，由于GPU相較于FPGA能夠?qū)崿F(xiàn)較高的頻譜高分辨率且具有更高的可重構(gòu)性和擴(kuò)展性[4]，基于CPU+GPU的異構(gòu)信號(hào)處理系統(tǒng)在射電天文[7-8]、雷達(dá)[9-10]、無(wú)線通信[11-12]等眾多領(lǐng)域成為研究熱點(diǎn)。

隨著基于GPU的高性能計(jì)算(High Performance Computation,HPC)技術(shù)的快速發(fā)展，基于GPU的PFB系統(tǒng)研究逐漸受到研發(fā)人員的重視。2011年，麻省理工學(xué)院Haystack天文臺(tái)的Mark[13]使用NVIDIA Tesla C2050 GPU設(shè)計(jì)了一款用于VLBI 數(shù)字采集系統(tǒng)RDBE的PFB系統(tǒng)，驗(yàn)證了用GPU替代現(xiàn)有的FPGA板信號(hào)處理系統(tǒng)的可能性。測(cè)試結(jié)果顯示，系統(tǒng)實(shí)時(shí)數(shù)據(jù)處理速率達(dá)到890 MB/s，而且隨著GPU技術(shù)的進(jìn)一步發(fā)展，系統(tǒng)的處理能力將進(jìn)一步提升和擴(kuò)展。2014-2016年，馬里蘭大學(xué)的Scott C. Kim等研究了基于GPU的多載波系統(tǒng)低延遲多速率重信道化器[14-16]和寬帶接收機(jī)[12]，該團(tuán)隊(duì)利用GPU多層次線程結(jié)構(gòu)和存儲(chǔ)結(jié)構(gòu)對(duì)PFB的實(shí)現(xiàn)進(jìn)行了優(yōu)化，采用時(shí)域卷積和高維線程模型，在數(shù)十兆赫茲帶寬內(nèi)實(shí)現(xiàn)了2G、3G和4G無(wú)線電通信信號(hào)的高效信道化，系統(tǒng)具有設(shè)計(jì)靈活、軟件重構(gòu)、低延遲和高數(shù)據(jù)吞吐量等優(yōu)點(diǎn)。2017-2018年，Simon Faulkner等基于GPU開(kāi)發(fā)了一款射頻頻譜感知截獲系統(tǒng)，并為該系統(tǒng)研發(fā)了信道化設(shè)備[17-18]。該系統(tǒng)采用多相濾波結(jié)構(gòu)實(shí)現(xiàn)寬帶信號(hào)的信道化接收，采用CUDA stream對(duì)流程進(jìn)行了并發(fā)優(yōu)化。最終，系統(tǒng)處理帶寬達(dá)到500 MHz帶寬，可實(shí)現(xiàn)1.333 Gs/s采樣率的實(shí)時(shí)信道化處理能力。該文獻(xiàn)提出的信道化器是一種比較成熟的PFB結(jié)構(gòu)。2017年，新疆天文臺(tái)[8]為其觀測(cè)系統(tǒng)設(shè)計(jì)了基于GPU的多相濾波系統(tǒng)，為GPU在國(guó)內(nèi)數(shù)字后端應(yīng)用進(jìn)行了有效的探索。

在多相濾波系統(tǒng)中，F(xiàn)IR濾波環(huán)節(jié)由于涉及大量數(shù)據(jù)的乘加運(yùn)算，是影響系統(tǒng)效率的關(guān)鍵因素。當(dāng)前，在GPU平臺(tái)上FIR濾波算法的實(shí)現(xiàn)方式主要分為頻域?yàn)V波算法[19]和時(shí)域?yàn)V波算法[20]2種。文獻(xiàn)[21-22]給出了基于FFT的頻域FIR濾波方法的實(shí)現(xiàn)方式，得到了較好的加速比,同時(shí)對(duì)比了GPU、Intel-ipps和FFTW三個(gè)平臺(tái)上長(zhǎng)序列頻域?yàn)V波效率，證明Intel-ipps性能略?xún)?yōu)于FFTW。而Scott C. Kim[12]、Mark[13]及Jayanth Chennamangalam[23]等均從時(shí)域?qū)崿F(xiàn)了多相濾波過(guò)程，給出了時(shí)域FIR在連續(xù)信號(hào)濾波中的應(yīng)用。以上文獻(xiàn)分別從時(shí)域和頻域給出了長(zhǎng)序列FIR濾波的優(yōu)勢(shì)和實(shí)現(xiàn)方法，卻沒(méi)有給出在不同的應(yīng)用場(chǎng)景下2種濾波算法的適用條件。另外，近年來(lái)，在基于FPGA的平臺(tái)上，基于DA法[24-25]和查找表法[26]等多種優(yōu)化方法的高效并行FIR濾波方法研究取得豐碩成果，但這些方法難以直接移植到GPU平臺(tái)。

本文首先介紹了基于多相濾波技術(shù)的并行信道化算法，并分析了其運(yùn)算效率；然后對(duì)運(yùn)算過(guò)程中耗時(shí)最長(zhǎng)的多通道并行濾波過(guò)程進(jìn)行了分析，基于CUDA流式架構(gòu)分別設(shè)計(jì)了基于時(shí)域卷積和頻域快速卷積的FIR濾波算法，并分析了2種算法在多相信道化結(jié)構(gòu)中的性能；最后基于GPU平臺(tái)設(shè)計(jì)了多相信道化實(shí)現(xiàn)方法并用實(shí)驗(yàn)驗(yàn)證了分析的正確性。

1 基于多相濾波技術(shù)的并行信道化算法

典型的K通道并行信道化算法低通濾波實(shí)現(xiàn)原理框圖如圖1所示[27]。

圖1 多通道基帶轉(zhuǎn)換原理框圖Fig.1 Block diagram of multi-channel baseband conversion

圖1中，各信道中心頻率為ωk=2πfk/fs，由于實(shí)際處理過(guò)程中所用均為實(shí)信號(hào)，因此本文重點(diǎn)針對(duì)實(shí)信號(hào)進(jìn)行分析。設(shè)經(jīng)過(guò)D倍抽取后,輸出采樣率為fs2=fs/D，采樣周期Ts2=DTs1，在實(shí)信號(hào)模式下，該結(jié)構(gòu)第k通道數(shù)學(xué)表達(dá)式為：

yk(mTs2)=x1(nTs1)e-jωkn*hLP(nTs1)|n=mD|=

(1)

為了提高運(yùn)算效率，對(duì)濾波過(guò)程進(jìn)行多相分解，將濾波器分解為并行的K路，每一路分支濾波器長(zhǎng)L=[N/K]。原型低通濾波器索引i可分解為：

i=qK+p，

(2)

式中，q為每一路分支濾波器內(nèi)部點(diǎn)數(shù)索引，q=0,1,2,…,L-1；p為分支濾波器索引，p=0,1,2,…,K-1?？傻茫?/p>

(3)

令hp(m)=hLP((mK+p)Ts1)表示第p路分支濾波器，xp(m)=x1[(mD-p)Ts1]表示并行分路后第p路輸入信號(hào)。式(3)可改寫(xiě)為：

(4)

令h=K/D，i=q,l=qh,則有：

(5)

(6)

此時(shí)，實(shí)際上按照多相分支濾波器架構(gòu)完成了中心頻率為ωk的射頻信號(hào)的單通道信道化濾波接收。由于各通道均勻劃分，可令ωk=2πk/D，濾波器截至頻率為π/D，則第k通道信道化后的信號(hào)可表示為：

(7)

式中，為了最大限度提高處理效率，降低每個(gè)通道的數(shù)據(jù)速率，令抽取倍數(shù)D等于通道數(shù)K，即可得到最大抽取的DFT濾波器組：

(8)

式(8)即為最大抽取DFT濾波器組的經(jīng)典多相信道化結(jié)構(gòu)的數(shù)學(xué)表示，也是干涉測(cè)量基帶轉(zhuǎn)換器、射電天文數(shù)字后端等系統(tǒng)常用信道化算法。該算法本質(zhì)是多通道并行分支濾波和多路DFT，DFT可用其快速算法FFT實(shí)現(xiàn)，因此多路分支濾波運(yùn)算效率成為制約該結(jié)構(gòu)實(shí)時(shí)性的主要因素。由式(8)可知，分支濾波過(guò)程本質(zhì)上仍然是多路FIR濾波過(guò)程，而對(duì)于長(zhǎng)度為L(zhǎng)的輸入序列x(n),n=0,1,…,L-1和長(zhǎng)度為M的FIR濾波器h(n),n=0,1,…,M-1，F(xiàn)IR濾波過(guò)程可用如下線性卷積關(guān)系表示：

(9)

式中，y(n)為線性卷積輸出，其長(zhǎng)度為L(zhǎng)+M-1。此時(shí)，問(wèn)題轉(zhuǎn)化為長(zhǎng)序列FIR濾波的高效實(shí)現(xiàn)。

2 基于CUDA的長(zhǎng)序列實(shí)時(shí)FIR濾波方法研究

在實(shí)際工程實(shí)踐中，由于輸入信號(hào)過(guò)長(zhǎng)無(wú)法一次處理，通常需要做分段處理，但由于卷積運(yùn)算本身特殊性，分段處理時(shí)在每一段數(shù)據(jù)兩端位置將出現(xiàn)結(jié)果異常。為解決這一問(wèn)題，研究者提出基于重疊相加法和重疊保留法[28]的長(zhǎng)序列FIR濾波連續(xù)方法。馬里蘭大學(xué)的Kim等[12]利用卷積算法和重疊保留法在CUDA平臺(tái)上實(shí)現(xiàn)了高效濾波運(yùn)算，設(shè)計(jì)了基于GPU的寬帶信道化接收機(jī)。為了提高卷積運(yùn)算的效率，文獻(xiàn)[21-22]也提出了重疊相加法和頻域卷積相結(jié)合的快速濾波方法，并給出了該方法在CUDA平臺(tái)上的應(yīng)用效果。本節(jié)將重點(diǎn)根據(jù)并行信道化運(yùn)算需求，設(shè)計(jì)適用于本文架構(gòu)的高效濾波方法。

2.1 典型長(zhǎng)序列濾波算法

為了實(shí)現(xiàn)長(zhǎng)序列的連續(xù)濾波，需要對(duì)輸入數(shù)據(jù)分段進(jìn)行運(yùn)算。而根據(jù)分段處理方法的不同，常用的長(zhǎng)序列濾波方法分為重疊相加法和重疊保留法[29]。2種算法實(shí)現(xiàn)流程如下。

2.1.1 基于重疊相加法的長(zhǎng)序列分段濾波方法

設(shè)長(zhǎng)為N，輸入數(shù)據(jù)x(n)每一段處理長(zhǎng)度為L(zhǎng)，濾波器h(n)階數(shù)為M-1，那么該算法的執(zhí)行流程為：

① 將輸入數(shù)據(jù)分段，每段長(zhǎng)度為L(zhǎng)；

② 計(jì)算第1段數(shù)據(jù)與h(n)的卷積，得到濾波結(jié)果y0(n)長(zhǎng)度為L(zhǎng)+M-1；

③ 計(jì)算第2段數(shù)據(jù)與h(n)的卷積，得到濾波結(jié)果y1(n)長(zhǎng)度為L(zhǎng)+M-1；

④ 將y1(n)與y0(n)拼接作為輸出，并使得y1(n)的前M-1點(diǎn)與y0(n)的后M-1點(diǎn)重疊相加；

⑤ 以拼接結(jié)果作為新的y0(n)并重復(fù)以上兩步，直到分段數(shù)據(jù)處理完畢。

2.1.2 基于重疊保留法的長(zhǎng)序列分段濾波方法

設(shè)長(zhǎng)為N，輸入數(shù)據(jù)x(n)每一段處理長(zhǎng)度為L(zhǎng)，濾波器h(n)階數(shù)為M-1，那么該算法的執(zhí)行流程為：

① 將輸入數(shù)據(jù)分段，每段長(zhǎng)度為L(zhǎng)；

⑥ 以y0(n)，y1(n)，y2(n)，…的順序拼接結(jié)果即可得到濾波輸出。

重疊相加法的優(yōu)勢(shì)在于對(duì)輸入數(shù)據(jù)的操作簡(jiǎn)單，在分段方法確定后即可直接進(jìn)行卷積運(yùn)算，流程清晰，但完成卷積運(yùn)算后各段結(jié)果數(shù)據(jù)之間增加了加法操作。重疊保留法需要對(duì)輸入數(shù)據(jù)進(jìn)行M-1點(diǎn)的循環(huán)拷貝，增加了數(shù)據(jù)傳輸?shù)膲毫?，但其輸出?shù)據(jù)直接為最終結(jié)果，無(wú)需額外操作。2種方法雖略有區(qū)別，但其綜合運(yùn)算復(fù)雜度相當(dāng)[30]，本文選擇重疊保留法作為數(shù)據(jù)處理方法。

2.2 基于GPU的多通道并行時(shí)域?yàn)V波方法

由式(9)可知，對(duì)一個(gè)時(shí)間序列做FIR濾波本質(zhì)上就是將該信號(hào)與濾波器單位沖激響應(yīng)做線性卷積。輸出信號(hào)實(shí)際上是以濾波器單位沖擊響應(yīng)為權(quán)值，對(duì)輸入信號(hào)滑動(dòng)求取加權(quán)和。在實(shí)時(shí)信號(hào)處理中，這一過(guò)程可以通過(guò)重疊保留法予以實(shí)現(xiàn)。GPU架構(gòu)下基于重疊保留法的多通道并行FIR實(shí)現(xiàn)過(guò)程如圖2所示。

圖2 基于重疊保留法的時(shí)域FIR濾波流程Fig.2 Time-domain FIR filtering process based on overlap preservation method

GPU架構(gòu)下基于重疊保留法的長(zhǎng)序列時(shí)域FIR濾波流程如下：

① 設(shè)參與濾波的數(shù)據(jù)為k路，定義并啟動(dòng)k路CUDA stream，保證每一個(gè)stream對(duì)應(yīng)處理一路數(shù)據(jù)；

② 在每一個(gè)stream內(nèi)，將輸入采樣信號(hào)x(n)用重疊保留法分塊為xi(n)，分塊后每一段數(shù)據(jù)大小為L(zhǎng)，濾波器長(zhǎng)度為M，然后將第1段數(shù)據(jù)前向擴(kuò)展M-1個(gè)點(diǎn)，使得其中擴(kuò)展后數(shù)據(jù)長(zhǎng)度為L(zhǎng)+M-1；

③ 利用線性卷積分別計(jì)算第1段前L次滑動(dòng)卷積過(guò)程，將L點(diǎn)結(jié)果輸出至顯存；

④ 將上一段后M-1點(diǎn)數(shù)據(jù)拷貝到下一段數(shù)據(jù)頭部重新組成L+M-1點(diǎn)數(shù)據(jù)；

⑤ 重復(fù)③、④兩步，計(jì)算每一段分塊數(shù)據(jù)時(shí)域?yàn)V波結(jié)果yi(n)；

⑥ 將yi(n)各段拼接整合，提取出整段數(shù)據(jù)最終的濾波結(jié)果。

2.3 基于GPU的多通道并行頻域?yàn)V波方法

由傅里葉變換的原理可知，2個(gè)序列的DFT的乘積相當(dāng)于該序列時(shí)域做循環(huán)卷積(或圓周卷積)。而根據(jù)文獻(xiàn)[33]利用圓周卷積無(wú)混疊計(jì)算線性卷積的條件是，輸出信號(hào)y(n)至少需要N點(diǎn)DFT，N≥L+M-1，即：

Y(k)=X(k)H(k),k=0,1,...,N-1，

(10)

式中，X(k)，H(k)對(duì)應(yīng)于x(n)和h(n)的N點(diǎn)DFT。由于輸出信號(hào)y(n)的N點(diǎn)DFT一定可以在頻域表示該信號(hào)，故利用DFT先求的輸入信號(hào)和濾波器系數(shù)的N點(diǎn)DFT，然后在頻域逐點(diǎn)相乘得到乘積序列Y(k)，最后對(duì)Y(k)求N點(diǎn)DFT得到的圓周卷積，該圓周卷積結(jié)果等于x(n)和h(n)的線性卷積，即：

y(n)=x(n)*h(n)=IDFT[X(k)H(k)]=

IDFT[Y(k)]，

(11)

式中，k=0,1,...,N-1；n=0,1,...,N-1。

以上方法給出了利用DFT實(shí)現(xiàn)線性卷積的過(guò)程，從公式和數(shù)據(jù)處理流程上看，與時(shí)域直接卷積方法相比，該方法增加了信號(hào)和濾波器系數(shù)時(shí)域擴(kuò)展、DFT運(yùn)算、頻域乘法運(yùn)算以及IDFT運(yùn)算，流程更加復(fù)雜，但是與復(fù)雜的卷積運(yùn)算相比，DFT和IDFT可以通過(guò)快速算法得到所需信號(hào)，極大地提高了運(yùn)算效率?；谥丿B保留法的頻域FIR實(shí)現(xiàn)過(guò)程如圖3所示。

圖3 基于重疊保留法的頻域FIR濾波流程Fig.3 Frequency-domain FIR filtering process based on overlap preservation method

GPU平臺(tái)上基于重疊保留法的頻域快速卷積算法流程如下：

① 設(shè)參與濾波的數(shù)據(jù)為k路，則定義并啟動(dòng)k路CUDA stream和FFT句柄，并將Cufft句柄綁定到stream；

② 將輸入采樣信號(hào)x(n)分段，第i段為xi(n)；每一段數(shù)據(jù)長(zhǎng)度為L(zhǎng)，濾波器h(n)長(zhǎng)度為M;

③ 將每一段數(shù)據(jù)xi(n)和濾波器系數(shù)h(n)擴(kuò)展為L(zhǎng)+M-1位。其中數(shù)據(jù)向前擴(kuò)展M-1位，即在x0(n)前M-1位補(bǔ)零，此后在xi(n)前M-1位補(bǔ)xi-1(n)的后M-1位；在h(n)后向擴(kuò)展L-1位并補(bǔ)零；

④ 利用Cufft庫(kù)函數(shù)計(jì)算分段數(shù)據(jù)L+M-1點(diǎn)FFT，Xi(k)=FFT[xi(n)]，同時(shí)對(duì)濾波器也做相同點(diǎn)數(shù)的FFT，Hi(k)=FFT[h(n)]；

⑤ 將經(jīng)過(guò)FFT運(yùn)算的數(shù)據(jù)Xi(k)和濾波器系數(shù)Hi(k)對(duì)應(yīng)相乘，得到每一段數(shù)據(jù)的濾波結(jié)果的頻譜，Yi(k)=Xi(k)Hi(k)；

⑥ 對(duì)濾波結(jié)果的頻譜做L+M-1點(diǎn)FFT，yi(n)=IFFT[Yi(k)]，即可得到分段數(shù)據(jù)長(zhǎng)度為L(zhǎng)+M-1的線性卷積結(jié)果yi(n)；

⑦ 將yi(n)中后L點(diǎn)數(shù)據(jù)取出并按順序拼接，即可得到整段數(shù)據(jù)的濾波結(jié)果。

2.4 2種濾波方法對(duì)比分析及實(shí)現(xiàn)

基于重疊保留法的頻域?yàn)V波過(guò)程有效利用了FFT算法優(yōu)勢(shì)，在長(zhǎng)序列高階濾波過(guò)程中加速效果明顯。然而，該方法將原本簡(jiǎn)單的乘加運(yùn)算關(guān)系變成了變量擴(kuò)展、分段FFT、相乘、分段IFFT等多個(gè)步驟，增加了流程的復(fù)雜度，給流程的調(diào)度帶來(lái)了額外的開(kāi)銷(xiāo)，而且在濾波器階數(shù)較少的情況下，在重疊部分需要額外增加大量的內(nèi)存操作和運(yùn)算，這些操作將直接影響算法性能；另外，為了使用FFT算法加速運(yùn)算過(guò)程，數(shù)據(jù)長(zhǎng)度也需滿(mǎn)足2的整數(shù)次冪的要求。與之相比，時(shí)域?yàn)V波方法由于簡(jiǎn)單直接，在低階運(yùn)算中將具有更優(yōu)的應(yīng)用效果，但卷積運(yùn)算本身的復(fù)雜性將導(dǎo)致時(shí)域算法在高階濾波運(yùn)算中效率急劇下降。下面用仿真方法對(duì)2種算法的運(yùn)算復(fù)雜度進(jìn)行分析。

設(shè)每一路輸入數(shù)據(jù)分段長(zhǎng)度為L(zhǎng)，濾波器系數(shù)長(zhǎng)度M，基于重疊保留法的時(shí)域和頻域?yàn)V波算法的計(jì)算復(fù)雜度根據(jù)卷積和FFT算法的不同而有較大差異。下面以2種算法中耗時(shí)較長(zhǎng)的乘法計(jì)算次數(shù)為參考，分析2種算法的運(yùn)算復(fù)雜度。

2.4.1 時(shí)域卷積算法

在基于重疊保留法的時(shí)域?yàn)V波算法中，直接卷積運(yùn)算將需要2×L×M次乘法運(yùn)算；由于數(shù)據(jù)前向擴(kuò)展M-1點(diǎn)，共增加約M(M-1)次乘法運(yùn)算。所以采用基于重疊保留法的時(shí)域?yàn)V波算法需要的乘法運(yùn)算次數(shù)為：

N=2ML+M(M-1)=M(2L+M-1)。

(12)

當(dāng)數(shù)據(jù)長(zhǎng)度與濾波器系數(shù)相當(dāng)，即M≈L時(shí)：

N≈3L2。

(13)

當(dāng)數(shù)據(jù)長(zhǎng)度遠(yuǎn)大于濾波器系數(shù)，即M?L時(shí)：

N≈2ML。

(14)

對(duì)比式(13)、式(14)可知，當(dāng)濾波器階數(shù)較少時(shí)，時(shí)域方法乘法運(yùn)算量與數(shù)據(jù)長(zhǎng)度近似成線性關(guān)系。隨著濾波器階數(shù)的增加，運(yùn)算復(fù)雜度急劇增大，當(dāng)濾波器階數(shù)與數(shù)據(jù)長(zhǎng)度相當(dāng)時(shí)，乘法運(yùn)算量與數(shù)據(jù)長(zhǎng)度近似成平方關(guān)系。由此可見(jiàn)，時(shí)域卷積方法在低階條件下將具有更大的優(yōu)勢(shì)。

2.4.2 頻域卷積算法

當(dāng)采用頻域?yàn)V波算法時(shí)，需要首先對(duì)數(shù)據(jù)和濾波器系數(shù)進(jìn)行擴(kuò)展，使得2路數(shù)據(jù)長(zhǎng)度均達(dá)到L+M-1點(diǎn)，然后對(duì)2路擴(kuò)展數(shù)據(jù)做FFT運(yùn)算，則FFT運(yùn)算實(shí)際乘法運(yùn)算次數(shù)為：

N=4×[(L+M-1)/2×lb(L+M-1)+(L+M-

1)/2×lb(L+M-1)]=

4×[(L+M-1)lb(L+M-1)]。

(15)

2個(gè)序列相乘需要經(jīng)過(guò)L+M-1次乘法，得到L+M-1個(gè)結(jié)果，再對(duì)L+M-1點(diǎn)數(shù)據(jù)進(jìn)行IFFT運(yùn)算，需要2×[(L+M-1)lb(L+M-1)]次乘法。由此可知，采用頻域卷積需要的乘法次數(shù)為：

N=6×(L+M-1)lb(L+M-1)+(L+M-1)=

(L+M-1)[6lb(L+M-1)+1] 。

(16)

當(dāng)數(shù)據(jù)長(zhǎng)度與濾波器系數(shù)相當(dāng)，即M≈L時(shí)：

N= 2L[6lb(2L)+1]。

(17)

當(dāng)數(shù)據(jù)長(zhǎng)度遠(yuǎn)大于濾波器系數(shù)，即M?L時(shí)：

N≈L[6lbL+1]。

(18)

2.4.3 2種算法對(duì)比

從以上分析可以看出，時(shí)域?yàn)V波方法過(guò)程簡(jiǎn)單，當(dāng)濾波器階數(shù)較少時(shí)，運(yùn)算效率較高；但受卷積運(yùn)算影響，當(dāng)濾波器階數(shù)較大時(shí)，單次運(yùn)算量過(guò)大，難以滿(mǎn)足要求。而頻域?yàn)V波方法充分利用了FFT運(yùn)算的高效結(jié)構(gòu)，能夠明顯降低乘法運(yùn)算次數(shù)，而且隨著濾波器階數(shù)的增加和處理數(shù)據(jù)點(diǎn)數(shù)的增加，這種優(yōu)勢(shì)將更加明顯；但頻域?yàn)V波方法在運(yùn)算過(guò)程中需要對(duì)參與運(yùn)算的濾波器系數(shù)和輸入數(shù)據(jù)進(jìn)行擴(kuò)充，并在擴(kuò)充的基礎(chǔ)上開(kāi)展2次FFT運(yùn)算、一次逐點(diǎn)乘法運(yùn)算和一次IFFT運(yùn)算，給運(yùn)算過(guò)程帶來(lái)額外的負(fù)擔(dān)，而且這種負(fù)擔(dān)在濾波器階數(shù)較低時(shí)將更為顯著。

為了對(duì)時(shí)域和頻域?yàn)V波方法在不同數(shù)據(jù)條件下的性能進(jìn)行定量分析，本文通過(guò)實(shí)驗(yàn)對(duì)2種方法的性能進(jìn)行了驗(yàn)證。輸入仿真數(shù)據(jù)長(zhǎng)度為L(zhǎng)為1 024點(diǎn)，濾波器系數(shù)序列點(diǎn)數(shù)M分別選擇4，8，16，32，64，128，256，512，1 024。如此參數(shù)設(shè)置將確保仿真過(guò)程能夠涵蓋從M?L到M≈L的整個(gè)范圍，從而對(duì)2種算法的性能進(jìn)行全面的分析。

仿真結(jié)果如圖4所示，設(shè)N為運(yùn)算量，M為濾波器系數(shù)個(gè)數(shù)，L為輸入數(shù)據(jù)總長(zhǎng)度。圖中藍(lán)色曲線表示時(shí)域卷積算法在不同濾波器系數(shù)M條件下運(yùn)算量N隨著輸入數(shù)據(jù)總點(diǎn)數(shù)L的變化關(guān)系。綠色曲線表示頻域卷積算法在不同M條件下N隨著L的變化關(guān)系。圖中ConvM和fftM分別表示濾波器系數(shù)為M時(shí)，時(shí)域卷積和頻域?yàn)V波算法的運(yùn)算量。

圖4 時(shí)域卷積和頻域卷積性能對(duì)比Fig.4 Performance comparison between time-domain convolution and frequency-domain convolution

從仿真結(jié)果可知：

① 隨著濾波器階數(shù)的增加，時(shí)域?yàn)V波方法的運(yùn)算量隨著處理數(shù)據(jù)點(diǎn)數(shù)的增加而迅速增大，與之相反，頻域?yàn)V波方法在M值增加的過(guò)程中運(yùn)算量并未出現(xiàn)劇烈變化。由此可知，在濾波器階數(shù)較大的條件下，頻域?yàn)V波算法具有較為明顯的優(yōu)勢(shì)，而且這種優(yōu)勢(shì)隨著M值的增加而更加明顯，這與上文分析一致。

② 在濾波器階數(shù)與數(shù)據(jù)長(zhǎng)度相當(dāng)時(shí)，頻域卷積算法得優(yōu)勢(shì)最為明顯。

③ 在M值與L值相比明顯較小時(shí)，頻域?yàn)V波方法的不足隨著M值的減小而越發(fā)顯著。在本文仿真條件下，當(dāng)M值等于或小于16時(shí)，頻域?yàn)V波方法在不同數(shù)據(jù)長(zhǎng)度時(shí)的運(yùn)算量和運(yùn)算量增加速率均明顯高于時(shí)域算法。

④ 當(dāng)L值在1 024量級(jí)，M=32時(shí)，時(shí)域和頻域2種卷積方法的計(jì)算量相當(dāng)，且頻域算法計(jì)算量增加速率更快。因此，在一次處理的輸入數(shù)據(jù)較大而M不大于32時(shí)，時(shí)域算法的運(yùn)算量將低于頻域算法。

綜合以上仿真結(jié)果分析，在長(zhǎng)序列實(shí)時(shí)濾波運(yùn)算中可得出如下結(jié)論：

① FIR濾波器階數(shù)是影響卷積運(yùn)算效率的主要因素，數(shù)據(jù)分段方法在一定條件下可以影響運(yùn)算效率；

② 更高的濾波器階數(shù)只有與之相當(dāng)?shù)臄?shù)據(jù)長(zhǎng)度時(shí)才能體現(xiàn)頻域?yàn)V波的高效率優(yōu)勢(shì)，即當(dāng)濾波器階數(shù)較高時(shí)，數(shù)據(jù)段的長(zhǎng)度最佳選擇是與濾波器階數(shù)相當(dāng)；

③ 當(dāng)濾波器階數(shù)較小，小于或等于16時(shí)，無(wú)論分段長(zhǎng)度多長(zhǎng)，時(shí)域?yàn)V波效率均優(yōu)于頻域?yàn)V波效率。且在此條件下，為了減少分段處理時(shí)額外的數(shù)據(jù)拷貝消耗，更長(zhǎng)的數(shù)據(jù)分段長(zhǎng)度將更有利于濾波運(yùn)算效率的提升。

在信道化運(yùn)算過(guò)程中，各路濾波器為原型濾波器多相分解之后的序列，其長(zhǎng)度將是M/D，其中D為多相分支的路數(shù)。在寬帶數(shù)據(jù)處理中，為了最大限度降低信號(hào)速率，一般會(huì)選擇較大的D，這將導(dǎo)致每一路分支濾波器階數(shù)的降低。若原型濾波器長(zhǎng)度為256，分路數(shù)D為16，則每一路分支濾波器長(zhǎng)度即為16，在此條件下選擇時(shí)域?yàn)V波算法將是最優(yōu)選擇，而每次處理的信號(hào)長(zhǎng)度盡可能大將有助于提升數(shù)據(jù)處理效率。

3 實(shí)驗(yàn)驗(yàn)證

基帶轉(zhuǎn)換器是VLBI系統(tǒng)數(shù)字信號(hào)采集和處理的核心裝備，其主要功能是為后端處理系統(tǒng)提供多通道的基帶信號(hào)。由于干涉測(cè)量系統(tǒng)需要利用帶寬綜合技術(shù)對(duì)大帶寬內(nèi)的多個(gè)子帶進(jìn)行處理，因此要求基帶轉(zhuǎn)換器能夠?qū)拵盘?hào)進(jìn)行信道化處理，為后端處理系統(tǒng)輸出所需的子帶信號(hào)。多相信道化算法由于其高效優(yōu)勢(shì)在基帶轉(zhuǎn)換器的實(shí)現(xiàn)過(guò)程中得到了廣泛應(yīng)用，目前國(guó)際主流基帶轉(zhuǎn)換器均采用該算法進(jìn)行寬帶信號(hào)的實(shí)時(shí)處理。

為了驗(yàn)證以上分析結(jié)果的正確性，按照式(8)的分析結(jié)果，設(shè)計(jì)了基于GPU的VLBI多相信道化器。設(shè)每個(gè)分支濾波模塊輸入信號(hào)x(n)長(zhǎng)度為L(zhǎng)，對(duì)應(yīng)的FIR分支濾波器系數(shù)序列h(n)長(zhǎng)度為M，則基于CUDA的分支濾波算法流程如圖5所示。

圖5 分支濾波算法流程Fig.5 Block diagram of branch filtering algorithm flow

圖中分支濾波模塊是整個(gè)信道化模塊中運(yùn)算量最大的部分，每一個(gè)分支濾波器是一個(gè)獨(dú)立的濾波運(yùn)算單元，完成輸入數(shù)據(jù)和濾波器系數(shù)的卷積運(yùn)算。在多項(xiàng)濾波實(shí)現(xiàn)中，各通道輸入數(shù)據(jù)和濾波器系數(shù)均已經(jīng)通過(guò)多相分解實(shí)現(xiàn)了并行化，將各通道數(shù)據(jù)與運(yùn)算流程解耦合，使得各通道數(shù)據(jù)在處理過(guò)程中相互獨(dú)立。而在各通道內(nèi)部，輸入數(shù)據(jù)與分支濾波器通過(guò)卷積運(yùn)算實(shí)現(xiàn)濾波操作。

由圖5可以看出，輸入數(shù)據(jù)被分成并行的D路，原型濾波器也做相應(yīng)的多相分解，使得分支路數(shù)也為D，然后根據(jù)分支路數(shù)啟動(dòng)D路并行的CUDA stream。每一路數(shù)據(jù)的運(yùn)算過(guò)程發(fā)生在各自的CUDA stream內(nèi)，這種設(shè)計(jì)一方面充分利用了各路數(shù)據(jù)的并行性，使得各路數(shù)據(jù)處理流程在stream之間完全獨(dú)立，消除了因串行等待造成的時(shí)間延遲；另一方面采用這種方法有效合并了各路數(shù)據(jù)的內(nèi)存訪問(wèn)，解決了因內(nèi)存讀取不連續(xù)造成的效率降低下的問(wèn)題。

由式(9)可知，在每一個(gè)CUDA stream內(nèi)，實(shí)際上進(jìn)行的是M點(diǎn)濾波器系數(shù)和L點(diǎn)輸入數(shù)據(jù)的線性卷積運(yùn)算，而對(duì)于一個(gè)輸出點(diǎn)而言，卷積運(yùn)算本質(zhì)上是2個(gè)M點(diǎn)長(zhǎng)度序列的逐點(diǎn)對(duì)應(yīng)乘加。在圖5所示的算法流程中，每個(gè)線程負(fù)責(zé)完成M次乘加運(yùn)算并最終實(shí)現(xiàn)一個(gè)卷積結(jié)果的輸出，所有L個(gè)線程并行啟動(dòng)即可完成所有輸出點(diǎn)的運(yùn)算。為了降低因運(yùn)算過(guò)程中對(duì)濾波器系數(shù)反復(fù)讀取造成的運(yùn)算效率損失，利用線程塊的共享內(nèi)存存儲(chǔ)參與運(yùn)算的數(shù)據(jù)和系數(shù)。

本文VLBI基帶轉(zhuǎn)換器基于CUDA10.2和VS2015開(kāi)發(fā)，GPU采用NVIDIA Tesla V100，信道化軟硬件開(kāi)發(fā)設(shè)計(jì)采用圖2、圖3所示多相信道化結(jié)構(gòu)。輸入數(shù)據(jù)在信道化入口即按通道數(shù)進(jìn)行串并轉(zhuǎn)換，劃分為均勻的多路，然后按照通道數(shù)啟動(dòng)CUDA stream對(duì)每一之路數(shù)據(jù)異步并行處理。在時(shí)域卷積模式下每一個(gè)流內(nèi)按照?qǐng)D5所示流程進(jìn)行分支濾波操作，每一個(gè)流內(nèi)數(shù)據(jù)處理順序進(jìn)行，各流之間數(shù)據(jù)傳輸和核函數(shù)的執(zhí)行異步并發(fā)。在各流完成各自操作后，對(duì)各通道數(shù)據(jù)按通道數(shù)進(jìn)行D點(diǎn)FFT，得到并行的多路基帶數(shù)據(jù)輸出。輸出數(shù)據(jù)按要求將寬帶信號(hào)均勻劃分為D路。

輸入數(shù)據(jù)采樣率為1 024 Ms/s，信道化輸出16路復(fù)信號(hào)，每一路帶寬64 MHz。取0.156 25 s原始數(shù)據(jù)進(jìn)行測(cè)試，即總共數(shù)據(jù)點(diǎn)數(shù)為1 677 216。多相濾波原型濾波器階數(shù)255，則每路分支濾波器系數(shù)點(diǎn)數(shù)為16。分別按照每段數(shù)據(jù)長(zhǎng)度1 024，4 096，65 536，262 144點(diǎn)對(duì)原始數(shù)據(jù)進(jìn)行分段處理。分別統(tǒng)計(jì)每種分段長(zhǎng)度條件下時(shí)域卷積算法和頻域?yàn)V波算法條件下信道化所用時(shí)長(zhǎng)，結(jié)果如圖6所示。

圖6 時(shí)域卷積和頻域卷積性能對(duì)比Fig.6 Performance comparison of time domain convolution and frequency domain convolution

由圖6可以看出，當(dāng)濾波器階數(shù)為15時(shí)，各種數(shù)據(jù)分段長(zhǎng)度下，時(shí)域?yàn)V波性能均優(yōu)于頻域?yàn)V波性能，且隨著數(shù)據(jù)分段長(zhǎng)度的增加，2種濾波方法性能均明顯增強(qiáng)，這與本文的預(yù)測(cè)完全一致。但是，受GPU片上資源限制，單次處理的數(shù)據(jù)長(zhǎng)度不可能無(wú)限制增加，在數(shù)據(jù)長(zhǎng)度超過(guò)65 536點(diǎn)時(shí)，算法性能增加速度變慢甚至出現(xiàn)下降。因此，在本文所示信道化運(yùn)算條件下，以單通道每次處理65 536點(diǎn)數(shù)據(jù)為最佳設(shè)置。另外,從圖6放大部分可看出，2種算法中每一段數(shù)據(jù)的處理耗時(shí)均隨著處理數(shù)據(jù)長(zhǎng)度的增加而近似線性增加，且頻域算法的耗時(shí)增加速率明顯高于時(shí)域算法，這與本文預(yù)測(cè)一致。

4 結(jié)束語(yǔ)

以干涉測(cè)量多相信道化基帶轉(zhuǎn)換器為應(yīng)用背景，以基于GPU的通用計(jì)算平臺(tái)為應(yīng)用目標(biāo)，提出一種基于CUDA的多相信道化實(shí)現(xiàn)方法，為了提高FIR濾波環(huán)節(jié)的運(yùn)算效率，針對(duì)信道化算法中運(yùn)算量最大的分支濾波模塊，利用重疊保留法設(shè)計(jì)了基于CUDA的頻域?yàn)V波算法和時(shí)域?yàn)V波算法。然后，利用仿真平臺(tái)對(duì)2種算法的運(yùn)算復(fù)雜度和適用條件進(jìn)行了分析，結(jié)果顯示：

① 時(shí)域?yàn)V波方法簡(jiǎn)單高效，但受卷積運(yùn)算復(fù)雜度影響，單次運(yùn)算量大，在濾波器階數(shù)較低(≤16)且數(shù)據(jù)速率較高時(shí)適用；

② 頻域?yàn)V波方法流程復(fù)雜，但FFT算法的應(yīng)用使得運(yùn)算量大大降低，適用于濾波器階數(shù)較大(≥32)的場(chǎng)合。

利用以上分析結(jié)果設(shè)計(jì)了基于GPU的16通道多相信道化算法實(shí)現(xiàn)結(jié)構(gòu)，并對(duì)該結(jié)構(gòu)進(jìn)行了測(cè)試，結(jié)果顯示在各分支濾波器系數(shù)長(zhǎng)度為16條件下，隨著處理數(shù)據(jù)點(diǎn)數(shù)的增加，采用多通道并行時(shí)域卷積算法能夠達(dá)到更高的運(yùn)算效率，與分析結(jié)果一致。未來(lái)，該分析結(jié)果有望為干涉測(cè)量系帶轉(zhuǎn)換器的高效實(shí)現(xiàn)提供有效的技術(shù)支持。