基于BWDSP100的高性能FFT實現

(中國電子科技集團公司第三十八研究所孔徑陣列與空間探測安徽省重點實驗室， 安徽合肥 230088)

0 引言

隨著數字技術的發展，基于時頻域二維處理的算法以其優越的性能逐漸成為信號處理算法的主流，時域-頻域間的相互轉換成為數字信號處理中重要的一環。自從20世紀60年代庫利圖基發展出快速傅里葉變換(FFT)以來，FFT一直是時頻域轉換的首要工具。正是基于對FFT運算在諸多領域中廣泛應用的重視， Freescale, TI, ADI等國際主流DSP芯片廠家紛紛采用FFT運行時間作為表征芯片性能的重要參數。對于一款高性能DSP而言，如何實現高效FFT運算是芯片走向實用的重要一環[1-2]。

1 芯片介紹

BWDSP100芯片(即“魂芯一號”)是中國電子科技集團公司第三十八研究所自主研發的一款高性能、超標量浮點運算DSP芯片。該處理器融合了多種現代芯片設計理念并加以擴展，充分發揮并行處理的優勢和特點，可以滿足大數據量、實時性要求高的數字信號處理、數據處理等工程應用要求。基于該芯片的雷達信號處理機已經完成設計，并在雷達數字陣列的大規模數字波束形成(DBF)中得到應用[3-5]。

BWDSP100采用宏作為基本處理模塊。其片內擁有4個處理宏(X/Y/Z/T)，每個宏擁有4個乘法累加器、8個ALU運算器、2個移位器、1個特定運算單元及64個32位寄存器，內嵌28 Mbit SRAM(24 Mbit數據空間，4 Mbit程序空間)，單宏單周期內能夠完成8次浮點加減法及4次浮點乘法運算。同時它也支持32位和16位定點運算。考慮到取數效率，BWDSP100內部擁有3個結構相同、相互獨立工作的地址發生器(U/V/W)，可以通過來自譯碼器的立即數或者來自通用寄存器組的數據，對地址發生器的內部寄存器進行初始化。地址發生器利用這些初始化信息作為地址計算所需的基地址和偏移量。整個芯片內部架構如圖1所示。

圖1 BWDSP100內部架構

BWDSP100內部24 Mbit數據空間采用SRAM搭建而成，不存在Cache命中問題。整個存儲部件分為3組(memory block)，分別稱為B0，B1，B2，每組又由8個獨立存儲器組成，每個存儲器的存儲容量為32 K×32 bit。因而器件內部數據存儲器總數為3×8=24個。這樣的配置可以保證存儲單元在數據傳輸高峰情況下能夠同時提供兩組256 bit的數據給執行核，與此同時接收一組從執行內核送過來的256 bit的數據。BWDSP100與ADI公司TS201設計性能比較如表1所示。不考慮讀寫帶寬的前提下，BWDSP100的性能是TS201的8倍。

表1 BWDSP100與TS201的性能比較

2 算法選擇

對于FFT算法來說，其數學表達雖然一致，但具體實現架構卻千差萬別。針對硬件特點選擇合適的處理算法會讓實際程序在處理性能、實用性等方面得到較大的提升。以ADI公司的TigerSHARC系列芯片為例，雖然TS201與TS101保持了指令向下兼容，但是同樣的1 024點FFT，TS101的運行周期數為9 872，TS201的運行周期數卻為12 170。在針對TS201的內部DRAM存儲新特性改用SingLeton算法優化后，ADI提供的新程序1 024點FFT運行周期數為9 419，運算性能得到了明顯提升。

因此，為了滿足設計指標，必須針對BWDSP100結構特點對FFT算法架構進行選擇，這樣才可能貼近工程實際需要，充分發揮芯片特點，給出真正反映芯片性能的測試結果。

2.1 頻域分解與時域分解

FFT算法可以分為頻域分解(DIF)與時域分解(DIT)兩種結構形式。兩種結構的區別在于：頻域分解是在兩點DFT之后有旋轉因子，時域分解是在兩點DFT的節點之前有旋轉因子。很顯然，如果在第一級蝶形運算中加入適當的旋轉因子，時域分解法的乘法系數可以完成對輸入數據實部虛部的不等加權；如果在最后一級蝶形運算中加入適當的旋轉因子，頻域分解法的乘法系數可以實現對輸出數據實部虛部的不等加權。從工程上看，時域分解法的不等加權既可以實現對前方數字化通道不一致的補償，也可以實現對輸入數據取共軛從而方便復用FFT程序得到IFFT的共軛結果。此時只要對輸出結果取共軛就可以得到正常IFFT結果。相較于另外編寫一套IFFT函數，復用FFT程序可以進一步減小程序體積，因此，基于BWDSP100的FFT標準函數采用時域分解FFT算法實現。

2.2 原位算法與非原位算法

原位(in-place)FFT運算是指在不附加緩沖的情況下，FFT的輸入及輸出結果可以在同一地址空間。非原位(not-in-place)FFT運算則指在不附加緩沖的情況下，輸入及輸出結果必須在不同地址空間。很明顯，相較于非原位FFT運算，原位FFT占用的存儲空間更少。

原位FFT必須通過位反序才能獲得正確結果。2的n階FFT需要進行n位的位反序。記某數s的二進制表現形式為(Sn-1Sn-2…S0)，則其對應的逆序數S為(S0S1…Sn-1)。相應的s稱為原序數。將按原序排列的某個矢量改為按逆序排列的過程稱為位反序操作。顯然，對于連續原序數系列來說，其對應的逆序數序列必定非連續。在原位FFT實現中，如果單獨增加一級位反序處理必然會降低程序效率，最好的辦法是將位反序操作與蝶形運算結合起來，在運算中完成位反序。非原位FFT可以通過調整運算架構避免位反序操作。這種調整的代價是要么輸入地址的數據在運算中被改寫，要么必須增加額外緩沖。

位反序對程序效率的影響在于兩點：計算逆序數導致的計算延遲和逆序本身引發取數不連續所導致的讀寫速度下降。現代DSP(包括BWDSP100芯片)大都具有逆序取數功能，逆序數的計算延遲可忽略不計；逆序循環可以保證取數的局部連續性，改善逆序讀寫帶寬。逆序循環的實質在于提取s的最低m個比特與最高m個比特單獨進行逆序交換。由于最低m個比特的順序標號具有連續性，表現在存儲上即具有地址臨近的特點，所以采用逆序循環算法后可以連續讀寫以充分利用芯片帶寬。同時由于逆序數與原序數間互為逆序，因此這種部分逆序交換的結果必然具有封閉性(即循環性)[6-7]。如果處理芯片內部寄存器足夠，可以對循環體涉及到的數據進行暫存和重組，從而實現原位FFT。實際編程中為提升效率而引入了多級合并，此時BWDSP100芯片內部寄存器數量不足以支持這種暫存和重組操作。基于效率優先的考慮，最終選擇了基于逆序循環的帶位反序操作的非原位算法。在不增加額外緩沖的前提下，輸入地址的數據在運算后也可以保留下來。其算法架構如圖2所示。

圖2 8點時域分解原位FFT算法架構

3 算法實現

3.1 運算能力與IO能力的匹配

如前所述，BWDSP100芯片具有非常強大的浮點運算功能，其單核單周期內能夠完成8次浮點加減法及4次浮點乘法運算。與之相比，BWDSP100芯片的數據吞吐能力略顯不足，其單核單周期可以讀取兩個浮點復數(64 bit×2)，存儲一個浮點復數(32 bit×1)。

FFT基本運算為形如(a+bj)±(c+dj)×(Wr+Wij)的蝶形運算。單個蝶形運算的運算量是4次浮點乘法和6次浮點加減法，對應芯片I/O能力的要求是讀取3個浮點型復數(兩個運算數、一個旋轉因子)，存儲兩個浮點型復數(兩個運算結果)。與BWDSP100芯片能力相比，BWDSP100可以在單周期內完成單個蝶形計算，但是無法在單周期內完成相應的IO操作。為了充分發揮芯片的運算性能，必須進行多階合并。如果進行兩階合并的話，兩個周期內要求芯片的I/O能力降低為：讀取4個復數(兩個運算數、兩個旋轉因子)，存儲2個復數(兩個運算結果)，BWDSP100芯片可以滿足這一要求。由于此時乘法運算能力已經全部應用完畢，更多階數的合并并不能帶來更大的得益，反而造成程序復雜度增加。因此，兩階合并是比較合理的選擇。

兩階合并帶來的問題是FFT長度被強制要求為2的偶數次方，這在很多應用場合使用十分不便。觀察圖2可以發現，時域分解FFT算法的前兩級均為特殊旋轉因子(第一級為全1，第二級為1和-j)，這樣單個蝶形運算只需要進行4次加減法。前后兩級蝶形運算合并后共8次加減法，即使考慮到不等加權帶來的4次乘法，兩級運算一個周期即可完成。于是再次遇到了數據吞吐能力不足的問題。繼續延用前面所提的多階合并辦法解決這個問題。多階合并的最大限制由芯片運算核所擁有的內部存儲器個數給出。由于BWDSP100芯片單核擁有64個內部存儲器，合并階數最高可以達到四階。采用判斷處理的辦法。對于偶數階，采用四階合并；對于奇數階，采用三階合并。這樣既解除了FFT長度必須為2的偶數次方限制，又提高了程序的運行效率。

3.2 充分發揮IO帶寬

2.2節已經指出采用逆序循環可以充分發揮芯片的IO帶寬。本節將對其具體實現作出說明。

在第1節已經指出，BWDSP100內部數據存儲器分為3組，每組由8個32 bit獨立存儲器組成，合并為一組256 bit數據提供給4個執行核，每個核從中可以獲得64 bit數據。8個32 bit獨立存儲器地址分布如圖3所示。

圖3 8個32 bit獨立存儲器地址分布圖

如果拼成256 bit數據的幾個地址落到了同一個存儲器就會發生地址沖突。此時，程序流水會被打斷，只有當沖突的數據按周期全部讀寫至各獨立存儲器的緩存后，程序才會繼續執行。

流水打斷會對讀寫帶寬產生極大的影響，為了解決這個問題，在程序編制中應該盡量保證一組所取數據存放在不同存儲器中。BWDSP100的位反序取數指令如下所示[6]：

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

{x,y,z,t}Rs+1:s=br(C)[Vn+=Vm,Vk]

{x,y,z,t}Rs+1:s=br(C)[Wn+=Wm,Wk]

指令中的br前綴代表位反序尋址，前綴br后的立即數C代表位反序的位數，它決定了指令中的哪些位元需要反序。例如指令br(7)[U0+=U1,U2]=R1:0，反序控制參數值為7，表明將地址的[7:0]位進行反序操作。它產生的4 對地址是：[反序(U0)]、[反序(U0)+1]；[反序(U0)+2×U2]、[反序(U0)+2×U2+1]；[反序(U0)+2×2×U2]、[反序(U0)+2×2×U2+1]；[反序(U0)+3×2×U2]、[反序(U0)+3×2×U2+1]。4對地址所取得的數據依次存放到4個執行核X/Y/Z/T的R1和R0中去。而基地址寄存器U0 將被修改為[U0+U1] ，用不反序的U0 值進行基地址的自增。

可以看出，若U2設置為1，則指令所涉及到的4對地址連續，可以保證分布在8個獨立存儲器中，不會發生地址沖突。如果限制輸入數據地址為8的邊界(即低3位為000)，再將增量U1設置為8，則連續的4對取數地址可以等效為FFT逆序數的最低2 bit單獨進行逆序交換的逆序循環，如圖4所示。

圖4 最低2 bit單獨逆序交換取數圖

圖中*號表示任意取值。可以看出，將s的最低2 bit轉化為地址基址U0后，采用4組地址發生器即可實現對所有逆序地址的連續內存訪問，完全規避了地址沖突的風險。完成核間數據交換后，就可以得到正確順序的處理結果。根據逆序循環原理，該處理結果同樣占據輸出地址空間的4段連續內存，保證了輸出存儲時的數據帶寬。需要注意的是，由于復數占據了2個地址，因此實際逆序地址計算中逆序位數為FFT階次n加2。要求輸入FFT數據地址的第n+1位必須為0以保證取數首址的正確。

在編程實現時發現，由于在這種三級/四級合并的編程中，已經盡可能地發揮了BWDSP100的處理能力，再加入核間數據交換操作，其單個周期的并發指令數、ALU資源等都將成為瓶頸。而如果為此增加額外指令周期則影響了程序效率。

仔細研究圖2所示算法架構，可以看出：在后面的處理中，除最后兩級外，其余各級順序讀入的4個運算核的數據剛好分處在4個不同的FFT組內。各組旋轉因子完全相同，運算步驟完全相同，編程十分方便。最后兩級中可以進行高兩位和低兩位的逆序交換，從而完成整個逆序操作。由于最后兩級BWDSP100的運算任務并不飽滿(芯片支持單周期完成8次加減法實際只需要6次加減法)，這種交換操作可以在不影響效率的前提下插入原程序中運行，問題得到了圓滿的解決。由于需要最后兩級完成交換操作，程序能夠實現的最小FFT點數為32點，基本可以滿足工程實現需要。

3.3 提高程序的通用性

FFT運算所需要的位反序位數由FFT的長度決定。設FFT長度N為2的n次階，則FFT運算所需的位反序位數為n，如果考慮到復數占據了2個地址，則BWDSP100位反序取數指令中位反序的位數需要設置為n+2。前面指出，BWDSP100的位反序取數指令中位反序的位數為固定常數，這就意味著一種FFT位反序只能處理對應點數的FFT。如果在實際應用中需要多種長度的FFT，就需要配備多種位反序操作。這樣既增大了程序體積也降低了程序效率。

記原序數s的二進制表現形式為(Sn-1Sn-2…S0)，如果對其尾部添加m個零，記為s*(Sn-1Sn-2…S00m-1…0100)，則其對應的逆序數S*為(0001…0m-1S0S1…Sn-1)。可以看出，對新的原序數s*進行n+mbit位反序后得到的新的逆序數S*與原來的逆序數S數值相等。采用這樣的辦法，可以將原本位數不同的位反轉操作統一擴充到固定點數，對應需要修改的是原序數的遞增量。而BWDSP100指令中該遞增量是可變的。設新的固定位反轉位數為n+m，則新的遞增量由原來的20改變為2m。BWDSP100芯片內部數據地址空間分配如圖5所示。

圖5 BWDSP100內部數據地址空間分配

從圖5可以看出，BWDSP100內部數據Cache分為了3塊，每塊256K字。從地址角度看共18 bit。對于位反序取數指令{x,y,z,t}Rs+1:s=br(C)[Un+=Um,Uk]，Uk取1，固定位反序位數C應該取17。考慮到一個復數數據需占有2個字地址空間，Um的值應該為2(17-r)，式中r為FFT階數。

由于芯片內部有X/Y/Z/T 4個運算核，所取數據將分落到這4個核中，為避免在前面的四級合并中因逆序而進行數據交換，考慮保留前后各2 bit地址到最后交換。將此2 bit考慮后，Um的值為2(19-r)。

這樣修改后，新的程序可以適應各種2的整數次冪FFT處理。同時，由于對原序數尾部添0，逆序數對應頭部填充了多個0，而BWDSP100指令中直接將計算得到的逆序數作為取數地址，因此帶來的限制是輸入FFT數據必須放在3個內部數據塊的某一首地址。在實際工程應用中，這個限制條件對于可以進行絕對地址操作的BWDSP100而言比較容易實現。

4 測試結果

程序編制完成后，在BWDSP100評估板上進行了基于硬件的實際性能測試，與TS201芯片比較結果如表2所示。

表2 BWDSP100與TS201 復數FFT周期數比較

從表2可以看出，BWDSP100的1 024點FFT性能達到TS201的8倍，滿足原設計指標。2 048點FFT由于程序前端只進行了三階合并，性能有所下降。從4 096點開始，BWDSP100采用SRAM的優勢開始體現，較TS201性能有了較大的提升。

5 結束語

基于BWDSP100硬件特點，本文有針對性地給出了FFT算法實現中遇到的問題及解決辦法。與國際成熟的DSP相比，最終FFT性能達到了設計理論要求。盡管本文基于特定平臺，但文中給出的方法對BWDSP100應用及其他芯片FFT實現均有一定的參考作用。

[1] 袁琪,楊康,周建江,等. 大點數FFT算法C6678多核DSP的并行實現[J]. 電子測量技術, 2015, 38(2):74-80.

[2] 張杰,顧乃杰,張明. 龍芯3B處理器上FFT算法向量化研究[J]. 小型微型計算機系統, 2015, 36(7):1639-1643.

[3] 史鴻聲,穆文爭,劉麗. 基于“魂芯一號”的雷達信號處理機設計[J]. 雷達科學與技術, 2012, 10(2):161-164.

SHI Hongsheng, MU Wenzheng, LIU Li. Design of Radar Signal Processor Based on HunXin-1 Chip[J]. Radar Science and Technology, 2012, 10(2):161-164.(in Chinese)

[4] 許德剛. 基于BWDSP100處理器的無源雷達信號處理系統[J]. 艦船電子對抗, 2015, 38(2):72-75.

[5] 賈光帥,洪一,劉小明,等. 基于“魂芯一號”的自適應截位浮點乘法實現[J]. 雷達科學與技術, 2015, 13(3):324-327.

JIA Guangshuai, HONG Yi, LIU Xiaoming, et al. Implementation of Adaptive Truncation of Floating-Point Multiplier Based on BWDSP100 Processor[J]. Radar Science and Technology, 2015, 13(3):324-327.(in Chinese)

[6] 方志紅,張長耀. 逆序循環在原位FFT中的應用[C]∥全國第21屆計算機技術與應用學術會議, 合肥:中國科學技術大學出版社, 2010:43-46.

[7] 方志紅,張長耀,俞根苗. 利用逆序循環實現FFT運算中倒序算法的優化[J]. 信號處理, 2004, 20(5):533-535.