GPU平臺二維快速傅里葉變換算法實現(xiàn)及應(yīng)用

張 全，鮑 華，饒長輝，彭真明( 1. 中國科學院自適應(yīng)光學重點實驗室，成都 61009；. 電子科技大學 光電信息學院，成都 610054; 3. 中國科學院光電技術(shù)研究所，成都 61009；4. 中國科學院大學，北京 100049 )

?

GPU平臺二維快速傅里葉變換算法實現(xiàn)及應(yīng)用

張 全1,2,3,4，鮑 華1,3，饒長輝1,3，彭真明2
( 1. 中國科學院自適應(yīng)光學重點實驗室，成都 610209；2. 電子科技大學 光電信息學院，成都 610054; 3. 中國科學院光電技術(shù)研究所，成都 610209；4. 中國科學院大學，北京 100049 )

摘要：NVIDIA在其GPU平臺上開發(fā)的FFT庫CUFFT經(jīng)過幾次升級，但在二維FFT實現(xiàn)上效率還有提升空間，而且對于特定不能與上下文的計算融合，導致多次對Global memory的訪問。本文分析合并內(nèi)存訪問事務(wù)大小與占用率之間的關(guān)系，優(yōu)化使用GPU存儲器資源，對小數(shù)據(jù)量2次冪二維復數(shù)FFT在GPU上的實現(xiàn)進行改進，加速比最高達到CUFFT 6.5的1.27倍。利用實數(shù)FFT結(jié)果的共軛對稱性，算法的效率比復數(shù)FFT算法運算量降低了40%。最后將FFT的改進應(yīng)用到光學傳遞函數(shù)(OTF)的計算中，采用Kernel 融合的方法，使得OTF的計算效率比CUFFT計算方法提高了1.5倍。

關(guān)鍵詞：快速傅里葉變換；CUDA；光學傳遞函數(shù)；圖形處理器

0 引 言

快速傅里葉變換(Fast Fourier Transform, FFT)是離散傅里葉變換的快速算法。1963年，J.W.Cooley和J.W.Tukey提出的Cooley-Tukey算法是FFT算法的早期版本，該算法以分治法為策略使DFT的運算量從減小到了，計算效率提高1～2個數(shù)量級。目前除了Cooley-Tukey算法以外，還涌現(xiàn)出許多高效的算法。其中包括素因子算法(Prime Factor)、分裂基算法，混合基算法等[1]。

在FFT實現(xiàn)方面，MIT開發(fā)了基于CPU的FFT算法庫FFTW(Fast Fourier Transform in the West，F(xiàn)FTW)，該庫函數(shù)具有很強的移植性和自適應(yīng)性，能夠自動配合所運行的硬件平臺，通過最優(yōu)化搜索，找到最優(yōu)的組合使運行效率達到最佳。2003年Kenneth Moreland與Edward Angel利用GPU的著色器編譯程序把FFT算法移植到了GPU平臺上[2]。2007年NVIDIA公司推出CUDA并行開發(fā)環(huán)境，同時也發(fā)布了基于CUDA 的FFT庫函數(shù)CUFFT，該庫進行FFT處理的速度大約是同期CPU進行FFT運行速度的20倍[3]。2008年，Vasily Volkov在GPU上對FFT算法進行改進使一維FFT算法效率比CUFFT 1.1提高3倍[4]。N.K. Govindaraju，B. Lloyd，Y. Dotsenko等人在Volkov的基礎(chǔ)上使得一維FFT效率比CUFFT 1.1提高4倍，并實現(xiàn)了二維FFT[5-6]。隨后，NVIDIA把Volkov等人改進的算法加入到下一版本的CUFFT中。

本文針對GPU多核處理器架構(gòu)，通過分析FFT的Cooley-Tukey算法框架，針對小數(shù)據(jù)量2的冪次方二維FFT在GPU上的實現(xiàn)進行分析，對列做一維FFT變換時，存在非合并內(nèi)存訪問的問題，使得訪存效率下降。通過分析調(diào)整Block的大小，使得同時讀入多列數(shù)據(jù)進行計算，盡量滿足128 Byte的Cache line，將非合并內(nèi)存訪問的影響降到最低。最后將該算法應(yīng)用到圖像處理中常用的光學傳遞函數(shù)(OTF)的計算中，通過Kernel融合的方法提高了OTF的計算效率。

1 算法理論

1.1 Cooley-Tukey DFT算法框架[7]

離散傅里葉變換DFT(Discrete Fourier Transform)將時域信號轉(zhuǎn)換成頻域信號。一維DFT計算式：

Cooley-Tukey算法主要包括兩次索引變換和兩次一維小點數(shù)DFT，其計算分為如下5步。

Step1: 輸入序列根據(jù)式(2)做輸入索引變換；

Step2: 計算N2個長度為N1點的一維DFT運算(第一級變換)；

Step4：計算N1個長度為N2點的一維DFT運算(第二級變換)； Step5: 對第二級變換結(jié)果按照式(3)做輸出索引變換。

1.2 實數(shù)FFT算法

設(shè)x(n)是2N點的實序列，現(xiàn)將x(n)分為偶數(shù)組x1(n)和奇數(shù)組x2(n)，即：

然后將x1(n)及x2(n)組成一個復序列y(n)=x1(n)+jx2(n)，通過N點FFT運算可得到Y(jié)(k)=X1(k)+jX2(k)，根據(jù)前面的討論，得：

其中：k=0,1,2,..., N -1。為求2N點x(n)所對應(yīng)的X(k)，需求出X(k)與X1(k)，X2(k)的關(guān)系，即：

而：

1.3 二維傅里葉變換

二維離散傅里葉變換的定義為

由一維DFT性質(zhì)知，當x(n)為實序列時，滿足復共軛對稱性，所以行向量做一維FFT結(jié)果只保存M/2+1列數(shù)據(jù)，如圖2(b)所示。然后對M/2+1列作N點DFT得到圖2(c)結(jié)果。

圖1 二維復數(shù)傅里葉變換計算過程Fig.1 Calculating procedure of 2-D FFT on complex number

圖2 二維實數(shù)傅里葉變換計算過程Fig.2 Calculating procedure of 2-D FFT on real number

1.4 GPU編程環(huán)境

在CUDA編程模型中，必須要有足夠多的活動線程來隱藏訪存延時，提高計算性能?？梢岳谜加寐屎饬炕顒泳€程數(shù)。但是程序優(yōu)化的過程中，高占用率不代表高效率[8]。高占用率雖然提高了并行度，需要在占用率與資源使用之間做出一種折中。

此外，在GPU上進行程序優(yōu)化主要遵循三個準則，其優(yōu)先級由高到底。首先是Kernel函數(shù)內(nèi)部應(yīng)盡量減少全局存儲器的訪問操作，同時全局存儲器訪問一定要滿足合并訪問；其次是設(shè)計應(yīng)盡量使用共享內(nèi)存(Share Memory, SM)進行Block內(nèi)部通信，同時對于共享存儲器的訪問讀寫也要防止共享存儲器的Bank沖突；最后算法設(shè)計Block內(nèi)部應(yīng)盡可能保持高密集度的寄存器計算，但同時也要協(xié)調(diào)好如何在片內(nèi)寄存器資源有限的情況下選擇最佳的Thread數(shù)量[9]。

2 基于GPU平臺的FFT工程實現(xiàn)

根據(jù)第一節(jié)中Cooley-Tukey算法思想，對于2m點數(shù)的FFT可以由基-2、基-4、基-8、基-16等多種組合實現(xiàn)，基數(shù)越高，數(shù)據(jù)重排的次數(shù)就越少，同時寄存器需求量也越大，在寄存器資源有限的情況下，高基數(shù)導致并行thread數(shù)受影響，降低了占用率。綜合上述因素，對于2m可以按順序優(yōu)先分解為基-16、基-8、基-4、基-2的組合。

對于一維復數(shù)FFT，根據(jù)向量的大小，一維復數(shù)FFT的實現(xiàn)分為三種情況：當向量長度為8或者16時，不需要進行線程之間的數(shù)據(jù)交換，其計算可以在寄存器中完成；當向量的長度為64、256、512、1 024、2 048和4 096時，線程之間數(shù)據(jù)交換需要在Share Memoy中完成；當向量的長度超過8 192時，必須通過兩次Kernel調(diào)用, 數(shù)據(jù)交換在Global Memory中完成。本文主要針對二維FFT計算進行分析，整個算法執(zhí)行流程如圖3所示，算法分為行方向一維傅里葉變換和列方向一維傅里葉變換。一維傅里葉變換又分為多級，相鄰的級間需要同步并且通過Share Memory交換數(shù)據(jù)。

圖3 二維傅里葉變換GPU實現(xiàn)過程Fig.3 Implementation of 2-D FFT on GPU

2.1 行方向一維FFT變換

對于小數(shù)據(jù)量的二維FFT算法，二維矩陣每一行數(shù)據(jù)可以放入Share Memory，不需要通過Global Memory進行數(shù)據(jù)交換。假設(shè)數(shù)據(jù)長度為N的一維向量，將其排列成N1×N2的二維矩陣。計算過程中首先選擇合適的N2，使其能夠在Share Memory中完成FFT計算。然后確定N1=N/N2，如果N1不能夠在Share Memory中完成計算，則對N1進一步分解為N1=N11×N12，以此類推。另外一種解釋就是相當于把長度為N的向量，重新排列為多維矩陣，每一維的長度為FFT的一個固定基，其FFT計算稱作整個N點FFT的一級。對N的每一次分解都要做一次矩陣轉(zhuǎn)置。以N=512為例，首先，512可以分解為64×8，先做64行8點一維FFT；其次，在Share Memory中進行數(shù)據(jù)轉(zhuǎn)置，然后，再做8行64點一維FFT，考慮到64點做FFT比較大，將其進一步分解為8×8。通常N存在多種分解例如512還可以分解為2×16×16，或者4×8×16，不同的分解效率也會有差別，基越小運算量越大，盡量避免基-2或者基-4的分解。分解的級越多同步操作越多，所以盡量在滿足共享存儲要求的前提下采用大基。依據(jù)以上原則，N=128的最優(yōu)分解為16×8而非4×4×8；N=256的最優(yōu)分解為16×16，而非8×8×4。

512點FFT存儲器訪問如圖4所示，完成一行512點一維FFT需要的線程數(shù)為64。為了滿足合并內(nèi)存訪問條件，512個數(shù)據(jù)被分為8組，每組64個數(shù)據(jù)。每個線程從8組數(shù)中各取1個數(shù)據(jù)存入8個寄存器中。為提高GPU的占用率，Block的大小取64的倍數(shù)，使一個Block同時完成多行512點數(shù)據(jù)的一維FFT。

2.2 列方向一維FFT變換

從GPU硬件的角度來看，GPU上的內(nèi)存延遲被不同線程束間的切換所隱藏。如果同一線程束的線程訪問相鄰內(nèi)存位置，并且內(nèi)存區(qū)域的開始位置是對齊的，這些訪問請求會自動組合或者合并。計算能量1.x設(shè)備上，合并內(nèi)存事務(wù)最小為128 Byte。如果被合并線程訪問的數(shù)據(jù)比較小，會導致內(nèi)存帶寬迅速下降。費米架構(gòu)的設(shè)備支持32 Byte和128 Byte的合并內(nèi)存事務(wù)。開普勒架構(gòu)的設(shè)備支持32 Byte、64 Byte和128 Byte的合并內(nèi)存事務(wù)。

列方向一維FFT變換與行方向一維FFT變換計算過程一樣。區(qū)別在于對全局內(nèi)存的訪問方式不一樣，在行方向一維FFT變換中線程訪問Global Memory是滿足合并內(nèi)存訪問的，而列方向一維FFT變換則不是。為了提高訪存效率，需要對存儲器訪問方式進行調(diào)整。以double數(shù)據(jù)類型512×512二維矩陣列FFT變換為例，矩陣中每個數(shù)據(jù)占16 Byte。本文實驗采用的設(shè)備為計算能力為3.5的開普勒計算卡，其支持32 Byte、64 Byte、128 Byte三種緩存行模式，為了滿足這三種內(nèi)存合并模式，每個Block分別至少處理2列、4列、8列一維FFT計算，經(jīng)過實驗測試每個Block處理4列一維FFT時效率最佳。512×512 double類型二維矩陣按列FFT變換存儲器訪問如圖5所示，線程數(shù)為64×4的二維矩陣，每一列線程訪問Global memory中對應(yīng)的一列數(shù)據(jù)，每個線程訪問列方向間隔為64的8個數(shù)據(jù)。

圖4 512×512行方向一維FFT變換全局存儲器訪問Fig.4 Global memory access of 1-D FFT on row for array size of 512×512

圖5 512×512列方向一維FFT變換全局存儲器訪問Fig.5 Global memory access of batched 1-D FFT on column for array size of 512×512

2.3 FFT在OTF計算中應(yīng)用[10]

光學成像系統(tǒng)在像面的復振幅可表示為

式中：F-1{}表示傅里葉逆變換，(u, v)和(x, y)分別為瞳面和像面坐標；P(u, v)為廣義光瞳函數(shù)，可以表示為

式中：j為虛數(shù)單位，φ(u, v)為波前相位，p(u, v)為孔徑函數(shù)。對應(yīng)成像系統(tǒng)的點擴散函數(shù)可表示為

而光學傳遞函數(shù)定義為點擴散函數(shù)的傅里葉變換，即

由上所述，OTF的計算分為6步：1) 廣義光瞳函數(shù)計算；2) 對光瞳函數(shù)列方向批量一維IFFT；3) 行方向批量一維IFFT；4) 對取模平方運行得到點擴散函數(shù)；5) 對點擴散函數(shù)列方向批量一維FFT；6) 行方向批量一維FFT。在GPU上實現(xiàn)時，由于1)和2)計算過程中每個數(shù)據(jù)相互獨立，可以合并，由一個Kernel完成計算；二維FFT計算中，交換行方向一維FFT和列方向一維FFT的計算次序不影響計算結(jié)果，交換5)和6)的次序，從而3) ～5)行與行之間數(shù)據(jù)相互獨立，三步可以合并由一個Kernel 完成。整個合并過程使得Kernel數(shù)由6個減小為3個。

3 實驗結(jié)果分析

3.1 實驗環(huán)境

本文實驗所采用的硬件配置：Intel(R) core(TM) i7-3930k主頻3.2 GHz CPU，內(nèi)存8 G，NVIDIA Telsa K20c GPU。軟件配置：Windows 7 64-bit操作系統(tǒng)，VS2010+CUDA6.5編程環(huán)境。

3.2 FFT并行設(shè)計性能比較

二維FFT的GPU實現(xiàn)分兩步執(zhí)行，首先進行行方向批量一維FFT，然后再列方向批量一維FFT，行方向計算時Global memory滿足合并內(nèi)存訪問，而列方向時不利于合并內(nèi)存訪問，根據(jù)不同的數(shù)據(jù)類型，盡量滿足相鄰線程訪問事務(wù)大小的數(shù)據(jù)。對列方向一維FFT時，一次讀入多列可以增加合并內(nèi)存數(shù)據(jù)塊，但是會消耗過多的共享內(nèi)存和寄存器資源從而降低了占用率，所以需要在內(nèi)存訪問效率與GPU占用率之間做出折中的選擇，針對不同的配置其占用率和內(nèi)存訪問效率都不一樣，表1分別針對Float2(Float型復數(shù)) 和Double2(Double型復數(shù))數(shù)據(jù)類型的列方向FFT在不同線程配置下做了測試。本實驗采用的Telsa K20c GPU，支持32 Byte、64 Byte、128 Byte的合并內(nèi)存事務(wù)。從實驗數(shù)據(jù)我們知道占用率不代表高性能，也不是內(nèi)存合并事務(wù)越大性能越高，需要在這兩者之間進行折中。對于256×256、512×512的二維FFT，64 Byte的合并內(nèi)存事務(wù)使得占用率不會最低，同時內(nèi)存訪問效率也不是最低，是一種較優(yōu)的配置。

表1 Float 2型和Double 2型數(shù)據(jù)Kernel 2不同線程配置性能比較Table 1 Performance comparison of Kernel2 for Float2 and Double2 under different thread configuration

對CPU平臺下FFTW庫函數(shù)，以及GPU平臺下CUFFT和本文方法，分別做1 000次二維FFT計算其平均值，三種方法都不考慮數(shù)據(jù)拷貝時間和內(nèi)存分配時間，另外FFTW和CUFFT庫函數(shù)時間只包括執(zhí)行時間，不包括配置時間。表2分別為數(shù)據(jù)大小為128×128、256×256、512×512的二維FFT執(zhí)行時間。從表中可以看出，數(shù)據(jù)較小時，CUFFT與FFTW的加速比較??；數(shù)據(jù)量變大時，加速比隨之增加，其范圍在2.5～53.7之間。對于128×128大小的二維FFT GPU資源并沒有充分利用，但是較之FFTW還是具有一定優(yōu)勢。本文的方法相對于CUFFT也有提高，雙精度512×512大小二維復數(shù)FFT計算效率是CUFFT的1.24倍。利用實數(shù)FFT計算結(jié)果的共軛對稱性，512×512的二維雙精度實數(shù)FFT計算效率是CUFFT的1.36倍。

表2 128×128、256×256、512×512二維FFT性能比較Table 2 Performance comparison of 2-D FFT for array of 128×128, 256×256, 512×512 μs

針對前面提到的OTF的計算，采用CUFFT庫函數(shù)方法和本文FFT分別對其實現(xiàn)，其性能如圖6所示，由于本文方法采用Kernel合并提高了GPU指令并行度，降低了Global Memory的訪問，對于512×512的雙精度二維OTF計算效率是采用CUFFT方法的1.51倍。

圖6 使用庫函數(shù)與本文方法計算OTF執(zhí)行時間比較Fig.6 Comparison of execution time between library method and this paper’ method

4 結(jié)論及展望

本文分析了Cooley-Tukey FFT計算框架，完成了基于CUDA編程模型的小數(shù)據(jù)量2的冪次方二維FFT 在GPU上的實現(xiàn)，在合并內(nèi)存訪問與GPU占用率之間進行折中處理，實驗結(jié)果表明，該方法效率是CUFFT 的1.24倍。OTF的GPU實現(xiàn)利用二維FFT計算特性，交換了一維行變換和列變換的次序，通過Kernel合并的方法，使OTF的計算效率相比CUFFT方法提高1.5倍。下一步工作主要針對大數(shù)據(jù)量的二維FFT并行化和非基-2 FFT的實現(xiàn)。

參考文獻：

[1] Rao K R，Kim D N，Hwang J J. 快速傅里葉變換：算法與應(yīng)用 [M]. 北京：機械工業(yè)出版社，2012：1-33.

[2] Moreland K，Angel E. The FFT on a GPU [C]// Proceedings of the ACM Siggraph/Eurographics Conference on Graphics Hardware，San Diego，California，July 26-27，2003：112-119.

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

ZHAO Lili，ZHANG Shengbing，ZHANG Meng，et al. High performance FFT computation based on CUDA [J]. Application Research of Computers，2011，28(4)：1556-1559.

[4] Volkov V，Kazian B. Fitting FFT onto the G80 Architecture [M]. University of California，2008，E63(40)：1-12.

[5] Naga K Govindaraju，Brandon Lloyd，Yuri Dotsenko，et al. High performance discrete Fourier transforms on graphics processors [C]// Proceedings of the 2008 ACM/IEEE conference on Supercomputing，Austin，Texas，November 15-21，2008：1-12.

[6] Brandon L D，Boyd C，Govindaraju N. Fast computation of general Fourier Transforms on GPUs [C]// IEEE International Conference on Multimedia and Expo，Hannover，Germany，June 23-26，2008：5-8.

[7] 楊麗娟，張白樺，葉旭楨. 快速傅里葉變換FFT及其應(yīng)用 [J]. 光電工程，2004，31(增)：1-3.

YANG Lijuan，ZHANG Baihua，YE Xuzhen. Fast Fourier transform and its applications [J]. Opto-Electronic Engineering，2004，31(Suppl)：1-3.

[8] Volkov V. Better performance at lower occupancy [C]// Proceedings of the GPU Technology Conference，San Jose，California，September 20-23，2010.

[9] Rob Farber. 高性能CUDA 應(yīng)用設(shè)計與開發(fā)：方法與最佳實踐 [M]. 北京：機械工業(yè)出版社，2013：1-27.

[10] Jim Schwiegerling. Relating wavefront error apodization，and the optical transfer function on-axis case [J]. Journal of the Optical Society of America A(S1520-8532)，2014，31(11)：2476-2483.

Realization and Application of Two-dimensional Fast Fourier Transform Algorithm Based on GPU

ZHANG Quan1, 2, 3,4，BAO Hua1,3，RAO Changhui1,3，PENG Zhenming2
( 1. Key Laboratory on Adaptive Optics, Chinese Academy of Sciences, Chengdu 610209, China; 2. School of Optoelectronic Information, University of Electronic Science and Technology of China, Chengdu 610054, China; 3. Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu 610209, China; 4. University of Chinese Academy of Sciences, Beijing 100049, China )

Abstract:NVIDIA as the inventor of the GPU provides a library function CUFFT for computing Fast Fourier Transform (FFT). After several generations update of CUFFT, there is still promotion space and it is not suit for kernel fusing on GPU to reduce the memory access and increase the Instruction Level Parallelism (ILP). We develop our own custom GPU FFT implementation based on the well-known Cooley-Tukey algorithm. We analyze the relationship of coalesce memory access and occupancy of GPU and get the optimal configuration of thread block. The results show that the proposed method improved the computational efficiency by 1.27 times than CUFFT 6.5 for double complex data 512×512. And then it is used to the computation of OTF with kernel fusing strategy, and it improved the efficiency of computation about 1.5 times than conventional method using CUFFT.

Key words:FFT; CUDA; optical transfer function (OTF); graphic processing unit (GPU)

作者簡介：張全(1985-)，男(漢族)，甘肅武威人。博士研究生，主要研究工作是GPU高性能計算，圖像處理。E-mail: quanzhang100@126.com。

基金項目：國家自然科學基金(11178004)；中國科學院實驗室創(chuàng)新基金(YJ14K018)

收稿日期：2015-02-04； 收到修改稿日期：2015-06-12

文章編號：1003-501X(2016)02-0069-07

中圖分類號：TP391

文獻標志碼：A

doi：10.3969/j.issn.1003-501X.2016.02.012