基于OpenCL的FFT算法研究

賈 格 彭先蓉 左顥睿

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

基于OpenCL的FFT算法研究

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

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

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

傅里葉變換 OpenCL GPU 并行加速

0 引 言

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

1 OpenCL簡介

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

圖1 OpenCl平臺模型

圖2 OpenCl應用執行模型

2 FFT簡介

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

3 FFT算法原理及其改進

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

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

(1)

(2)

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

(3)

(4)

(5)

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

(6)

(7)

同理可以得到：

(8)

(9)

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

3.2 算法復雜度分析

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

圖3 常規算法和本文算法處理流程

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

表1 算法復雜度對比

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

(10)

(11)

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

3.3 算法仿真分析

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

表2 二維FFT處理時間對比

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

圖4 實驗結果圖

4 OpenCL實驗結果及其分析

4.1 OpenCL算法實現

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

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

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

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

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

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

4.2 存儲層次的優化

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

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

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

④ 內存的寫優化。實驗中發現在寫的過程中，同樣是按行寫，但是按行順序寫和按行亂序寫，仍然存在著速度的差異，因此在寫的過程中，預先對數據排序后，然后依次往全局內存寫數據，實驗結果有一定的提升。

4.3 實驗結果

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

表3 本文的平臺配置

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

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

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

實驗數據表明：當圖像數據小的時候，GPU和CPU之間數據相互拷貝的時間會抵消很大一部分并行加速的時間，所以加速比很小，當數據規模增大的時候，這種加速效果就逐步提升。本文中AMD Radeon HD 7450支持最大工作組大小為256，所以只針對最大圖像數據為512×512。從實驗數據可以看出，文獻[4]基于CUDA的方法和本文方法結果基本一致，相比較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 結 語

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

[1] AMD上海研發中心.跨平臺的多核與眾核編程講義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].計算機應用研究,2011,28(4):1556-1559.

[5] 張櫻,張云泉,龍國平,等.基于OpenCL的圖像模糊化算法優化研究[J].計算機科學,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正反變換實時處理方法及硬件實現[J].計算機應用研究,2011,28(11):4376-4379.

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

[9] 范進,金生震,孫才紅.超高速FFT處理器的設計與實現[J].光學精密工程,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。賈格，博士生，主研領域：并行加速，圖像盲復原。彭先蓉，研究員。左顥睿，副研究員。

TP302

A

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