魂芯DSP高效訪存并行FFT算法研究

楊振浩，鄭啟龍，鄧文齊，王向前

1(中國科學技術大學 計算機科學與技術學院，合肥 230027)2(中國電子科技集團公司第三十八研究所，合肥 230088)

1 概 述

FFT(Fast Fourier Transform)是信號處理、音頻分析等領域的基本研究工具，其性能的優劣影響著這些領域的發展.自1965年由Cooler-Tukey提出FFT算法以來[1]，FFT算法的研究取得了很大的進展，主要包括設計專門的處理器應對FFT計算[2-4]、改進FFT存儲結構[5]、針對具體處理器的硬件特性，改進FFT計算過程[6-8]等.其中，He S，Torkelson M等人[2，3]和Yang Z X等人[4]分別針對FFT和IFFT設計單獨的處理器；Ma Y等人[5]和馬瀟等人[9]對FFT的訪存過程進行優化；張杰等人[8]和李軍等人[10]在FFT計算過程中引入SIMD指令等方法都取得了很好的效果.

BWDSP100是一款自主設計的支持VLIW(Very Long Instruction Word，超長指令字)和SIMD技術的數字信號處理器，主要面向高性能計算、高端嵌入式等領域.BWDSP提供了向量化訪存、向量化計算、復數指令、位反序尋址指令等高效指令.由于通用CPU上的FFTW[11](Fast Fourier Transform in the West)等軟件包很難充分利用BWDSP處理器的硬件資源，通用向量化FFT算法未能在BWDSP處理器上取得令人滿意的性能表現.

本文結合BWDSP處理器的體系結構特征對FFT算法進行研究，利用BWDSP支持向量指令、位反序尋址指令、復數指令等高效指令的特性，提出了有效減少旋轉因子重復計算或訪存，并支持向量計算、向量訪存的高效FFT算法.實驗結果表明，優化后的FFT算法，平均性能達到7.61Gflops，分別是相同環境下串行FFT的16.54倍，普通向量化FFT的4.03倍.

2 背景介紹

2.1 魂芯DSP處理器

BWDSP系列處理器由CETC38研究所研制，可廣泛運用于各種高性能計算領域.作為本系列的第一款產品，BWDSP100是一款32位浮點DSP，采用VLIW、SIMD(Single Instruction Multiple Data，單指令多數據流)架構，具有強大的并行處理能力.BWDSP100處理器內核結構稱為eC104，核內部包含4個基本執行宏(Element operation Macro，簡稱宏)，每個宏由8個ALU、4個MUL、2個SHF、1個SPU和1個通用寄存器組組成.內部數據讀總線位寬為512bit，內部數據寫總線位寬為256bit，工作主頻500MHz，體系結構如圖1所示.可以把BWDSP看作具有256位向量寄存器及其計算能力的向量處理器.

圖1 BWDSP100 的結構框圖Fig.1 Architecture of BWDSP100

BWDSP系列處理器提供了高效的復數操作指令，大大減少了復數操作的指令周期.如：

浮點復數向量加：

{x，y，z，t}CFRs+1：s = CFRm+1：m + CFRn+1：n (減法操作類似)

浮點復數向量乘：

{x，y，z，t}QFRm+1：m_n+1：n = CFRm+1：m * CFRn+1：n

{x，y，z，t}Rs+1 = Rm+1 - Rn+1 || {x，y，z，t}Rs = Rm + Rn

位反序尋址讀：

{x，y，z，t}Rs+1：s br(C) = [Un+=Um，Uk]

位反序尋址寫：

br(C)[Un+=Um，Uk] = {x，y，z，t}Rs+1：s

其中位反序尋址是專門為FFT運算設定的一種特殊尋址方式，它將發出的存儲器地址的若干位前后顛倒，形成存儲器實際所需要的地址.與普通的位反序指令相比BWDSP處理器位反序尋址指令計算存儲器實際所需要的地址以完成向量訪存操作，而不是單純計算某數值的位反序值.另外需特別指出的是，反序操作僅僅是指產生地址的基地址進行反序操作，而參與地址修改的基地址并不做反序操作，這為向量化訪存提供了保證，也是本文高效訪存并行算法的立足點.

2.2 FFT算法

離散傅里葉變換(Discrete Fourier Transform，DFT)是數字信號領域的基石之一，有如下定義[12]：

(1)

1965年由Cooley和Tukey提出了DFT的一種快速算法，將傅里葉變換的時間復雜度從O(N2)降為O(Nlog2N)[1].目前FFT算法基本上可以分成兩大類，即按時間抽選(DIT)法和按頻率抽選(DIF)法.本文針對基-2 DIT算法展開.

圖2 倒位序輸入FFT算法Fig.2 Inverted sequence input FFT algorithm

3 FFT算法的優化

3.1 FFT算法向量化分析

對于任何按時間抽選的基-2 FFT流圖，如圖3所示，只要保持各節點所連的支路及其傳輸系數不變，則不論各蝶形運算單元計算順序如何，所得最后結果都是正確的[12]，即每一級中任意兩個蝶形單元運算是相互獨立的，可以并行執行.鑒于這種特性，我們可以充分利用魂芯DSP對SIMD指令的支持，實現蝶形單元的向量化運算.

圖3 8點倒位序輸入，自然序輸出FFT流程圖Fig.3 Inverted order input，natural sequence output FFT flow chart of eight points

魂芯DSP100處理器對單個block塊訪存帶寬為256位，采用將旋轉因子與輸入數據存放在不同block塊以充分利用其帶寬，4個簇可一次存取4個單精度復數.圖4所示的蝶形運算向量單元總共進行了4個蝶形運算.

圖4 向量運算單元Fig.4 Vector computing unit

算法2.首先向量雙字讀取4個蝶形單元第一個數據到4個簇相鄰寄存器(不妨設為r1：r0)，讀取蝶形單元第二個數據到寄存器(不妨設為r3：r2)，分4次雙字讀取旋轉因子到各個簇，然后運用3條復數向量指令(即圖4中的OP操作)即可完成4個蝶形單元運算，最后由兩條向量雙字寫指令將結果保存在原地址.向量化前，每個蝶形單元需要3條復數運算(加減乘)指令，5條雙字訪存指令，經VLIW(Very Long Instruction Word)合成優化后需要9個時鐘周期，4個蝶形單元至少需要36個時鐘周期.向量化后，每4個蝶形單元需要3條SIMD復數運算(加減乘)指令，4條SIMD訪存指令，4條雙字訪存指令，經VLIW合成優化后需要10個時鐘周期.運算指令條數下降了75%，訪存指令條數下降了60%，指令周期數下降了72%.其中訪存指令條數下降不明顯是由于不能向量化讀取4個不同的旋轉因子，但通過將旋轉因子放在不同block塊，通過VLIW指令合成，其中兩個旋轉因子讀能與兩次輸入數據讀并行執行，故指令周期數仍有較大幅度下降.

3.2 SIMD位反序尋址多階合并FFT算法

上述算法受限于訪存帶寬的限制，未能充分發揮BWDSP計算性能.BWDSP單簇單周期內能夠完成8次浮點加減法及4次浮點乘法運算.單個形如(a+bj)±(c+dj)*(Wr+Wj)的蝶形運算只需要4次乘法和6次加減法即可完成，可見BWDSP處理器能在單周期內實現4次浮點蝶形運算.而上述算法只能在4個周期內實現4次蝶形運算，即單周期1次蝶形運算.主要原因如下：首先，蝶形運算向量單元的旋轉因子為4個不同復數，只能分4次分別讀取4個旋轉因子，此處只用到正是帶寬的1/4，故需要新的數據排列方式，同化旋轉因子；其次，在每階計算完成后，需要將結果保存，并在下階迭代中取出，此操作將產生大量不必需的訪存，因此有必要實現多階合并的FFT算法[9]，減少數據IO操作；最后，三層循環FFT算法中，最內層循環迭代次數是動態變化的，編譯器不容易進行軟流水優化，故將內兩層循環轉換為一層固定迭代次數的循環，便于軟流水優化.

void_br_write8(int*，int*，int*，int*，Complex_Float4);

void_br_read8(Complex_Float4，int*，int*，int*，int*);

圖5 8點自然序輸入，倒位序輸出FFT流程圖Fig.5 Natural sequence input，Inverted order output FFT flow chart of eight points

采用m階合并時，訪存操作由原來的2Nlog2N減少為2Nlog2N/m，但受限于寄存器壓力的限制，需保證最內層循環每次迭代保存輸入數據的寄存器個數Ninput、保存旋轉因子寄存器個數Nw、保存臨時數據寄存器個數Ntmp之和小于等于寄存器總個數Nregisters，即滿足下式：

Ninput+Nw+Ntmp≤Nregisters

(2)

m與寄存器使用數量關系滿足下式：

Ninput=2m+1，Nw=2m，Ntmp=2m+1

(3)

代入(2)式可得：

2m+2+2m≤Nregisters

(4)

已知魂芯DSP單簇寄存器個數64個，當m≥4時，最內層循環單次迭代單簇使用的寄存器數量超過了其寄存器總量，故m最大等于3.當m=3時，使用的寄存器數量為40個，超過了DSP寄存器總數的一半，阻礙了軟流水優化的進行，因此我們最終采用了2階合并的FFT算法，此時寄存器使用量為20個，可以較好的進行軟流水優化，且計算資源使用足夠充分.Input:A,N

圖6 高效訪存并行FFT算法Fig.6 Efficient access to parallel FFT algorithm

圖6算法3.詳細描述了高效訪存并行FFT算法流程.

3.3 邊界處理

上述算法中我們每次選取的蝶形單元中，必須保證至少4個旋轉因子是相同的.而對于原序輸入FFT流圖的最后兩階分別有2個、1個相同的旋轉因子，對此我們可以在此處進行位反序，利用逆序輸入的向量化FFT算法進行計算.由于用到了兩階合并，必須保證在原序輸入的位反序尋址FFT算法中階數為偶數階.DIT算法中，前兩級為特殊旋轉因子，第一級全部為1，第二級為1和-j交替.對于這兩種特殊旋轉因子，都不需要進行乘法運算，單個蝶形運算的運算量僅為4次加減法.為此我們可通過判斷FFT輸入序列的奇偶性，若是奇數階，則進行前3階合并；若是偶數階，則進行前4階合并.

4 實驗及分析

本文將原序輸入的位反序尋址向量化FFT算法在BWDSP100處理器上進行實現和測試，并與相同實驗環境下的逆序輸入向量化FFT算法和標量FFT算法進行了對比.實驗結果表明，高效訪存并行FFT算法較標量FFT算法和逆序輸入向量化FFT算法都有很大提升.

測試所用的硬件平臺為魂芯DSP100處理器，工作頻率為500Mhz，編譯器為OCC—為充分支持BWDSP體系結構而開發的一款高效編譯器.算法1表示標量FFT算法，算法2表示普通向量化FFT，算法3表示高效訪存并行FFT算法.

圖7 算法1、算法 2 和算法 3 性能對比Fig.7 Algorithm 1，algorithm 2 and algorithm 3 performance comparison

圖7描述了BWDSP100處理器下，標量版FFT算法、普通向量化FFT算法和高效訪存并行FFT算法的運行結果.從圖7中可以看出，高效訪存并行FFT算法整體性能比標量FFT算法和普通向量FFT算法都有明顯提升.其中標量FFT算法完全沒有利用BWDSP的分簇架構性能最差，平均性能僅為460Mflops.普通向量化FFT算法作為一種通用的優化算法，在BWDSP上的平均性能為1.89Gflops.高效訪存并行FFT算法性能表現最好，平均性能為7.61Gflops，是標量版FFT算法的16.54倍，是普通向量化FFT算法的4.03倍.高效訪存并行FFT算法通過原序輸入的方式，減少了3/4旋轉因子的訪存，通過兩階合并的計算方式，減少了不必要的中間結果的存取，并在三層循環合并為兩層之后，充分發揮了BWDSP的VLIW特性，使軟流水優化成為可能.

5 總 結

FFT算法在通用CPU上的優化工作已經取得了很大的進展，但是在BWDSP上的相關研究工作較少.通用CPU上的優化算法往往不能發揮BWDSP的VLIW和分簇等體系結構特性.隨著BWDSP處理器的迅速發展，這一問題變得尤為突出.本文結合BWDSP體系結構和FFT算法特性，提出了基于位反序的并行FFT算法.實驗表明，該算法在BWDSP上取得了很大的性能提升，平均性能達到7.61Gflops.

接下來的工作將在本文的理論基礎上，針對多核魂芯DSP，構建高效的并行FFT算法的多核求解方式.由于FFT算法對數據訪問過于頻繁和變步長訪問，未能發揮BWDSP多塊同時訪存特性，故本文所提算法仍有提升空間，未來還需尋求一種更高效的數據布局方式.