FFT的數據并行計算方法研究

楊 琳，鐘 升，張家田

(1.西安石油大學 光電油氣測井與檢測教育部重點實驗室，陜西 西安 710065；2.西安微電子研究所，陜西 西安 710065)

FFT的數據并行計算方法研究

楊 琳1，鐘 升2，張家田1

(1.西安石油大學 光電油氣測井與檢測教育部重點實驗室，陜西 西安 710065；2.西安微電子研究所，陜西 西安 710065)

為滿足G(Gigabytes)級像素幀的實時性處理需求，針對信號處理系統中處理計算量大、實時性要求高的特點，剖析了解算過程內在的數據并行特性，深入研究了基于計算陣列的譜圖解算數據并行算法。提出了一種基于MPP(Massively Parallel Processor)計算機SIMD PE陣列的FFT的數據并行計算實現方法。首先根據FFT架構中的數據交互一致性，給出了數據并行計算的表達式。提出一種基于PE標識，進行條件操作的SIMD PE陣列數據并行實現方法。該方法不但省去了并行處理中的數據尋址時間開銷，而且使得數據并行操作更為規則、簡潔，滿足了陣列操作規則性強的處理要求，大幅度地提高了MPP計算機并行計算處理速度。該方案是一種簡潔有效的PE自治問題解決方案，以更合理的方法和更高的效率實現了常規經典算法，在數據并行計算領域中，無疑具有重要的理論意義和應用價值，將在嵌入式信號處理中發揮愈來愈重要的作用。

快速傅里葉變換;SIMD PE陣列;映射語言；MPP計算機

0 引 言

在數字圖像變換處理中，DFT(Discrete Fourier Transform)是重要的變換算法，被廣泛地用于圖像匹配、紋理分析、去噪濾波等圖像處理領域。FFT是采用蝶式交換的處理方式，實現了DFT變換的快速算法。針對G級像素幀的FFT處理，由于具有較高的實時性要求，使得基于MPP計算機處理陣列的數據并行實現方式成為研究熱點[1-4]。在并行處理計算機上如何結合體系結構高效率地實現典型算法是并行計算的核心技術[5-10]。基于陣列的數據并行實現就是針對確定的并行算法，基于數據在陣列中的位置分布特性，擬定合適的并行計算處理步驟。在各個并行步驟已確定的前提下，各個并行步中，被操作數在陣列中的存儲位置、傳送距離及處理方式也是確定、已知的。對各并行步中存儲被操作數PE單元進行事先標識，并將標識碼預置于相應PE的PSR寄存器[11-14]，使得基于PE標識的數據并行操作，不但具有了PE“自治”(automations)能力，而且省去了數據尋址時間。各并行步的數據處理是規范，符合算法的陣列實現特點，能夠很好地滿足MPP計算機處理陣列操作規則性強的要求。

1 FFT的數據并行計算

FFT是采用多次蝶式變換和位序反轉來實現DFT變換的快速計算方法[1-4]。為了方便起見，首先給出8點的FFT蝶式變換架構和相應的數據并行計算公式，此基礎上再給出N點的FFT蝶式變換的數據并行計算公式。

圖1 8點的FFT蝶式計算架構

從圖1可直觀看出，各次蝶式變換中的數據交互和相應的計算處理是規則的，采用不同基色(帶填充色的點和無填充色)的點來標識，在各次蝶式變換中的數據交互總是發生在以相同間隔分布的兩類基色點間，且無填充色的點只進行“+”運算，填充色點進行“-”和與變換系數的乘法運算。兩類代表不同操作數據類型的填充色原點的分布以及它們間的間隔距離，在每次蝶式解算中均呈現顯著規律，數據并行操作便是基于這些規律來實現的。為便于描述，基于不同的運算操作來標識不同的數據類型，參照提升小波變換的命名方式，簡稱為數據分裂步，將具體的數據傳送和計算操作稱為數據計算步。這樣一來，8點的蝶式變換就劃分為分裂1、分裂2和分裂3以及計算步1、計算步2和計算步3。式(1)～(9)采用向量操作方式，給出對應于分裂步和數據計算步的并行計算公式。

分裂步1:

A0:={a(0),a(1),a(2),a(3)};B0:={a(4),a(5),a(6),a(7)}

(1)

計算步1:

A1:={s1(0),s1(1),s1(2),s1(3)}=A0+B0

(2)

(3)

分裂步2:

(4)

計算步2:

(5)

(6)

分裂步3:

(7)

計算步3:

(8)

(9)

表1 N點FFT蝶式變換對應的分裂步與計算步的數據并行計算公式

至此，FFT蝶式變換的數據并行計算，就由L個分裂步和L個計算步給出。下面將主要討論如何在SIMD PE陣列上實現FFT蝶式變換。基于SIMD PE陣列的蝶式變換數據并行實現，就是在SIMD PE陣列中實現L個分裂步和L個計算步。

2 基于SIMD PE陣列的FFT蝶式變換的數據并行實現方法

為簡化起見，對8個通道以信號長度均為8的信號并行進行FFT蝶式變換為典型示例，揭示FFT變換是如何在特定處理元陣列中實現數據級并行處理的。假定處理元陣列的大小與待處理數據的規模一致，即也為8×8，且數據編號與陣列位置一致，即Signal[i][j]是放在PE陣列的處理元PE[i][j]中。陣列中各PE單元之間的互連關系采用mesh結構，如圖2所示，PE間僅具備相鄰通信能力。

對多個通道進行FFT蝶式變換，可看成是對一幅二維信號圖進行FFT變換。在SIMD PE陣列上，就是對存放像素值的所有PE，按行方向對各PE中的像素值進行數據并行計算，計算后的結果依然存入相應的PE單元中。至此，一幅經FFT變換后的譜上的圖像就存放于PE陣列中。以下先以在8×8的PE陣列中對所有的行同步執行8點FFT蝶式變換為例，具體說明FFT變換在SIMD PE陣列中的數據并行實現方法，對所有列的處理同理。在此基礎上給出N點FFT蝶式變換的SIMD PE陣列實現方法。

圖2 SIMD PE陣列體系結構(LS-MPP)

在SIMD PE陣列中，指令執行將結合PE單元具體的坐標位置進行相應的條件操作。由于各個計算步中各個分量的序號是與像素值無關的，是依據表1所給各個分裂步事先確定的，或者說，可以在執行計算操作前就可以確定。因此，不采用求解坐標位置的辦法，而是利用各個PE單元內的Program Status Register (PSR)，即程序狀態寄存器，并利用PSR中相應的狀態位來提前標識各個分裂步的結果，或者說執行分裂步的工作。在8×8的PE陣列中的各個PE的PSR狀態位中相應bit位將按如下原則設置：對于列坐標j=0,1,2,3的PE單元，其PSR寄存器的b0位置為0，其他PE單元內PSR寄存器的b0位設置為1；對于列坐標j=0,1,4,5的PE單元，其PSR寄存器的b1位設置為0，其他PE單元內PSR寄存器的b1位設置為1；對于列坐標j=0,2,4,6的PE單元，其PSR寄存器的b2位設置為0，其他PE單元內PSR寄存器的b2位設置為1。這樣一來就可通過判別PSR的狀態位來確定各個計算步中相應向量的各個分裂，或執行指定運算操作的各個PE[i,j]單元，該方法使得后續的陣列操作簡潔、規則，而且沒有了尋址時間開銷。此外，FFT變換系數依據表1中的計算公式確定并作為立即數被提前預置于對應位置的PE內的立即數寄存器中。在SIMD PE陣列中，各個計算步的并行實現，就是基于各個PSR設置，按表1給定的計算公式，在陣列中執行相應的條件傳送與條件計算操作。各個計算步的并行執行是以狀態位為條件，采用映射語言M(Mapping language/Middle language)[12]的條件傳送語句與條件算術語句等描述的。

//在數據加載階段，初始數據A0與B0已加載至各PE單元的DATA_registeR0寄存器中，變換系數加載于相應PE單元的立即數寄存器LITERAL_register_1與LITERAL_registe_2中，DATA_registeR1

//在計算階段計算步1的并行實現

1 If DATA_registeR1=(b0= =0)? DATA_registeR0[j+4]:NOP;//PE[j]←PE[j+4]；

2 If DATA_registeR1=(b0= =1)? DATA_registeR0[j-4]:NOP;//PE[j]←PE[j-4]；

3 If DATA_registeR0=(b0= =0)? {DATA_registeR0+DATA_registeR1}:NOP;

4 If DATA_registeR0=(b0= =1)? {DATA_registeR1-DATA_registeR0}:NOP;

5 If DATA_registeR0=(b0= =1)? {DATA_registeR0×LITERAL_register_1}:NOP;

//在計算階段計算步2的并行實現

1 If DATA_registeR1=(b0= =0)? DATA_registeR0[j+2]:NOP;//PE[j]←PE[j+2]；

2 If DATA_registeR1=(b0= =1)? DATA_registeR0[j-2]:NOP;//PE[j]←PE[j-2]；

3 If DATA_registeR0=(b0= =0)? {DATA_registeR0+DATA_registeR1}:NOP;

4 If DATA_registeR0=(b0= =1)? {DATA_registeR1-DATA_registeR0}:NOP;

5 If DATA_registeR0=(b0= =1)? {DATA_registeR0×LITERAL_register_2}:NOP;

//在計算階段計算步3的并行實現

1 If DATA_registeR1=(b0= =0)? DATA_registeR0[j+1]:NOP;//PE[j]←PE[j+1]；

2 If DATA_registeR1=(b0= =1)? DATA_registeR0[j-1]:NOP;//PE[j]←PE[j-1]；

3 If DATA_registeR0=(b0= =0)? {DATA_registeR0+DATA_registeR1}:NOP;

4 If DATA_registeR0=(b0= =1)? {DATA_registeR1-DATA_registeR0}:NOP;

對應于IDFT的IFFT蝶式處理算法，由于各處理階段中的待處理數據在陣列中的分布位置，傳送距離及處理方式與FFT基本一致，唯一不同的僅僅是各計算步中的變換系數而已，所以僅需將FFT中各計算步的變換系數更換為其倒數即可。例如，將ωk更換為ω-k，并存儲于相同的PE單元中。其他設置及處理均與FFT一致，在此不再贅述。

3 性能分析與仿真驗證

在運算量確定的情況下，PE間的通信開銷成為衡量性能的主要指標。以下針對上述實現方法，對計算量和通信量進行討論。對于FFT變換，計算量主要取決于算法實現中的乘、加次數和判別運算的次數。通信量主要取決于蝶式運算和位序倒置變換中的數據交互次數，表2給出了相關的統計值。

表2 FFT的計算量及通信次數

在FFT變換運算中，涉及大量的復數加法和乘法運算。1次復數加法由2次實數加法構成，1次復數乘法由4次實數乘法和2次實數加法構成。以N=64×64點的FFT計算為例，在單處理器環境下，完成蝶式變換，共需64×(64×6)次復數加法和64×(32×5)次復數乘法；完成位序倒置，共需64×56次數據交換操作。為簡化起見，約定實數的加法和乘法運算、位判別運算以及PE單元間的相鄰通信均可在一個指令周期內完成，并且訪存時間也為2個指令執行周期。則對于單處理器而言，按上述約定，大約需要372 736個指令周期。對于文中方案，由表2可簡單地計算出共需634個指令周期。其加速比約為587.92，很可觀。文中采用的仿真軟件為Parallaxis-Ⅲ，完成PE陣列的規模和互聯結構的定義，其中PE間的數據交互是使用系統函數MOVE模擬實現。

如圖3所示，通過仿真驗證了實現方法的正確性。

圖3 2D-FFT變換仿真圖像

4 結束語

針對G級像幀實時處理的需要，研究了FFT在MPP計算機處理SIMD PE陣列上的數據并行計算實現方法。該方法的并行度不受算法限制，而且省去了尋址計算量和尋址時間，各個并行操作步驟規則、簡潔且通用性強。該方法亦可應用于其他算法的并行實現[12-14]。

從硬件上講，并行度僅受PE陣列規模的限制。由于SIMD PE陣列具有較強的可剪裁性，其大小是容易隨應用的并行性要求而變化的；特別是隨著芯片集成度的提高，陣列規模急速上升，容易滿足大規模數據的處理需求。

[1] Tong C,Swartztrauber P N.Ordered fast Fourier transform on an massively parallel hyper cube multiprocessor[J].Journal of Parallel and Distribute Computing,1991,12(1):50-59.

[2] Swartztrauber P N.Multiprocessor FFTs[J].Parallel Computing,1987,5(1-2):197-210.

[3] Jamieson L H,Muller L P,Siegal H.FFT algorithm for SIMD parallel processing system[J].Journal of Parallel and Distributed Computing,1986,3(1):48-71.

[4] Berland G D,Wilson D.A fast Fourier transform algorithm for a global,highly parallel processor[J].IEEE Transactions on Audio and Electroacoustics,1969,17(2):125-127.

[5] Khailany B,Dally W J,Kapasi U J,et al.Imagine:media processing with streams[J].IEEE Micro,2001,21(2):35-46.

[6] Rixner S. Stream processor architecture[M].Boston,MA:Kluwer Academic Publishers,2002.

[7] Kapasi U J.The imagine stream processor[C]//Proceedings of IEEE 2002 international conference on computer design.[s.l.]:IEEE,2002:282-288.

[8] Serebrin B.A stream processor development platform[C]//Proceedings of IEEE 2002 international conference on computer design.[s.l.]:IEEE,2002:303-308.

[9] Crowley P. Network processor design issues and practices[M].[s.l.]:Morgan Kaufmann Publishers,2003.

[10] 陳國良.并行算法的設計與分析[M].北京：高等教育出版社，2002:380-409.

[11] 沈緒榜,劉澤響,王 茹.計算機體系結構的統一模型[J].計算機學報，2007，30(5)：729-736.

[12] 鐘 升,沈緒榜,鄭江濱,等.提升小波變換的數據并行計算方法研究[J].計算機學報,2011,34(7):1323-1331.

[13] 鐘 升.基于SIMD PE陣列的DCT數據并行實現方法研究[J].電子學報,2009,37(7):1546-1553.

[14] 鐘 升.基于小波變換的圖像bit糾錯數據并行實現研究[J].系統仿真學報,2008,20(8):2137-2141.

ResearchonDataParallelComputationMethodofFFT

YANG Lin1，ZHONG Sheng2，ZHANG Jia-tian1

(1.Key Laboratory of Photo-electricity Gas-oil Logging and Detecting of Ministry of Education, Xi’an Shiyou University,Xi’an 710065,China; 2.Xi’an Microelectronic Technique Institute,Xi’an 710065,China)

In order to satisfy the real-time processing requirement of G-level pixel frame,considering the intensive and real time computing requirement of signal processing in embedded signal system,the inner data parallelism of the calculation process is analyzed,and the data parallel algorithms of spectrogram calculation based on computing array is also researched.A data parallel computation method implemented on SIMD PE array of MPP (Massively Parallel Processor) computer for FFT transform is presented.Based on the data computing consistency of FFT,the expression of data parallel computing is given firstly.Then a method of data parallel computing based on SIMD PE array to execute conditional operations by using of PE identifier is proposed,which not only omits the time cost of addressing,but makes the data parallel operation more regular and compact (only in computation statements and move statements).It meets the features of high regularity required by SIMD and greatly improves MPP computer processing speed,which is also a simple and effective PE autonomy solution,realizing conventional classic algorithms with more rational method and higher efficiency.It has important theoretical significance and application value in the area of data parallel undoubtedly,which will play a more and more important role in embedded signal processing.

FFT;SIMD PE array;mapping language；MPP computer

TP301

A

1673-629X(2017)10-0091-05

2016-10-27

2017-02-14 < class="emphasis_bold">網絡出版時間

時間：2017-07-19

陜西省科技統籌創新工程計劃(2015KTTSGY04-05)

楊 琳(1982-)，女，碩士研究生，研究方向為電子與通信工程；鐘 升，研究員，研究方向為計算機體系結構、信號處理、微工藝系統；張家田，教授，研究方向為通信工程與信號處理。

http://kns.cnki.net/kcms/detail/61.1450.TP.20170719.1109.028.html

10.3969/j.issn.1673-629X.2017.10.020