基于FFTW庫分步傅里葉變換算法并行方案研究

劉 帥， 李 智， 王 晶

（1.裝備學院 研究生管理大隊，北京101416； 2.裝備學院 重點實驗室，北京101416； 3.91336部隊）

在求解拋物型波動方程時，Hardin和Tappert于1973年提出了一種分步傅里葉變換算法［1］，該算法是一種采用FFT技術的頻域算法，其計算不受波長限制，在求解大區域電波傳播問題中得到了廣泛的應用。該算法屬于一種步進迭代算法，當計算區域采樣點數或迭代次數增加時，計算速度便會下降。

FFTW是由麻省理工學院（MIT）的Frigo和Johnson基于C語言開發的程序庫［2］，能夠高效地進行任意維數的離散傅里葉變換（DFT）、離散正弦變換（DST）以及離散余弦變換（DCT），是目前公認的最有效的變換工具。它能夠根據用戶的底層硬件屬性如緩存、內存、寄存器大小等，自動調整算法參數，以獲得最優方案。FFTW2.1.5和最新的FFTW3.3Alpha測試版還能夠支持共享存儲多線程并行和分布式存儲MPI（消息傳遞接口）并行［3］。

為解決適應波動方程計算量大、計算速度慢的問題，作者借助于FFTW庫研究了SSFT算法的并行方案，測試結果表明，所設計方案是有效的。

1 SSFT算法介紹

為闡述方便，先對SSFT算法的計算過程做簡要的介紹。無論是三維還是二維情況，采用Feit－Fleck近 似 的 拋 物 型 波 動 方 程［4－6］都 可 以 化為如下形式

這樣拋物方程的解可以近似為

稱B為折射算子，A為繞射算子。在SSFT算法中，折射部分的計算在空間域內進行，繞射部分，也就是exp（ΔxA）u部分的計算在變換域即頻域內進行，最終結果為

在一個計算步長Δx內，計算被分解成4個步驟：

1）計算初始場分布的傅氏變換U0＝F（u）；

2）繞射步計算T1＝exp（ΔxA）F（u）；

3）對T1計 算 其 傅 氏 反 變 換T2＝F－1（exp（ΔxA）F（u））；……