使用融合乘加加速快速傅里葉變換計算的向量化方法*

劉 仲，陳海燕，向宏衛(wèi)

(國防科技大學(xué)計算機(jī)學(xué)院，湖南長沙 410073)

現(xiàn)有的許多處理器體系結(jié)構(gòu)提供融合乘加(Fused Multiply-Add，F(xiàn)MA)指令，如 Intel的Itanium，IBM的Power處理器等。給定3個輸入操作數(shù)a，b，c，使用一條FMA指令即可實現(xiàn)y=±a±(b×c)中的任何一個操作。FMA指令對于數(shù)值計算亦具有重要的意義。因為，在軟件上，一條FMA指令在執(zhí)行時間上和一條乘法或加法指令幾乎一樣快;在硬件上，F(xiàn)MA單元通常比分開的乘法器和加法器耗費少;在計算上，F(xiàn)MA指令通過減少舍入操作可以提高數(shù)值計算精度。在提供FMA指令的處理器平臺上，對于有些數(shù)值計算應(yīng)用，轉(zhuǎn)換到FMA指令是比較直接的，如普通的矩陣和向量乘法、矩陣和矩陣乘法以及成對出現(xiàn)的乘法和加法操作計算。但是，對更多的其他應(yīng)用就沒有那么簡單了，需要平衡加法和乘法操作，改造原有的計算方法或流程以適合FMA指令，才能充分發(fā)揮FMA指令的作用，加速計算的性能。

向量處理器在單個芯片上集成多個向量處理單元(Vector Processing Element，VPE)，每個 VPE包含相同的MAC，ALU，BP等多個功能部件，能夠在降低面積和功耗的情況下大幅提高整芯片的計算能力，適合處理高密集運算的任務(wù)，如矩陣計算、FFT運算、濾波運算等。然而，現(xiàn)有的大量程序和算法是基于單核處理器設(shè)計的，很多高密集運算的任務(wù)由于算法本身的特性，向量化處理困難，如何針對向量處理器的多運算部件的向量計算的體系結(jié)構(gòu)特點，充分開發(fā)各個層次的并行性，高效地向量化這些算法是當(dāng)前面臨的主要困難［1-2］。

快速傅里葉變換(Fast Fourier Transform，F(xiàn)FT)作為時域和頻域的基本變換工具，在現(xiàn)代信號處理系統(tǒng)領(lǐng)域應(yīng)用廣泛，如雷達(dá)、聲吶、地震信號分析、頻譜分析、語音識別、3G通信和圖像處理等。FFT能顯著減少離散傅里葉變換(Discrete Fourier Transform，DFT)計算的復(fù)雜度，對于實時性要求很高的信號處理應(yīng)用來說，DFT計算效率越高，信號處理的實時性就越好，因此FFT算法的優(yōu)化方法一直是國內(nèi)外的研究熱點。Pease［3］采用數(shù)學(xué)工具Kronecker積描述FFT算法，方便將FFT算法高效地映射到并行計算機(jī)上;Linzer［4］，Goedecker［5］和 Karner［6］等針對 FMA 結(jié)構(gòu)的處理器，分別研究了如何利用FMA指令優(yōu)化FFT計算;Voronenko［7］依據(jù)有向無環(huán)圖和具有一定結(jié)構(gòu)的矩陣因子化算法，提出一種將離散傅里葉變換、離散余弦變換等線性變換轉(zhuǎn)化為FMA算法的一般方法;FFTW(Fastest Fourier Transform in theWest)［8］是通用CPU平臺上廣泛使用的FFT數(shù)學(xué)庫，移植性好，可以為任意長度和維數(shù)的FFT自動生成一種最有效的實現(xiàn)方式;Loberiras［9］等針對GPU提出小點數(shù)FFT的優(yōu)化技術(shù)，受限于紋理存儲器的大小，大點數(shù) FFT的性能不理想;Li［10］等針對 GPU提出一種 FFT自適應(yīng)優(yōu)化框架，當(dāng)點數(shù)大于1K時，獲得接近CUDA快速傅里葉變換(CUDA Fast Fourier Transform，CUFFT)的90%性能;文獻(xiàn)［11］針對GPU的大點數(shù)FFT優(yōu)化方法達(dá)到CUFFT性能的2．1倍。

在傳統(tǒng)的FFT計算中，蝶形單元的最終計算往往轉(zhuǎn)化為實數(shù)的乘法和加法操作，如時域抽取基2 FFT算法的一個蝶形單元計算需要4次實數(shù)乘法和6次實數(shù)加法，即需要10次實數(shù)乘(加)操作;時域抽取基4 FFT算法的一個蝶形單元計算需要12次實數(shù)乘法和22次實數(shù)加法，即需要34次實數(shù)乘(加)操作。顯然，F(xiàn)FT計算中的乘、加操作是不平衡的，不能直接利用FMA指令來實現(xiàn)蝶形單元的計算，所以傳統(tǒng)的FFT計算方法不能有效發(fā)揮具有FMA指令的處理器的計算效率。劉仲等針對FMA結(jié)構(gòu)的向量處理器，提出一種基于FMA加速FFT計算的向量化方法，通過變換FFT的蝶形單元運算流程，將傳統(tǒng)計算方式中的獨立的乘法和加法操作組合成次數(shù)更少的FMA操作，使得計算一個DIT基2 FFT算法的蝶形單元所需要的浮點操作減少到6次FMA操作，計算一個DIT基4 FFT算法的蝶形單元所需要的浮點操作減少到24次FMA操作。實驗結(jié)果表明，提出的方法能夠顯著加速FFT的計算性能。

1 向量處理器Matrix的體系結(jié)構(gòu)

如圖1所示，Matrix是一款面向高密度計算應(yīng)用的高性能浮點向量處理器，內(nèi)核是超長指令字(Very Long Instruction Word，VLIW)體系結(jié)構(gòu)，包括標(biāo)量處理部件(Scalar Processing Unit，SPU)和向量處理部件(Vector Processing Unit，VPU)。SPU負(fù)責(zé)標(biāo)量任務(wù)計算和流控，SPU和VPU可通過共享寄存器交換數(shù)據(jù)。Matrix每時鐘周期發(fā)射11條指令，包括5條標(biāo)量指令和6條向量指令。指令派發(fā)單元對執(zhí)行包進(jìn)行識別，并將其中的指令派發(fā)到相應(yīng)的功能單元中執(zhí)行。VPU負(fù)責(zé)向量計算，包括16個向量處理單元(Vector Processing Element，VPE)，每個VPE含一個局部寄存器文件，以及3個浮點乘加單元(Float Multiply and Accumulate，F(xiàn)MAC)、1個BP和2個L/S共6個并行功能部件，3個FMAC均支持FMA指令。局部寄存器文件包含64個64位寄存器，所有VPE的同一編號的局部寄存器在邏輯上又組成一個1024位的向量寄存器。功能部件支持定點和浮點操作，向量指令在各個VPE上同時獨立運行。向量數(shù)據(jù)訪問單元支持向量數(shù)據(jù)的Load/Store，提供大容量陣列向量存儲器(Array Memory，AM)，每周期同時支持2個Load/Store指令。

圖1 Matrix的體系結(jié)構(gòu)Fig．1 Architecture of Matrix

2 基于FMA的FFT向量化優(yōu)化方法

序列x(n)(n=0，…，N-1)的離散傅里葉變換 X(k)(k=0，…，N-1)定義為:

依據(jù)離散傅里葉變換公式計算的DFT復(fù)數(shù)乘加運算的計算復(fù)雜度是N2。當(dāng)N比較大時，運算量增長得非?？?，使得直接采用DFT計算方法難以滿足很多實際應(yīng)用的需求。1965年Cooley和Turkey提出一種快速傅里葉變換算法［12］，計算復(fù)雜度由原來的O(N2)降到O(N log2N)，可顯著地減少運算量。FFT算法分為時間抽取法(Decimation In Time，DIT)和頻率抽取法(Decimation In Frequency，DIF)兩種，現(xiàn)以 DIT 方法為例闡述FFT的向量化方法。

2．1 傳統(tǒng)的FFT蝶形單元計算方法

2．1．1 DIT基2FFT的蝶形單元計算方法

當(dāng)N是2的整數(shù)次方時，DIT基2FFT將輸入數(shù)據(jù)序列x(n)進(jìn)行奇、偶分組:

由旋轉(zhuǎn)因子的周期性特性易知:

令 a(l)=x(2l)，b(l)=x(2l+1)，則序列X(k)劃分為2個長度為N/2的子序列:

圖2是DIT基2FFT的蝶形單元運算流程圖，DIT基2FFT的每次蝶形單元運算需要1次復(fù)數(shù)乘法，2次復(fù)數(shù)加法，轉(zhuǎn)變實數(shù)計算即為4次實數(shù)乘法和6次實數(shù)加法，即需要10次實數(shù)乘(加)操作。

圖2 DIT基2 FFT的蝶形單元運算流程圖Fig．2 Radix-2 DIT FFT butterfly diagram

2．1．2 DIT基4FFT的蝶形單元計算方法

當(dāng)N是4的整數(shù)次方時，DIT基4FFT將輸入數(shù)據(jù)序列x(n)按模4后的余數(shù)分組:

由旋轉(zhuǎn)因子的周期性特性易知:

令 a(l)=x(4l)，b(l)=x(4l+1)，c(l)=x(4l+2)，d(l)=x(4l+3)，則序列 X(k)劃分為4個長度為N/4的子序列:

圖3是DIT基4FFT的蝶形單元運算流程圖，DIT基4FFT的每次蝶形單元運算需要3次復(fù)數(shù)乘法，8次復(fù)數(shù)加法，轉(zhuǎn)變實數(shù)計算即為12次實數(shù)乘法和22次實數(shù)加法，即需要34次實數(shù)乘(加)操作。

圖3 DIT基4 FFT的蝶形單元運算流程圖Fig．3 Radix-4 DIT FFT butterfly diagram

2．2 FMA優(yōu)化的FFT蝶形單元運算

2．2．1 FMA優(yōu)化的DIT基2FFT蝶形單元計算

假定DIT基2FFT的一個蝶形單元的兩個輸入分別為A和B，輸出分別為Y1和Y2，蝶形因子為W，下標(biāo)r和i分別表示復(fù)數(shù)的實部和虛部，根據(jù)式(1)，則有:

根據(jù)式(3)，DIT基2 FFT的蝶形單元運算可以分解成兩步完成:

1)計算中間結(jié)果ˉA，ˉB:

2)計算輸出結(jié)果Y1和Y2:

其中，步驟1的計算由2條FMA指令完成，步驟2的計算由4條FMA指令完成。從而計算一個DIT基2 FFT蝶形單元僅需要6條FMA指令操作，相比傳統(tǒng)的計算方式減少了4次。在傳統(tǒng)的計算方法下，計算N點的DIT基2 FFT需要的浮點操作次數(shù)為5N log2N;FMA優(yōu)化后需要的浮點操作次數(shù)減少到3N log2N，減少了40%。

2．2．2 FMA優(yōu)化的DIT基4FFT蝶形單元計算

假定DIT基4 FFT的一個蝶形單元的4個輸入分別為 A，B，C，D，輸出分別為 Y1，Y2，Y3，Y4，3個蝶形因子分別為W1，W2，W3，下標(biāo)r和i分別表示復(fù)數(shù)的實部和虛部，根據(jù)式(2)，則有:

根據(jù)式(4)，DIT基4 FFT的蝶形單元運算可以分解成兩步完成:

1)計算中間結(jié)果ˉA，ˉB，ˉC，ˉD:

其中，步驟1的計算由12條FMA指令完成，步驟2的計算由12條FMA指令完成。從而計算一個DIT基4FFT蝶形單元僅需要24條FMA指令操作，相比傳統(tǒng)的計算方式減少了10次。在傳統(tǒng)的計算方法下，計算N點的DIT基4 FFT需要的浮點操作次數(shù)為8．5N log4N;FMA優(yōu)化后需要的浮點操作次數(shù)減少到6N log4N，減少了29．4%。

2．3 混合基4和基2的FFT計算方法

如圖4所示，對于任意的N=2n，可采用混合基4和基2的FFT方法加速FFT的計算。若n是偶數(shù)，令n=2m，則N=22m=4m;若n是奇數(shù)，令n=2m+1，則 N=22m+1=4m×2。因此，N 點的FFT計算可轉(zhuǎn)化為m級基4FFT或者m級基4FFT和最后一級的基2 FFT計算。

圖4 混合基4和基2的FFT計算流程Fig．4 FFT computation diagram ofmixed radix-2/4

2．4 旋轉(zhuǎn)因子向量訪問的優(yōu)化

傳統(tǒng)的基4 FFT計算方法中，每個蝶形單元運算需要乘以3個不同的旋轉(zhuǎn)因子:和FMA優(yōu)化后的蝶形單元運算只需計算2個不同旋轉(zhuǎn)因子:WiN和W2Ni，并且預(yù)先計算出旋轉(zhuǎn)因子虛部/實部(Wi/Wr)的值，然后把旋轉(zhuǎn)因子按實部(Wr)、虛部/實部(Wi/Wr)交叉存儲放置。這樣，能減少旋轉(zhuǎn)因子存儲的個數(shù)。由于向量處理器只能連續(xù)存取向量數(shù)據(jù)，所以每級FFT運算的不同旋轉(zhuǎn)因子都要單獨放置。

同樣，DIT基2 FFT算法需要預(yù)先計算旋轉(zhuǎn)因子WiN和虛部/實部(Wi/Wr)的值，并且按實部(Wr)、虛部/實部(Wi/Wr)交叉的方式存放。

以1024點的DIT基4 FFT為例，算法由5級蝶形運算實現(xiàn)，且每一級不同旋轉(zhuǎn)因子所需的數(shù)目分別為:1個、4個、16個、64個和256個。由此可知，第1級和第2級的不同旋轉(zhuǎn)因子的數(shù)目少于VPE的個數(shù)，并且第1級需要乘的旋轉(zhuǎn)因子的值為1，所以可以省略。為了提高運算的并行度，可以把第2級所需的旋轉(zhuǎn)因子進(jìn)行冗余存放4次。

表1 1024點DIT基4 FFT的旋轉(zhuǎn)因子Tab．1 Twiddle factors of 1024-point radix-4 DIT FFT

3 性能測試與分析

在向量處理器Matrix上對不同點數(shù)的FFT計算性能進(jìn)行了測試(稱 MatrixFFT)，并與 Intel，AMD，IBM 平臺上的 MKL，ACML，ESSL 和 FFTW 3．0算法庫的性能進(jìn)行比較。實驗平臺中，Matrix的主頻為1GHz(雙精度峰值性能為96GFLOPS)，MatrixFFT的測試結(jié)果為RTL級仿真環(huán)境下的測試性能數(shù)據(jù);比對實驗的 CPU中，Intel選擇主頻3．0GHz Intel Xeon Core Duo(Woodcrest)(峰值性能為12GFLOPS);AMD選擇主頻2．2 GHz Dual Core AMD Opteron(峰值性能為4．4GFLOPS);IBM選擇主頻2GHz PowerPC 970(峰值性能為10GFLOPS)，其性能數(shù)據(jù)選自FFTW網(wǎng)站報告的評測數(shù)據(jù)［8］。

圖5和圖6分別給出了在Matrix上測試的不同點數(shù)的基2和基4FFT的單精度和雙精度計算性能。從圖5、圖6中可以看出，在點數(shù)較小時，F(xiàn)FT的性能較低，隨著點數(shù)的增大，F(xiàn)FT的計算性能顯著提高，16K點的單精度基2和基4FFT性能分別達(dá)到74GFLOPS和98GFLOPS，計算效率分別為38．5%和51．04%。16K點的雙精度基2和基4FFT 性能分別達(dá)到48．83GFLOPS 和49．05GFLOPS，計算效率分別為50．08%和51．09%。

圖5 基2FFT的計算性能Fig．5 Performance of radix-2 FFT

圖6 基4FFT的計算性能Fig．6 Performance of radix-4 FFT

圖7 對比了不同點數(shù)的雙精度FFT在 FMA優(yōu)化前后的計算性能，從圖7中可以看出，無論是DIT基2FFT，還是 DIT基4FFT，經(jīng)過 FMA優(yōu)化后，F(xiàn)FT的計算性能均顯著提高，其中 DIT基2FFT的計算性能平均提高30．4%，DIT基4FFT的計算性能平均提高30．3%。

圖7 FMA優(yōu)化的FFT性能對比Fig．7 Performance comparison of FMA optimized FFT

表2 不同處理器平臺下FFT性能(MFLOPS)Tab．2 FFT performance of different processors(MFLOPS)

表3 不同處理器平臺下FFT計算效率Tab．3 FFT efficiency of different processors

表2和表3分別從計算性能和效率兩方面給出了不同點數(shù)的雙精度FFT分別在Matrix，Intel，AMD，IBM 平 臺 上 的 MKL，ACML，ESSL 和FFTW3．0算法庫測試結(jié)果。

從絕對計算性能上看，MatrixFFT的性能遠(yuǎn)超其他幾種算法庫，主要原因是Matrix的峰值性能遠(yuǎn)超對比CPU的峰值性能;另一方面也體現(xiàn)出提出的基于FMA優(yōu)化的FFT計算效率較高，能夠充分發(fā)揮Matrix的計算性能。這一點從表3可以看出，在點數(shù)較小時，MatrixFFT的計算效率較低，這是因為Matrix是向量處理器，點數(shù)較小時，軟件流水中的循環(huán)填充和排空的開銷在整個計算中的占比較高，影響計算效率。在點數(shù)超過1024點以后，MatrixFFT的計算效率顯著提高，在4K點以后，計算效率都是最高的，其中16K點效率達(dá)50．86%，而對比的其他算法庫效率分別是47．88%，36．69%，26．83%，29．27%，12．28%。

4 結(jié)論

本文針對 FMA結(jié)構(gòu)的向量處理器，提出FMA加速FFT計算的向量化方法。通過優(yōu)化和重組FFT算法的蝶形單元運算流程，將原本不平衡的乘法和加法操作組合成融合乘加操作，利用FMA指令減少FFT計算的浮點操作指令次數(shù)，進(jìn)而提高FFT算法的計算性能和向量處理器的計算效率，提高硬件資源的利用率，加速FFT算法的計算性能。

References)

［1］ 劉仲，陳躍躍，陳海燕．支持任意系數(shù)長度和數(shù)據(jù)類型的FIR濾波器向量化方法［J］．電子學(xué)報，2013，41(2):346-351．LIU Zhong，CHEN Yueyue，CHEN Haiyan．A vectorization of FIR filter supporting any length and data types of coefficients［J］．Acta Electronica Sinica，2013，41(2):346-351．(in Chinese)

［2］ 劉仲，邢彬朝，陳躍躍．一種面向多核處理器的高效并行PCA-SIFT算法［J］．國防科技大學(xué)學(xué)報，2012，34(4):103-107．LIU Zhong，XING Binchao，CHEN Yueyue．An efficient parallel PCA-SIFT algorithm for multi-core processor［J］．Journal of National University of Defense Technology，2012，34(4):103-107．(in Chinese)

［3］ Pease M C．An adaptation of the fast Fourier transform for parallel processing［J］．Journal of the ACM，1968，15(2):252-264．

［4］ Linzer EN，F(xiàn)eig E．Implementation ofefficient FFT algorithms on fused multiply-add architectures［J］．IEEE Transactions on Signal Processing，1993，41(1):93-107．

［5］ Goedecker S．Fast radix 2，3，4，and 5 kernels for fast Fourier transformations on computers with overlapping multiply-add instructions［J］．SIAM Journal on Scientific Computing，1997，18(6):1605-1611．

［6］ Karner H，Auer M，Ueberhuber CW．Multiply-add optimized FFT kernels［J］．MathematicalModels and Methods in Applied Sciences，2001，11(1):105-117．

［7］ Voronenko Y，Puschel M．Mechanical derivation of fused multiply-add algorithms for linear transforms［J］． IEEE Transactions on Signal Processing，2007，55(9):4458-4473．

［8］ Frigo M，Johnson S G．BenchFFT［EB/OL］．［2014-03-15］．http://www．fftw．org/benchfft/．

［9］ Lobeiras J，Amor M，Doallo R．Influence of memory access patterns to small-scale FFT performance［J］． Journal of Supercomputing，2013，64(1):120-131．

［10］ Li Y，Zhang Y Q，Liu Y Q，et al．MPFFT:An auto-tuning FFT library for OpenCL GPUs［J］．Journal of Computer Science and Technology，2013，28(1):90-105．

［11］ 何濤，朱岱寅．大點數(shù)一維FFT的GPU設(shè)計實現(xiàn)［J］．計算機(jī)工程與科學(xué)，2013，35(11):34-41．HE Tao，ZHU Daiyin．Design and implementation of largepoint1D FFT on GPU［J］．Computer Engineering ＆ Science，2013，35(11):34-41．

［12］ Cooley JW，Turkey JW．An algorithm for the machine calculation of complex Fourier series［J］．Mathematics of Computation，1965，19:297-301．