基于嵌入式移動(dòng)GPU的離散傅里葉變換并行優(yōu)化※

曾寶國(guó)，楊斌

(1.成都工業(yè)職業(yè)技術(shù)學(xué)院，成都 610213； 2.西南交通大學(xué))

?

基于嵌入式移動(dòng)GPU的離散傅里葉變換并行優(yōu)化※

曾寶國(guó)1，楊斌2

(1.成都工業(yè)職業(yè)技術(shù)學(xué)院，成都 610213； 2.西南交通大學(xué))

摘要:GPGPU能夠針對(duì)計(jì)算密集型的計(jì)算問(wèn)題進(jìn)行大規(guī)模的并行加速，為DFT在嵌入式平臺(tái)上的高效實(shí)現(xiàn)提供了一種新的方式?；贛ali-T604嵌入式GPU實(shí)現(xiàn)了針對(duì)DFT和FFT的并行加速方案，并進(jìn)行了實(shí)際測(cè)試。實(shí)驗(yàn)結(jié)果證明，所設(shè)計(jì)的并行方案能夠在ARM嵌入式平臺(tái)上有效加速DFT和FFT，可大大提升移動(dòng)設(shè)備進(jìn)行數(shù)字信號(hào)處理的實(shí)時(shí)性。

關(guān)鍵詞:DFT；FFT；GPGPU；Mali-T604 GPU；數(shù)字信號(hào)處理；ARM嵌入式系統(tǒng)

引言

GPGPU(General purpose GPU)技術(shù)近年來(lái)在嵌入式領(lǐng)域廣泛應(yīng)用，移動(dòng)GPU對(duì)于輸入量龐大的計(jì)算密集型、數(shù)據(jù)可并行化的通用計(jì)算問(wèn)題有顯著的加速能力[1]。DFT以及FFT的運(yùn)算復(fù)雜度極高，嵌入式平臺(tái)上計(jì)算能力有限的CPU難以對(duì)其進(jìn)行快速處理。

在實(shí)時(shí)信號(hào)分析場(chǎng)景下，需要高效的計(jì)算方案，DFT和FFT的數(shù)學(xué)模型適合使用GPU對(duì)其進(jìn)行并行加速。本文基于Mali-T604 GPU對(duì)DFT和FFT的并行計(jì)算方案進(jìn)行設(shè)計(jì)，提供了實(shí)現(xiàn)方法和實(shí)際測(cè)試結(jié)果，為使用GPGPU技術(shù)在嵌入式平臺(tái)上進(jìn)行數(shù)字信號(hào)處理的研究人員提供了參考和借鑒。

1離散傅里葉變換并行化解析

離散傅里葉變換在頻譜分析、數(shù)字通信、圖像處理、遙感遙測(cè)等領(lǐng)域有著廣泛的應(yīng)用，有多種算法以不同的時(shí)間復(fù)雜度對(duì)其進(jìn)行實(shí)現(xiàn)，常規(guī)的DFT方式和Cooley-Tukey快速變換方式是最常見(jiàn)的兩種離散傅里葉變換實(shí)現(xiàn)方式。

1.1常規(guī)DFT數(shù)學(xué)模型分析

數(shù)字信號(hào)處理中應(yīng)用的離散傅里葉變換通常用于將時(shí)域采樣信號(hào)轉(zhuǎn)換至頻域進(jìn)行分析，其輸入數(shù)據(jù)通常是實(shí)數(shù)，對(duì)于長(zhǎng)度為N的實(shí)數(shù)輸入信號(hào)序列x[0DK∶N-1]，其一維DFT變換公式為：

(1)

從式(1)可以看出，從實(shí)數(shù)序列x[0∶N-1]到復(fù)數(shù)序列X[0∶N-1]的變換過(guò)程實(shí)際上可轉(zhuǎn)換為矩陣乘法的形式：

(2)

(3)

可見(jiàn)對(duì)于輸出序列X[0∶N-1]的每個(gè)元素而言，其變換過(guò)程是相互獨(dú)立的，其中涉及的是大量的一維向量?jī)?nèi)積計(jì)算，對(duì)于數(shù)字信號(hào)而言，輸入序列x[0∶N-1]為實(shí)數(shù)浮點(diǎn)序列，式(2)和式(3)涵蓋的是多次浮點(diǎn)乘法和三角函數(shù)計(jì)算，嵌入式的CPU(如ARM處理器)對(duì)于這樣的計(jì)算,速度是比較慢的，如果每個(gè)元素的變換過(guò)程能夠獨(dú)立分解到不同的任務(wù)中并行執(zhí)行,且浮點(diǎn)運(yùn)算能夠快速完成，則常規(guī)的DFT就能夠快速完成，而MaliGPU正好就是符合這樣計(jì)算特征的協(xié)處理器。

Cooley-Tukey在1965年提出了快速傅里葉變換算法，該算法高效簡(jiǎn)單，在數(shù)字信號(hào)處理領(lǐng)域經(jīng)久不衰，該算法利用了離散傅里葉變換的疊加性質(zhì)、移位性質(zhì)和延展性質(zhì),簡(jiǎn)化了DFT變換的運(yùn)算量，能夠得出和常規(guī)DFT變換相同的運(yùn)算結(jié)果。

Cooley-TukeyFFT算法可處理長(zhǎng)度為2次冪的任意序列的變換，其核心思想是歸并，不斷利用短序列的變換結(jié)果歸并產(chǎn)生長(zhǎng)序列的變換結(jié)果。在計(jì)算變換結(jié)果X[0∶N-1]的時(shí)候，首先分別計(jì)算子序列x[0∶N/2-1]和x[N/2,N-1]的變換結(jié)果X[0∶N/2-1]和X[N/2∶N-1]，再歸并產(chǎn)生X[0∶N-1]。根據(jù)原始序列x[0∶N-1]構(gòu)建序列x1[0∶N-1]={x(0),0,x(1),0,…,x(N/2-1), 0}和序列x2[0∶N-1]= {0,x(N/2),0,x(N/2+1),…,0,x(N-1)}，根據(jù)DFT的延展性質(zhì)和移位性質(zhì)，有下式成立：

(4)

(5)

(6)

根據(jù)DFT變換的疊加性質(zhì)，有下式成立：

(7)

4個(gè)元素的 FFT相對(duì)簡(jiǎn)單，算法從4個(gè)元素的變換開(kāi)始，將原始序列的每4個(gè)元素分解為一組，每組元素首先用蝶形運(yùn)算進(jìn)行DFT變換，然后利用式(4)～(7)的算法不斷歸并，得到更大長(zhǎng)度的FFT變換結(jié)果，直至歸并長(zhǎng)度等于原始序列長(zhǎng)度N為止。

Cooley-Tukey FFT算法的時(shí)間復(fù)雜度為O(N×log2(N))，相對(duì)常規(guī)的DFT算法的O(N2)復(fù)雜度而言，運(yùn)算量有所下降，F(xiàn)FT算法中涉及的4元素分組的蝶形運(yùn)算可以獨(dú)立進(jìn)行，互不干擾，后期每個(gè)子段的歸并過(guò)程也是互不依賴的，涉及的也多為浮點(diǎn)乘法運(yùn)算，適合使用Mali GPU進(jìn)行性能優(yōu)化。

2Mali GPU并行計(jì)算模型簡(jiǎn)介

Mali-T600系列的GPU可以支持OpenCL 1.1 Full Profile標(biāo)準(zhǔn)，OpenCL是真正意義上的跨平臺(tái)異構(gòu)并行框架，能夠真正挖掘出Mali GPU的并行計(jì)算特性。在OpenCL框架下，設(shè)計(jì)者可將預(yù)定數(shù)目的計(jì)算任務(wù)下載至Mali GPU端，以多線程形式實(shí)現(xiàn)全局的任務(wù)并行，其運(yùn)行過(guò)程略——編者注。

OpenCL通過(guò)主程序來(lái)定義上下文，并管理內(nèi)核程序在GPU設(shè)備的執(zhí)行[2]，應(yīng)用程序通過(guò)host 提交命令, 驅(qū)動(dòng)設(shè)備上的PE 執(zhí)行內(nèi)核程序[3]。GPU可以同時(shí)處理成百上千的線程，大量晶體管用于ALU[4]，每個(gè)Mali-T604 GPU的著色器核心最多可以同時(shí)容納256個(gè)線程[5]。在單個(gè)線程內(nèi)部，可以利用Mali GPU內(nèi)置的128位寄存器和ALU實(shí)現(xiàn)向量數(shù)據(jù)的SIMD(Single Instruction Multiple Data)局部并行加速，Mali-T604 GPU的4個(gè)著色器核心以及其中的2條向量運(yùn)算流水線和1條向量讀寫(xiě)流水線能夠高效完成128位向量的讀寫(xiě)和運(yùn)算。

3并行DFT以及并行FFT變換實(shí)現(xiàn)

傳統(tǒng)的DFT以及Cooley-Tukey FFT算法均適合使用Mali GPU進(jìn)行并行優(yōu)化，下文分別針對(duì)這兩種算法介紹并行化的方法。

3.1常規(guī)DFT算法并行化實(shí)現(xiàn)

常規(guī)的DFT算法中，變換結(jié)果X[0∶N-1]中的每個(gè)元素都是由一維向量的內(nèi)積得到的，因此每一個(gè)元素的計(jì)算過(guò)程可以分解到單個(gè)線程中完成，對(duì)于實(shí)數(shù)序列的DFT變換而言，X(0)和X(N/2)的虛數(shù)部分始終為0，對(duì)于其他的元素而言，X[k]和X[N-k]具有共軛復(fù)數(shù)關(guān)系，所以真正的計(jì)算中只需要計(jì)算X[0∶N/2]。因此GPU端總共分配N/2+1個(gè)線程，每個(gè)元素完成一次長(zhǎng)度為N的一維向量的內(nèi)積，得出X[0∶N/2]中的一個(gè)標(biāo)量計(jì)算結(jié)果，圖1顯示了單個(gè)線程的運(yùn)行流程。

圖1　并行DFT變換運(yùn)算核流程圖

OpenCL實(shí)現(xiàn)的并行DFT內(nèi)核源碼如下所示：

__kernel void dft(__global float *dft_src,

//變換之前的原始數(shù)據(jù)緩存，長(zhǎng)度為N

__global float *dft_dst, //變換結(jié)果的緩存，長(zhǎng)度為N

const int N //dft變換的長(zhǎng)度){

……

thread_id = get_global_id(0); two_pai_kN = 2.0f*PI*thread_id/N;

sum_real = sum_imag = 0.0f;

for(i=0; i<(N>>2); i++){

src_vec = vload4(i, dft_src); //讀取輸入向量

angle = (float4)(two_pai_kN) * (float4)( (i<<2)+0, (i<<2)+1, (i<<2)+2, (i<<2)+3);

real = cos(angle); imag = sin(angle);

sum_real += dot(real, src_vec); sum_imag -= dot(imag, src_vec);

}

……

}

3.2Cooley-Tukey FFT算法并行優(yōu)化實(shí)現(xiàn)

Cooley-Tukey FFT算法分為前后兩部分，前半部分每4個(gè)元素一組,利用蝶形運(yùn)算進(jìn)行FFT變換，4個(gè)元素進(jìn)行蝶形運(yùn)算的方式如圖2所示。

圖2　4個(gè)元素的FFT蝶形運(yùn)算

在FFT算法中，一開(kāi)始每4個(gè)元素的FFT變換是相互獨(dú)立的，因此分配N/4個(gè)線程對(duì)每一組數(shù)據(jù)進(jìn)行圖2所示的蝶形運(yùn)算，圖2中的X(0)～X(3)需要從原始數(shù)據(jù)緩存中讀取，其讀取位置和分組的組號(hào)相關(guān)，設(shè)線程i(i=0,1,2,…,N/4-1)處理第i組數(shù)據(jù)中的4個(gè)元素，則線程i讀取的4個(gè)元素的下標(biāo)應(yīng)該分別是4i、4i+1、4i+2、4i+3四個(gè)數(shù)字進(jìn)行位翻轉(zhuǎn)之后的結(jié)果。設(shè)log2(N)=b，則下標(biāo)index位翻轉(zhuǎn)的結(jié)果應(yīng)該是其二進(jìn)制第(b-1)位和第0位交換、第(b-2)位和第1位交換……直至較高數(shù)位和較低數(shù)位全部交換一遍之后的結(jié)果，能夠用Mali GPU進(jìn)行位運(yùn)算在每個(gè)線程中并行執(zhí)行，這樣能夠保證最終歸并后的結(jié)果是順序排列的。圖2中蝶形運(yùn)算的輸入量都是實(shí)數(shù)，其計(jì)算過(guò)程可以拆分成實(shí)數(shù)部分和虛數(shù)部分分別進(jìn)行，并且都能夠用Mali GPU中的4通道浮點(diǎn)向量高效完成，代碼如下：

__kernel void fft4(__global float *fft_src,

//變換之前的原始數(shù)據(jù)緩存，長(zhǎng)度是N

__global float *fft_dst_real,

//變換結(jié)果的緩存，長(zhǎng)度為N(保存實(shí)部)

__global float *fft_dst_imag,

//變換結(jié)果的緩存，長(zhǎng)度為N(保存虛部)

const int bits

//bits=log2(N),用于位翻轉(zhuǎn)計(jì)算

){

int id = get_global_id(0);

……//位翻轉(zhuǎn)部分省略

X_real = (float4)(x0);X_real += (float4)(x1, 0-x1, x1, 0-x1);

X_real += (float4)(x2, 0, 0-x2, 0);X_real += (float4)(x3, 0, 0-x3, 0);

X_imag = (float4)(0.0f, x3-x2, 0.0f, x2-x3);

vstore4(X_real, id, fft_dst_real);vstore4(X_imag, id, fft_dst_imag); //回存結(jié)果

}

Cooley-Tukey FFT算法的后半部分需要對(duì)前半部分的計(jì)算結(jié)果進(jìn)行歸并，對(duì)于長(zhǎng)度為N(N為2的冪)的FFT變換，在進(jìn)行了4元素分組的FFT蝶形運(yùn)算后，還需進(jìn)行l(wèi)og2N-2次歸并，第i次(i從1開(kāi)始)歸并過(guò)程的歸并長(zhǎng)度為2i+2，即將相鄰的兩個(gè)長(zhǎng)度為2i+1的FFT變換序列歸并為長(zhǎng)度為2i+2的序列。對(duì)式(4)～(7)進(jìn)行分析可知，雖然每次歸并的長(zhǎng)度翻倍后歸并的分組會(huì)減少，但總計(jì)算量是不變的，因此，將長(zhǎng)度為N的浮點(diǎn)復(fù)數(shù)序列每4個(gè)相鄰元素分成一組，則不論在第幾次歸并過(guò)程中，必然能夠讓每組數(shù)據(jù)兩兩配對(duì)，且分別屬于歸并運(yùn)算的兩個(gè)輸入子段，這樣每個(gè)線程可分別讀取這兩個(gè)長(zhǎng)度為4的子段，進(jìn)行式(4)～(7)的歸并計(jì)算，并可將實(shí)部和虛部的計(jì)算拆開(kāi)，4個(gè)浮點(diǎn)數(shù)剛好可以湊齊128位,用Mali GPU的向量運(yùn)算進(jìn)行處理，綜上所述，每一次歸并過(guò)程可以分配N/8個(gè)線程，每個(gè)線程進(jìn)行長(zhǎng)度為4的子段的歸并計(jì)算。

設(shè)Mali GPU端編號(hào)為i(i=0,1,2,…,N/8)的線程在參與歸并長(zhǎng)度為merge_len的歸并過(guò)程中，上一次歸并的結(jié)果為X_old[0DK∶N-1]，則線程i應(yīng)首先讀取X_old中的第一個(gè)長(zhǎng)度為4的子段vec0，其偏移地址為：

offset0=(i/(merge_len>>3))×merge_len

(8)

后一個(gè)子段vec1和vec0間的地址間隔是由歸并長(zhǎng)度決定的：

(9)

根據(jù)前文分析，由vec0和vec1兩個(gè)向量可以計(jì)算出新序列X_new[0DK∶N-1]中的兩個(gè)長(zhǎng)度為4的新子段X_new[offset0∶offset0+3]和X_new[of]fset1∶offset1+3]，下面是OpenCL歸并內(nèi)核的核心代碼略——編者注。

由于Cooley-Tukey FFT算法中每次歸并都依賴于上一次歸并的結(jié)果，因此在ARM CPU主機(jī)端需要將歸并內(nèi)核順序log2N次加入Mali GPU在OpenCL框架下的命令隊(duì)列中，并且每次傳入的歸并長(zhǎng)度參數(shù)merge_len翻倍。

4優(yōu)化效果測(cè)試

筆者在Arndale Board開(kāi)發(fā)板(核心為ARM Cortex-A15雙核CPU+Mali-T604 GPU的Exynos5250 SoC)上和嵌入式Linux操作系統(tǒng)上，對(duì)所述的并行DFT和并行FFT優(yōu)化方案進(jìn)行了測(cè)試，對(duì)比ARM Cortex-A15 CPU進(jìn)行的串行方案和Mali GPU進(jìn)行的并行方案，以不同長(zhǎng)度的隨機(jī)浮點(diǎn)實(shí)數(shù)序列進(jìn)行離散傅里葉變換效率測(cè)試，得到的測(cè)試數(shù)據(jù)如表1所列。

表1　離散傅里葉變換并行優(yōu)化效果對(duì)比

從表1的測(cè)試數(shù)據(jù)可以看出，在離散傅里葉變換的變換長(zhǎng)度較短的時(shí)候，Mali GPU端的運(yùn)算線程數(shù)量不足，其并行運(yùn)算能力未能發(fā)揮出來(lái)，故并行方案的計(jì)算效率不及串行方案。隨著變換長(zhǎng)度的增加，Mali GPU端的并發(fā)線程數(shù)量上升，并行方案的計(jì)算效率迅速上升，不論是常規(guī)DFT算法還是Cooley-Tukey FFT算法的并行方案,效率都遠(yuǎn)超串行方案，加速比達(dá)到了百倍以上，證明所設(shè)計(jì)的并行方案的加速效果穩(wěn)定而有效。對(duì)比并行的DFT方案和并行的FFT方案，在變換長(zhǎng)度較低的時(shí)候，DFT方案的效率更高，雖然FFT方案的總計(jì)算量低于常規(guī)的DFT算法，但是在歸并過(guò)程需要多次將歸并內(nèi)核入隊(duì)，帶來(lái)了內(nèi)核運(yùn)行之間的同步開(kāi)銷。在并發(fā)線程不足的情況下，并行化帶來(lái)的收益未能彌補(bǔ)同步開(kāi)銷，所以并行FFT方案的性能提升不高，在變換長(zhǎng)度較大的場(chǎng)景下，Mali GPU端的運(yùn)算并行度增加，同步開(kāi)銷成為影響效率的次要因素，F(xiàn)FT算法的低運(yùn)算量特征體現(xiàn)出其優(yōu)勢(shì)，使得后期并行FFT方案的效率遠(yuǎn)超并行DFT方案。在實(shí)際工程中，應(yīng)針對(duì)不同的嵌入式GPU硬件結(jié)構(gòu)和變換長(zhǎng)度，綜合選取兩者中效率更高的方案投入具體的應(yīng)用場(chǎng)景中。

結(jié)語(yǔ)

本文基于Mali-T604 GPU設(shè)計(jì)了針對(duì)常規(guī)DFT算法和Cooley-Tukey FFT算法的并行優(yōu)化方案，并在嵌入式Linux操作系統(tǒng)和OpenCL框架下進(jìn)行了實(shí)現(xiàn)。實(shí)際測(cè)試效果表明，該優(yōu)化方案效果明顯，嵌入式GPU的新興GPGPU技術(shù)在數(shù)字信號(hào)處理領(lǐng)域有著廣闊的應(yīng)用前景。

參考文獻(xiàn)

[1] 龔若皓,楊斌.基于移動(dòng)多核GPU的并行二維DCT變換實(shí)現(xiàn)方法[J] .成都信息工程學(xué)院學(xué)報(bào),2015,30(1):22-26.

[2] 向陽(yáng)霞,張惠民,王自強(qiáng).面向OpenCL模型的DCT并行化[J] .電腦知識(shí)與技術(shù),2013,9(26):6007-6011.

[3] 陳鋼,吳百鋒.面向OpenCL模型的GPU性能優(yōu)化[J] .計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2011,23(4):571-580.

[4] Owens J D,Houston M,Luebke D,et al.GPUComputing[J] .Proceedings of the IEEE,2008,96(5):879-899.

[5] ARM Company.ARM Mali-T600 Series GPU OpenCL Developer Guide Version 2.0，2012.

曾寶國(guó)(副教授)，主要研究方向?yàn)榍度胧綉?yīng)用開(kāi)發(fā)；楊斌(教授)，主要研究方向?yàn)榍度胧较到y(tǒng)應(yīng)用及異構(gòu)并行計(jì)算。

(責(zé)任編輯：薛士然收修改稿日期：2015-08-27)

Parallelization of DFT Based on Embedded Mobile GPU※

Zeng Baoguo1,Yang Bin2

(1.College of Chengdu Industrial Vocational Technical,Chengdu 610213;2.Southwest Jiaotong University)

Abstract：GPGPU can provide efficient parallel computing solution for the complex compute-intensive computing problem,which is a new way of the efficient implementation of DFT in the embedded platform.In the paper,the parallelization solution of DFT and FFT based on Mali-T604 GPU is proposed.The results of experiment show that the parallel scheme can effectively accelerate DFT and FFT on ARM embedded platform,which can greatly improve the real-time performance of digital signal processing.

Key words:DFT;FFT;GPGPU;Mali-T604 GPU;digital signal processing;ARM Embedded System

中圖分類號(hào):TP311

文獻(xiàn)標(biāo)識(shí)碼:A