賈 格 彭先蓉 左顥睿

1(中國科學(xué)院光電技術(shù)研究所 四川 成都 610209)2(中國科學(xué)院大學(xué) 北京 100039)

基于OpenCL的FFT算法研究

賈 格1,2彭先蓉1左顥睿1

1(中國科學(xué)院光電技術(shù)研究所 四川 成都 610209)2(中國科學(xué)院大學(xué) 北京 100039)

快速福利葉變換在圖像處理領(lǐng)域，尤其是在圖像復(fù)原算法中作為常用的計算工具，將時域計算轉(zhuǎn)變?yōu)轭l域計算，在工程應(yīng)用中有著非常重要的意義。采取多線程分塊以及并行的映射方法，可以使FFT算法最大程度并行。針對OpenCL的存儲層次特點和算法層次的優(yōu)化,在AMD GPU平臺上取得了明顯的加速效果。優(yōu)化后的算法性能比具有相同處理能力的CPU平臺提高了7倍，比具有相同處理能力的CUDA提高了4倍。

傅里葉變換 OpenCL GPU 并行加速

0 引 言

隨著GPU通用計算的高速發(fā)展，大量的算法被成功移植到GPU平臺，并且取得了良好的加速效果。而且GPU的應(yīng)用研究也越來越廣泛，通用計算適用領(lǐng)域也越來越大。但是已有的研究一般只限于單一的平臺，而OpenCL(開放式計算語言)是面向異構(gòu)計算平臺的通用編程框架[2,3,6]，為GPU通用計算的跨平臺移植提供了新的解決方案。由于GPU硬件結(jié)構(gòu)體系本身的差異，造成了GPU硬件平臺高效移植有很大難度，所以高效的移植不僅僅是功能的移植，高效的移植也就成為當(dāng)前可待研究的問題。

1 OpenCL簡介

OpenCL平臺模型定義了使用OpenCL設(shè)備的一個高層描述[1]。這個模型如圖1所示，OpenCL平臺通常包含一個宿主機(host)和多個OpenCl設(shè)備。宿主機同OpenCl設(shè)備程序外部的交互，主要包括程序和I/O的交互。宿主機與一個或多個OpenCL設(shè)備連接，設(shè)備就是執(zhí)行指令流(或kernel函數(shù))的地方，通常OpenCL設(shè)備被稱為計算設(shè)備。設(shè)備可以是GPU、CPU、DSP或者其他類型的處理器。圖2所示為OpenCl的抽象模型。OpenCL應(yīng)用執(zhí)行模型中使用主機端程序?qū)Χ鄠€支持OpenCL的計算設(shè)備進行統(tǒng)一的調(diào)度和管理。圖2中分別包括兩個OpenCl設(shè)備，一個為CPU作為宿主機，一個GPU作為執(zhí)行設(shè)備。然后定義了兩個命令隊列，一個是面向CPU的亂序命令隊列，在宿主機上執(zhí)行，另一個為面向GPU的有序命令隊列，在OpenCl設(shè)備上執(zhí)行。然后宿主機程序定義了一個程序?qū)ο?，將這個程序?qū)ο蠓謩e編譯后給兩個OpenCL設(shè)備提供內(nèi)核。接下來宿主機端定義程序所需的內(nèi)核對象，并將他們分別映射到內(nèi)核參數(shù)中。在OpenCL的內(nèi)存模型中，通常根據(jù)線程訪問權(quán)限的不同定義以下四種內(nèi)存空間：分別為全局內(nèi)存、局部內(nèi)存、常量內(nèi)存、和私有內(nèi)存，它們在程序中的使用方式很大程度上決定了程序的性能，一般來講，它們的大小和訪問速度都是不同的。

圖1 OpenCl平臺模型

圖2 OpenCl應(yīng)用執(zhí)行模型

2 FFT簡介

快速傅里葉變換FFT是離散傅利葉變換DFT的快速算法[7,9]，具有廣泛的應(yīng)用研究和工程價值，常被應(yīng)用于數(shù)字濾波和信號分解中。二維快速傅里葉變換(2D-FFT)是圖像處理領(lǐng)域，尤其是在圖像復(fù)原算法中作為常用的頻域計算工具，在工程應(yīng)用中有著非常重要的意義，典型的有圖像頻譜分析、卷積計算等。通常來說，在實際工程中會遇到高分辨率的圖像，但是對于這種圖像的二維傅里葉變換計算量非常大并且數(shù)據(jù)動態(tài)范圍寬，硬件設(shè)計實現(xiàn)難度大，一定程度上制約了其在工程中的應(yīng)用。通常的做法是利用二維傅里葉變換的數(shù)據(jù)可分性，圖像的二維FFT可分為兩次的一維FFT實現(xiàn)，分別是行方向和列方向。對于二維FFT的高速實現(xiàn)平臺，一般以高性能DSP、FPGA、通用GPGPU芯片為核心器件的三種常用平臺，并在系統(tǒng)設(shè)計中采用多核計算，多級流水線以及乒乓緩存等常用并行處理技術(shù)。本文利用二維傅里葉變換的共軛對稱性，實信號變換周期對稱性等特性減少存儲資源和計算量，并通過OpenCL工具在GPU平臺對算法加速。

3 FFT算法原理及其改進

3.1 利用周期對稱性減少FFT算法冗余計算

對于實信號的FFT，通常利用傅里葉變換的對稱性，可以將N點的實序列實現(xiàn)奇偶分離，變成一個N/2點的復(fù)數(shù)序列，將計算N點的實序列FFT變?yōu)橐粋€N/2點的復(fù)數(shù)序列FFT，然后利用對應(yīng)關(guān)系將N/2點的復(fù)數(shù)序列FFT輸出對應(yīng)到N點的實數(shù)序列結(jié)果輸出。這樣做將實信號FFT的計算量減少近一半。下面介紹它的基本原理：對于一個N點的實序列x(n)，我們根據(jù)奇偶序號，分別抽取奇數(shù)序列g(shù)(n)和偶數(shù)序列h(n)，分別將它們的離散傅里葉變換分別記為G(k)和H(k)。由FFT時間抽取基2算法得到如下關(guān)系式。

(1)

(2)

令y(n)=g(n)+j·h(n)，則y(n)為N/2點的復(fù)序列。y(n)的離散傅里葉變換Y(k)為：

(3)

(4)

(5)

將復(fù)序列Y(k)、G(k)、H(k)的實虛部分開，分別記為Yr(k)、Gr(k)、Hr(k)和Yi(k)、Gi(k)、Hi(k)，代入式(5)，則有：

(6)

(7)

同理可以得到：

(8)

(9)

從上面兩個式子可以看出，只需要得出構(gòu)造的N/2點復(fù)序列y(n)的傅里葉變換結(jié)果Y(k)，然后利用上述四個關(guān)系式就可以得到g(n)和h(n)的結(jié)果G(k)和H(k)；然后再代入式(1)和式(2)，就可以得到N點實序列x(n)的傅里葉變換結(jié)果X(k)。實驗中采取此種方法，可以將實信號FFT的計算量減少近一半，特別是針對于較高分別率的圖像的時候，在N值取值比較大，節(jié)約的計算量對于運算速度的提升非常明顯。

3.2 算法復(fù)雜度分析

如圖3上半部分表示常規(guī)的2D-FFT運算，經(jīng)過行變換后，再經(jīng)過列變換，即輸出結(jié)果；下半部分表示本文改進后的算法運算模塊，多了數(shù)據(jù)預(yù)處理模塊和數(shù)據(jù)恢復(fù)模塊。其中，數(shù)據(jù)預(yù)處理：將輸入分辨率為M×N的圖像f(x,y)，對于各行數(shù)據(jù)按奇偶序號分離；然后將每行的奇、偶序列組成N/2點的復(fù)數(shù)數(shù)據(jù)fc(x,y)。數(shù)據(jù)恢復(fù)：對于N/2點的復(fù)數(shù)中間數(shù)據(jù)結(jié)果Fc(u,v)，根據(jù)式(6)-式(9)即可以得到F(u,v)的數(shù)據(jù)結(jié)果。

圖3 常規(guī)算法和本文算法處理流程

本文改進方法大幅縮減了一維FFT的計算量，從N點變換轉(zhuǎn)化為N/2點變換，計算量縮減了50%，但是引入的其他模塊會帶來一定的成本增量，導(dǎo)致計算成本要小于50%的減少量。表1在計算量，包括乘加次數(shù)給出了常規(guī)FFT算法和本文改進方法的比較。

表1 算法復(fù)雜度對比

表1中，M,N代表圖像的尺寸大小,x,y代表額外等效的參數(shù)。在計算量方面，所有的計算等效為基本的乘加次數(shù)。對于無法轉(zhuǎn)換的操作，比如數(shù)據(jù)的轉(zhuǎn)存，移位開銷則間接地等效成乘加運算。進一步可以得到本文方法相對常規(guī)方法所減少的計算量表達式：

(10)

(11)

從上式可知，減少量與圖像的大小有關(guān)，圖像越大，減少量相對更大。理論上分析，忽略數(shù)據(jù)轉(zhuǎn)存以及內(nèi)存中移位的時間開銷，僅針對額外的數(shù)據(jù)計算，根據(jù)圖3中數(shù)據(jù)預(yù)處理模塊和數(shù)據(jù)恢復(fù)模塊可知，x=2，y=4。對于128×128的圖像，γ=39.286%,η=39.286%;對于512×512的圖像，γ=41.7%,η=41.7%。從理論上分析，可以減少40%左右的計算量。

3.3 算法仿真分析

為了驗證本文2D-FFT快速算法的正確性，以及計算效率。在個人PC平臺上，使用Matlab軟件環(huán)境，完成了4組對比試驗，分別對64×64、128×128、256×256、512×512的圖像進行了測試，測試結(jié)果如表2所示。

表2 二維FFT處理時間對比

可以看出，本文方法相對于常規(guī)方法速度有顯著的提升，速度提升30%左右，接近于上面給出的理論提升量。實驗表明改進算法具有快速處理的優(yōu)點。圖4為實驗得到的頻域數(shù)據(jù)圖，(a)和(d)為常用的測試圖Cameraman和Plane，(b)和(e)為本文算法得到的中間結(jié)果，(c)和(f)為恢復(fù)得到的最后結(jié)果。

圖4 實驗結(jié)果圖

4 OpenCL實驗結(jié)果及其分析

4.1 OpenCL算法實現(xiàn)

傳統(tǒng)算法如圖5所示[4]。一般的FFT采用三層循環(huán)，通常對于N點變換，可以分為M級，N=2M，每一級有N/2個蝶形單元。當(dāng)前組完成后，完成下一組；當(dāng)前級完成后，完成下一級；傳統(tǒng)的方法都是串行的。如圖5所示的8點的FFT變換，共分為3級，每一級間有4個蝶形運算，對于蝶形間的計算和級間的運算都是串行依次進行的。分析可知：在同一級中具有相同旋轉(zhuǎn)因子的蝶形單元，可以通過備份多個旋轉(zhuǎn)因子以達到并行執(zhí)行；在同一級中具有不同旋轉(zhuǎn)因子的蝶形單元，因為數(shù)據(jù)不存在相關(guān)依賴性，本身就是可以并行執(zhí)行的；而對于相鄰級使用的旋轉(zhuǎn)因子也并非完全不同，根據(jù)這個規(guī)律可以將旋轉(zhuǎn)因子存儲為二維查找表來實現(xiàn)，減少計算量。

圖5 N=8時，蝶形運算圖

本文旨在將傳統(tǒng)的算法通過OpenCl在GPU平臺上進行并行加速，因此需要找出算法中可以并行的部分，對于算法并行粒度的劃分如圖6所示。對于二維FFT，局部工作組通常是一行或者幾行的大小，局部工作組之間是并行的；由若干個線程組成一個線程塊，然后由若干個線程塊組成一個局部工作組，本文中最小的處理單元是一個線程，而一個線程對應(yīng)著一個蝶形計算單元，即一個線程處理兩個點的數(shù)據(jù)。線程塊之間是并行的，算法中每個蝶形單元也是并行的，對于每個數(shù)據(jù)點都是并行計算的，F(xiàn)FT算法實現(xiàn)了真正的數(shù)據(jù)并行。試驗中將旋轉(zhuǎn)因子放到全局內(nèi)存表中，將待處理的數(shù)據(jù)通過全局內(nèi)存讀取到局部內(nèi)存中。

圖6 算法并行粒度和存儲共享

在OpenCL模型中，要取得足夠高的加速性能，需要非常多的線程才能夠提高計算效能，隱藏計算的訪存延時，本文根據(jù)OpenCL多線程執(zhí)行模式的特點，對FFT算法的分析如圖7所示。將算法按級劃分成M級，N=2M，每一級有N/2個蝶形運算單元，因為每一級的蝶形單元數(shù)據(jù)沒有相關(guān)性，可以完全并行執(zhí)行；在每一級間是數(shù)據(jù)通信，實驗中采取原址計算，即將前級蝶形運算的結(jié)果繼續(xù)存放在原來位置，而次級運算的時候，按照索引的地址從新位置讀取數(shù)據(jù)，將上一級產(chǎn)生的結(jié)果輸出作為后一級的數(shù)據(jù)輸入，實現(xiàn)數(shù)據(jù)在各級操作級中流動起來。實驗中分為log2N級，每一級中有N/2個線程塊，每一個線程塊完成一個蝶形運算，當(dāng)每一級的蝶形運算完成后通過柵攔函數(shù)同步，數(shù)據(jù)流入下一級。然后下一級進行同類運算，直至完成M級計算。另外，每一個蝶形單元的兩個輸入操作數(shù)的地址根據(jù)線程號和當(dāng)前級數(shù)共同確定。同步柵攔barrier是為了確保線程間的同步，保證數(shù)據(jù)的準(zhǔn)確性。采用OpenCl提供的庫函barrier(CLK_LOCAL_MEM_FENCE)實現(xiàn)線程塊內(nèi)的線程同步，它保證線程塊中的線程執(zhí)行到barrier處才繼續(xù)執(zhí)行后面的語句。通過同步可以避免當(dāng)一個線程塊中的一些線程訪問共享存儲器的同一地址時可能發(fā)生的讀寫錯誤問題。

圖7 算法多線程并行示意圖

4.2 存儲層次的優(yōu)化

① 局部存儲器是GPU片上的高速存儲器，是能夠被相同線程塊中的線程同時訪問的可讀寫存儲器。在實驗中，按照圖7，將算法流程分為多級，上級的數(shù)據(jù)輸出作為下一級的輸入，數(shù)據(jù)完全在級間流動，數(shù)據(jù)操作采用同址方式。如果將輸入數(shù)據(jù)從全局內(nèi)存搬移到局部內(nèi)存中，就可以大大提高線程訪問數(shù)據(jù)的速度，通常來說全局內(nèi)存具有非常高的訪存延時。在內(nèi)核函數(shù)設(shè)計中，首先將顯存所需的數(shù)據(jù)搬遷到局部內(nèi)存中，最后處理完成后，再從局部內(nèi)存搬遷至顯存中，計算過程，只在局部內(nèi)存中完成，而訪問全局內(nèi)存也只需要兩次，大大提高了計算效率。

② 進程的充分利用。在實驗中針對N點FFT變換，實際上只用到了N/2個線程，將局部工作組劃分為一行的時候，就會造成N/2個空線程的資源浪費。對于一個線程完成一個蝶形運算單元，也就是完成兩點的數(shù)據(jù)變換，因此利用此點可以提高線程的利用率。實驗中利用N個線程完成2N點的變換，每N/2個線程分別完成一行的變換，充分利用了線程，提高了計算資源的利用率，實驗速度得到了顯著提升。

③ 內(nèi)存的讀寫優(yōu)化。針對2D-FFT算法，有大量的存儲數(shù)據(jù)讀取，因此提高速度，可以對數(shù)據(jù)的讀寫進行優(yōu)化。實驗中讀數(shù)據(jù)的讀寫速度進行測試，對于128×128的數(shù)據(jù)進行讀取，從全局存儲器讀到共享存儲器需要13μs，不論是按行讀取還是按列讀??；而從共享存儲器寫到全局存儲器則需要分為按行寫和按列寫，按行寫為32μs，按列寫則是85μs，可以發(fā)現(xiàn)寫的時候按行寫會大大節(jié)省時間。在2D-FFT變換中，行變換完成后需要進行一個轉(zhuǎn)置，接著進行列變換，對讀寫時間的分析得出最大耗時的地方是轉(zhuǎn)置，所以實驗中采取的做法是在行變換完成后仍然按行寫到全局內(nèi)存，在列變換的時候任然采取按列讀，這樣極大地減小了按列寫的時間和轉(zhuǎn)置的時間，實驗速度得到了近40%的提升。

④ 內(nèi)存的寫優(yōu)化。實驗中發(fā)現(xiàn)在寫的過程中，同樣是按行寫，但是按行順序?qū)懞桶葱衼y序?qū)?，仍然存在著速度的差異，因此在寫的過程中，預(yù)先對數(shù)據(jù)排序后，然后依次往全局內(nèi)存寫數(shù)據(jù)，實驗結(jié)果有一定的提升。

4.3 實驗結(jié)果

本文在基于Windows操作系統(tǒng)上和AMDGPU實現(xiàn)了該算法，驗證了OpenCL的FFT算法實際的性能。實驗中所用到的平臺和配置如表3所示，程序基于AMD的OpenCl開發(fā)包，并在OpenCl1.2的環(huán)境下編譯。為了驗證對比本文算法的性能，和文獻[4]基于CUDA版本的算法以及基于CPU版本的FFT算法做了對比，其中文獻[4]的處理器和顯卡配置如表4所示。

表3 本文的平臺配置

表4 文獻[4]的實驗硬件配置

現(xiàn)有的GPGPU通用計算的處理流程一般為：CPU端作為宿主機先處理少部分邏輯運算，然后通過調(diào)用API函數(shù)將大量的計算任務(wù)交給OpenCL設(shè)備(GPU)。但是內(nèi)存和顯存的實際位置不一樣，所以CPU端通過內(nèi)存對象將操作數(shù)據(jù)送到GPU顯存中，等到GPU處理完畢，然后再通過內(nèi)存對象送回到CPU內(nèi)存中。實驗中的計算時間為FFT算法的總運算時間，包括在CPU上執(zhí)行邏輯函數(shù)的時間和在GPU上執(zhí)行kernel函數(shù)的時間，以及從CPU拷貝數(shù)據(jù)到GPU顯存的時間。表5對比了本文和其他平臺的實驗結(jié)果。

表5 FFT算法的運行總時間 單位：us

實驗數(shù)據(jù)表明：當(dāng)圖像數(shù)據(jù)小的時候，GPU和CPU之間數(shù)據(jù)相互拷貝的時間會抵消很大一部分并行加速的時間，所以加速比很小，當(dāng)數(shù)據(jù)規(guī)模增大的時候，這種加速效果就逐步提升。本文中AMD Radeon HD 7450支持最大工作組大小為256，所以只針對最大圖像數(shù)據(jù)為512×512。從實驗數(shù)據(jù)可以看出，文獻[4]基于CUDA的方法和本文方法結(jié)果基本一致，相比較CPU版本有4～7倍的加速。但是考慮到NVIDIA GTX 465浮點運算峰值為1 277.7 GFLOPs，Inter i7 core 3770浮點運算峰值為282 GFLOPs，而AMD Radeon HD 7450的浮點運算峰值僅為300 GFLOPs，因此本文的算法相對于文獻[4]有4倍的加速，相對于CPU版本有真正意義上的4～7倍的加速。對于同樣浮點運算能力的處理器來說，這個加速效果還是很明顯的。

5 結(jié) 語

本文實現(xiàn)了基于OpenCL的高速FFT計算。通過分析2D-FFT算法的可并行性特點，將算法分級，多線程分塊，充分利用了數(shù)據(jù)間的非相關(guān)性，實現(xiàn)了數(shù)據(jù)的充分流動，極大地提高了并行度；針對OpenCL的存儲層次特點和算法層次的優(yōu)化，明顯地提高了運算的速度，和相同浮點運算能力的平臺相比，具有4倍左右的良好加速比。

[1] AMD上海研發(fā)中心.跨平臺的多核與眾核編程講義OpenCL的方式[M].2010.

[2] Aaftab Munshi,benedict R Gaster,Timothy G Mattson,et al.OpenCL Programming Guide[M].Pearson Education,2012.

[3] AMD Corporation.Accelerated Parallel Processing OpenCL TM[Z].Jan.2011.

[4] 趙麗麗,張盛兵,張萌,等.基于CUDA的高速FFT計算[J].計算機應(yīng)用研究,2011,28(4):1556-1559.

[5] 張櫻,張云泉,龍國平,等.基于OpenCL的圖像模糊化算法優(yōu)化研究[J].計算機科學(xué),2012,39(3):259-263.

[6] Khronos group.OpenCL-The open standard for parallel programming of heterogeneous systems[OL].http://www.khronos. org/opencl/.

[7] 溫博,張啟衡,張建林.高分辨圖像二維FFT正反變換實時處理方法及硬件實現(xiàn)[J].計算機應(yīng)用研究,2011,28(11):4376-4379.

[8] 左顥睿,張啟衡,徐勇,等.基于GPU的并行優(yōu)化技術(shù)[J].計算機應(yīng)用研究,2009,26(11):4115-4118.

[9] 范進,金生震,孫才紅.超高速FFT處理器的設(shè)計與實現(xiàn)[J].光學(xué)精密工程,2009,17(9):2241-2246.

RESEARCH ON FFT ALGORITHM BASED ON OPENCL

Jia Ge1,2Peng Xianrong1Zuo Haorui1

1(InstituteofOpticsandElectronics,ChineseAcademyofSciences,Chengdu610209,Sichuan,China)2(GraduateUniversityofChineseAcademyofSciences,Beijing100039,China)

Fast Fourier transform, as a commonly used computational tool in the field of image processing, especially in image restoration algorithm, transforms the time-domain computation into frequency-domain computation and has great significance for engineering applications. By adopting thread blocks and parallel mapping method, we can make FFT algorithm reach the maximum degree of parallelism. In view of the storage features of OpenCL and the optimisation of the algorithm, the AMD GPU platform has been significantly accelerated. Compared with CPU platform and CUDA with the same processing capability, the performance of the optimised algorithm has increased by 7 times and 4 times respectively.

Fast Fourier transform OpenCL GPU Parallel speedup

2016-01-20。賈格，博士生，主研領(lǐng)域：并行加速，圖像盲復(fù)原。彭先蓉，研究員。左顥睿，副研究員。

TP302

A

10.3969/j.issn.1000-386x.2017.03.042