FDTD 算法的FPGA 實現(xiàn)

徐奇澎，郭裕順

(杭州電子科技大學(xué)電子信息學(xué)院，杭州 310018)

FDTD 法即時域有限差分方法(Finite Difference in Time Domain)，是把Maxwell 方程式在時間和空間域進行差分，通過電場和磁場的交替計算，來模擬電磁波傳播的一種數(shù)值計算方法［1］。自上世紀六十年代出現(xiàn)這種算法以來，人們已對它作了大量研究，是目前電磁波數(shù)值仿真的主要方法，被廣泛地應(yīng)用于科學(xué)與工程中各種電磁場與電磁波問題。

FDTD 算法應(yīng)用中存在的一個主要問題是:當計算區(qū)域較大，即空間差分產(chǎn)生較多網(wǎng)格時，需要耗費大量的計算資源與時間［2］。為此降低計算代價，減少計算時間，一直是重要的研究課題。以往采取的措施不外乎兩個方面:一是改進算法，設(shè)法提高計算效率;二是在通常的計算平臺上采用并行計算。

FPGA 的發(fā)展為解決這類大規(guī)模計算問題提供了一種新途徑。FDTD 計算的特點是在固定的一個空間網(wǎng)格中進行重復(fù)的算術(shù)運算，數(shù)據(jù)流規(guī)則，控制流較為簡單，因此可利用FPGA 可編程的特點，設(shè)計專門的算術(shù)運算電路，實現(xiàn)FDTD 的計算。一塊FPGA 芯片可同時實現(xiàn)多個乘法器、加法器，輕易實現(xiàn)并行運算。這樣針對FDTD 算法專門設(shè)計的運算電路，可大大提高計算速度。最早用FPGA 實現(xiàn)FDTD 算法的是Schneider，他用定點整數(shù)格式實現(xiàn)了FDTD 算法的一維運算，但網(wǎng)格單元僅10個［3］。Gandhi、Durbano 分別用浮點數(shù)實現(xiàn)了FDTD 算法的二維、三維運算，但受限于當時FPGA 的性能，這些設(shè)計工作頻率較低，性能不高，未能體現(xiàn)硬件計算的優(yōu)勢［4］。后來Wu Shuguang 在Gandhi 的研究基礎(chǔ)上把數(shù)據(jù)存儲位置由主機系統(tǒng)內(nèi)存改為FPGA 片上SRAM，大大減少了數(shù)據(jù)傳輸耗費的時間［5］。最近Hideki Kawaguchi 等采用分配式并行存儲使數(shù)據(jù)讀寫速度進一步加快，當FPGA 工作頻率在50 MHz時，計算速度比軟件工作在2.8 GHz 的個人電腦上快了50倍［6－7］。因此可以預(yù)計隨著FPGA 的進一步發(fā)展，用FPGA 實現(xiàn)FDTD 算法的優(yōu)勢將更加明顯。

1 二維FDTD 算法

本文考慮二維FDTD 算法的實現(xiàn)。二維情形下所有電磁量與Z 坐標無關(guān)［8］，即?/?z=0。這時電磁場的直角分量可以分為獨立的兩組，即以Hx、Hy、Ez為一組的TM 波，及以Ex、Ey、Hz為一組的TE波［9］。以TM 波為例，二維TM 波的Yee 元胞如圖1所示，電場值根據(jù)周圍4個磁場值完成更新，磁場值根據(jù)周圍2個電場值完成更新，更新關(guān)系滿足右手螺旋法則。

圖1 二維TM 波的Yee 元胞

為了便于硬件實現(xiàn)，對上述關(guān)系中的1/2 按照一定規(guī)則做出處理［10］。處理后TM 波的計算公式如下:

下面介紹用FPGA 來完成上述計算的方法。

2 算法的硬件實現(xiàn)

用硬件實現(xiàn)式(1)～式(3)的計算首先要確定數(shù)字的表示方法。本文采用符合IEEE-754 32 位標準的單精度浮點數(shù)格式，即將一個數(shù)表示為V=(－1)S×1.M×2E-BIAS，其中最高位S 是符號位，隨后的8 位E是指數(shù)位，后面的23 位M 是尾數(shù)位，BIAS=28－1=127。

圖2 電路結(jié)構(gòu)

圖3 計算電路

實現(xiàn)一次FDTD 計算的工作流程如圖4所示，分為以下幾步:

(1)初始化:確定要計算的網(wǎng)格大小，所要計算的時間步，計算式系數(shù)。

(2)加入激勵后，整個計算電路開始工作。本文中因為激勵信號是添加在RAM Ez中，所以先計算Ez值。根據(jù)計數(shù)器計算出相應(yīng)的地址值，從RAM 中取得相鄰兩個單元的磁場值，Ez的計算電路完成一個點的Ez值計算，E值計數(shù)器會加1，然后檢查其是否數(shù)到最大值。如果計數(shù)器沒有達到最大值則再次計算Ez值，進行RAM Ez中下一個點的更新。如果是計數(shù)器達到最大值則表明RAM Ez中所有的點都完成一個時間步的更新，轉(zhuǎn)換更新算式開始進行磁場值的更新。

(3)Hx和Hy的更新同時進行，當一個點的Hx和Hy值更新完成之后，H值計數(shù)器會加1，然后檢查其是否數(shù)到最大值。如果計數(shù)器沒有達到最大值則再次計算Hx和Hy值，進行RAM Hx和RAM Hy中下一個點的更新。如果計數(shù)器達到最大值則表明RAM Hx和RAM Hy中所有的點都完成一個時間步的更新，同時也意味著網(wǎng)格中所有的點都完成了一個時間步的電場值和磁場值的更新。

(4)完成一個時間步更新后，時間步計數(shù)器T加1，然后檢查T 是否小于nstep，如果是則表示還未完成設(shè)定的時間步，回到算法第3 步，開始網(wǎng)格內(nèi)所有節(jié)點下一個時間步的更新。如果不是則表示計算完成設(shè)定的時間步，整個算法完成。

圖4 算法流程

在設(shè)計過程中我們可以根據(jù)并行原理［11－12］和FPGA 架構(gòu)特點，設(shè)計多個計算電路，讓這些計算電路同時運算，可以大大減少運算所需要的時間，提高算法的效率。

3 實驗結(jié)果

根據(jù)上述方案，用Verilog HDL 編寫了代碼，以ALTERA 公司的EP3C55F484C6為目標器件進行了編譯，用ModelSim 進行了仿真。設(shè)置計算網(wǎng)格大小為60 ×60，在網(wǎng)格中間加一個高斯脈沖g(t)=exp(－0.5(t－t0)2/δ2)充當激勵源，觀察3 600個網(wǎng)格點經(jīng)過60個步長時間的運算結(jié)果和運算所需要的時間。圖5 就是仿真所得部分結(jié)果，端口data_out 輸出的是EZ RAM 中的值，圖中可以看到輸出的值是符合IEEE 754 標準的32 位單精度浮點數(shù)，由于采用了流水線設(shè)計，每個時鐘輸出一個運算結(jié)果。圖6 是3 600個點在經(jīng)過60 步后仿真結(jié)果。

為了比較FPGA 計算的結(jié)果，我們對同樣計算用Visual C++6.0 編寫了軟件程序，結(jié)果如圖7，兩者誤差不超過0.11%。FPGA 設(shè)計的最高工作頻率可達118.39 MHz，完成整個運算過程需670 000個時鐘周期。假設(shè)FPGA 工作頻率為10 MHz，則完成全部運算需要0.067 s，而軟件程序在主頻2.2 GHz的個人電腦上運行時間是1.48 s，因此FPGA 計算要快幾乎22倍。當采用2倍并行計算時，完成全部運算需368 310個時鐘周期，加速比可達40倍。因此FPGA 計算對提高FDTD 算法速度的效果是十分明顯的。

圖5 ModelSim 仿真結(jié)果

圖6 FPGA 運算結(jié)果

圖7 C 語言運算結(jié)果

4 結(jié)論

本文對FDTD 算法的FPGA 實現(xiàn)進行了初步研究。以兩維情形為例，編寫了Verilog HDL 實現(xiàn)代碼，在Altera FPGA 上進行了編譯仿真，獲得了初步試驗結(jié)果。在設(shè)計中采用了流水線技術(shù)、并行運算等手段，同時采用了FPGA 片內(nèi)RAM 作為存儲單元，取得了較好的加速效果。隨著制造工藝的發(fā)展，F(xiàn)PGA 片內(nèi)資源會越來越豐富，基于FPGA 的計算可能成為提高FDTD 算法性能的一條有效途徑。

［1］Yee K S.Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equation in Isotropic Media［J］.IEEE Trans Antennas Propagate，1966，(14):302－307.

［2］丁偉.時域有限差分法關(guān)鍵技術(shù)及其應(yīng)用研究［D］.西安電子科技大學(xué)，2007.

［3］Schneider R N，Turner L E，Okoniewski M M.Application of FPGA Technology to Accelerate the Finite-Difference Time-Domain(FDTD)Method［C］//The 2002 ACM International Symposium on Field-Programmable Gate Arrays(FPGA'02)，2002:24－26.

［4］Durbano J P，Ortiz F E，Humphrey J R.Hardware Implementation of a Three-Dimensional Finite-Difference Time-Domain Algorithm［C］//IEEE Antennas and Wireless Propagation Letters，2003:2，54－57.

［5］Wu S.An FPGA Implementation of FDTD Codes for Reconfigurable High Performance Computing［D］.M.S.Thesis，Dept.of Electrical and Computer Engineering and Computer Science，University of Cincinnati，2005.

［6］Kawaguchi H，F(xiàn)ujita Y，F(xiàn)ujishima Y.Improved Architecture of FDTD/FIT Dedicated Computer for Higher Performance Computation［J］.IEEE Transactions on Magnetics，2008，44(6):1226－1229.

［7］Fujita Y，Kawaguchi H.Full Custom PCB Implementation of FDTD/FIT Dedicated Computer［J］.IEEE Transactions on Magnetics，2009，45(3):1100－1103.

［8］張煊，張明德.二維光子晶體波導(dǎo)及耦合器的特性模擬［J］.電子器件，2005，28(6):63－67.

［9］葛德彪，閏玉波.電磁波時域有限差分方法［M］.西安電子科技大學(xué)出版社，2002.

［10］Dennis M Sullivan.Electromagnetic Simulation Using the FDTD Method［M］.IEEE Microwave Theory and Techniques Society，2000.

［11］雷繼兆.PC和服務(wù)器集群下的并行FDTD 算法及其應(yīng)用研究［D］.西安電子科技大學(xué)，2009.

［12］劉瑜.FDTD 算法的網(wǎng)絡(luò)并行研究及其電磁應(yīng)用［D］.電子科技大學(xué)，2008.