楊 琳，鐘 升，張家田

(1.西安石油大學(xué) 光電油氣測(cè)井與檢測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710065；2.西安微電子研究所，陜西 西安 710065)

FFT的數(shù)據(jù)并行計(jì)算方法研究

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

(1.西安石油大學(xué) 光電油氣測(cè)井與檢測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710065；2.西安微電子研究所，陜西 西安 710065)

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

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

0 引 言

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

1 FFT的數(shù)據(jù)并行計(jì)算

FFT是采用多次蝶式變換和位序反轉(zhuǎn)來實(shí)現(xiàn)DFT變換的快速計(jì)算方法[1-4]。為了方便起見，首先給出8點(diǎn)的FFT蝶式變換架構(gòu)和相應(yīng)的數(shù)據(jù)并行計(jì)算公式，此基礎(chǔ)上再給出N點(diǎn)的FFT蝶式變換的數(shù)據(jù)并行計(jì)算公式。

圖1 8點(diǎn)的FFT蝶式計(jì)算架構(gòu)

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

分裂步1:

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

(1)

計(jì)算步1:

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

(2)

(3)

分裂步2:

(4)

計(jì)算步2:

(5)

(6)

分裂步3:

(7)

計(jì)算步3:

(8)

(9)

表1 N點(diǎn)FFT蝶式變換對(duì)應(yīng)的分裂步與計(jì)算步的數(shù)據(jù)并行計(jì)算公式

至此，F(xiàn)FT蝶式變換的數(shù)據(jù)并行計(jì)算，就由L個(gè)分裂步和L個(gè)計(jì)算步給出。下面將主要討論如何在SIMD PE陣列上實(shí)現(xiàn)FFT蝶式變換?；赟IMD PE陣列的蝶式變換數(shù)據(jù)并行實(shí)現(xiàn)，就是在SIMD PE陣列中實(shí)現(xiàn)L個(gè)分裂步和L個(gè)計(jì)算步。

2 基于SIMD PE陣列的FFT蝶式變換的數(shù)據(jù)并行實(shí)現(xiàn)方法

為簡(jiǎn)化起見，對(duì)8個(gè)通道以信號(hào)長(zhǎng)度均為8的信號(hào)并行進(jìn)行FFT蝶式變換為典型示例，揭示FFT變換是如何在特定處理元陣列中實(shí)現(xiàn)數(shù)據(jù)級(jí)并行處理的。假定處理元陣列的大小與待處理數(shù)據(jù)的規(guī)模一致，即也為8×8，且數(shù)據(jù)編號(hào)與陣列位置一致，即Signal[i][j]是放在PE陣列的處理元PE[i][j]中。陣列中各PE單元之間的互連關(guān)系采用mesh結(jié)構(gòu)，如圖2所示，PE間僅具備相鄰?fù)ㄐ拍芰Α?/p>

對(duì)多個(gè)通道進(jìn)行FFT蝶式變換，可看成是對(duì)一幅二維信號(hào)圖進(jìn)行FFT變換。在SIMD PE陣列上，就是對(duì)存放像素值的所有PE，按行方向?qū)Ω鱌E中的像素值進(jìn)行數(shù)據(jù)并行計(jì)算，計(jì)算后的結(jié)果依然存入相應(yīng)的PE單元中。至此，一幅經(jīng)FFT變換后的譜上的圖像就存放于PE陣列中。以下先以在8×8的PE陣列中對(duì)所有的行同步執(zhí)行8點(diǎn)FFT蝶式變換為例，具體說明FFT變換在SIMD PE陣列中的數(shù)據(jù)并行實(shí)現(xiàn)方法，對(duì)所有列的處理同理。在此基礎(chǔ)上給出N點(diǎn)FFT蝶式變換的SIMD PE陣列實(shí)現(xiàn)方法。

圖2 SIMD PE陣列體系結(jié)構(gòu)(LS-MPP)

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

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

//在計(jì)算階段計(jì)算步1的并行實(shí)現(xiàn)

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;

//在計(jì)算階段計(jì)算步2的并行實(shí)現(xiàn)

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;

//在計(jì)算階段計(jì)算步3的并行實(shí)現(xiàn)

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;

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

3 性能分析與仿真驗(yàn)證

在運(yùn)算量確定的情況下，PE間的通信開銷成為衡量性能的主要指標(biāo)。以下針對(duì)上述實(shí)現(xiàn)方法，對(duì)計(jì)算量和通信量進(jìn)行討論。對(duì)于FFT變換，計(jì)算量主要取決于算法實(shí)現(xiàn)中的乘、加次數(shù)和判別運(yùn)算的次數(shù)。通信量主要取決于蝶式運(yùn)算和位序倒置變換中的數(shù)據(jù)交互次數(shù)，表2給出了相關(guān)的統(tǒng)計(jì)值。

表2 FFT的計(jì)算量及通信次數(shù)

在FFT變換運(yùn)算中，涉及大量的復(fù)數(shù)加法和乘法運(yùn)算。1次復(fù)數(shù)加法由2次實(shí)數(shù)加法構(gòu)成，1次復(fù)數(shù)乘法由4次實(shí)數(shù)乘法和2次實(shí)數(shù)加法構(gòu)成。以N=64×64點(diǎn)的FFT計(jì)算為例，在單處理器環(huán)境下，完成蝶式變換，共需64×(64×6)次復(fù)數(shù)加法和64×(32×5)次復(fù)數(shù)乘法；完成位序倒置，共需64×56次數(shù)據(jù)交換操作。為簡(jiǎn)化起見，約定實(shí)數(shù)的加法和乘法運(yùn)算、位判別運(yùn)算以及PE單元間的相鄰?fù)ㄐ啪稍谝粋€(gè)指令周期內(nèi)完成，并且訪存時(shí)間也為2個(gè)指令執(zhí)行周期。則對(duì)于單處理器而言，按上述約定，大約需要372 736個(gè)指令周期。對(duì)于文中方案，由表2可簡(jiǎn)單地計(jì)算出共需634個(gè)指令周期。其加速比約為587.92，很可觀。文中采用的仿真軟件為Parallaxis-Ⅲ，完成PE陣列的規(guī)模和互聯(lián)結(jié)構(gòu)的定義，其中PE間的數(shù)據(jù)交互是使用系統(tǒng)函數(shù)MOVE模擬實(shí)現(xiàn)。

如圖3所示，通過仿真驗(yàn)證了實(shí)現(xiàn)方法的正確性。

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

4 結(jié)束語

針對(duì)G級(jí)像幀實(shí)時(shí)處理的需要，研究了FFT在MPP計(jì)算機(jī)處理SIMD PE陣列上的數(shù)據(jù)并行計(jì)算實(shí)現(xiàn)方法。該方法的并行度不受算法限制，而且省去了尋址計(jì)算量和尋址時(shí)間，各個(gè)并行操作步驟規(guī)則、簡(jiǎn)潔且通用性強(qiáng)。該方法亦可應(yīng)用于其他算法的并行實(shí)現(xiàn)[12-14]。

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

[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] 陳國良.并行算法的設(shè)計(jì)與分析[M].北京：高等教育出版社，2002:380-409.

[11] 沈緒榜,劉澤響,王 茹.計(jì)算機(jī)體系結(jié)構(gòu)的統(tǒng)一模型[J].計(jì)算機(jī)學(xué)報(bào)，2007，30(5)：729-736.

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

[13] 鐘 升.基于SIMD PE陣列的DCT數(shù)據(jù)并行實(shí)現(xiàn)方法研究[J].電子學(xué)報(bào),2009,37(7):1546-1553.

[14] 鐘 升.基于小波變換的圖像bit糾錯(cuò)數(shù)據(jù)并行實(shí)現(xiàn)研究[J].系統(tǒng)仿真學(xué)報(bào),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">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2017-07-19

陜西省科技統(tǒng)籌創(chuàng)新工程計(jì)劃(2015KTTSGY04-05)

楊 琳(1982-)，女，碩士研究生，研究方向?yàn)殡娮优c通信工程；鐘 升，研究員，研究方向?yàn)橛?jì)算機(jī)體系結(jié)構(gòu)、信號(hào)處理、微工藝系統(tǒng)；張家田，教授，研究方向?yàn)橥ㄐ殴こ膛c信號(hào)處理。

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

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