NEC

MathKeisanユーザーズガイド

FFTおよびPARFFT

概要

FFTは高速フーリエ変換を行うライブラリであり、離散フーリエ変換(DFT)を計算するアルゴリズムです。 $ n $ 個の複素数データ $ x_{j} $ に対するDFT $ y_{k} $ は以下の式で与えられます。

\[ y_k = \sum_{j=0}^{n-1} x_j e^{\frac{-2 \pi i}{n} jk} \; \; \; \; \; \; k = 0, 1, \ldots, n-1 \]

ここで、 $ i^{2} = -1 $ です。逆変換は以下の式で与えられます。

\[ x_j = \frac{1}{n} \sum_{k=0}^{n-1} y_k e^{\frac{ 2 \pi i }{n} jk} \; \; \; \; \; j = 0, 1, \ldots, n-1 \]

上記の式では、順変換に対しスケーリングを行なわず、逆変換を $ 1/n $ でスケーリングします。他に、 $ 1/n $ で順変換をスケーリングし逆変換にスケーリングを行わない方法、もしくは両方の変換を $ 1/\sqrt{n} $ でスケーリングする方法があります。

Itanium® Processor Family (IPF)版MathKeisanのFFTは、 Intel® Math Kernel Library(MKL)のFFTおよびDFTのインタフェースを備えています。 本章ではこれ以降、SX版のMathKeisanのFFTについてのみ記載しています。

ユーザインタフェース

FFTの機能はパッケージによってそれぞれ大きく異なっており、標準となるユーザインタフェースがありません。
SX版MathKeisanのFFTは、 Cray®のLibSciTM 3.1 および HP MLIBのVECLIB [3] と同じ機能を備えています。 MathKeisanのLibSciTMと Cray®のLibSciTM のインタフェースはワーク配列のサイズ以外は同じです。 FFT のリファレンス(man page)にワーク配列サイズの計算方法を記載しています。

MathKeisan FFTは、以下の変換に対応しています。

これらの変換はすべて、以下のデータ型のマッピングに利用できます。

ユーザインタフェース情報は、いくつかの箇所に記載されています。

PARFFT

PARFFTは、FFTのOpenMP版です。 FFTと同じユーザインタフェースを備えているので、 FFTにリンクされたコードであればPARFFTにリンクさせることができます。
環境変数 OMP_NUM_THREADS が np に設定されている場合、PARFFTは np 個のスレッド上で実行されます。OMP_NUM_THREADSが設定されていない場合には、 PARFFTは mp 個のスレッド上で実行されます。この時、 mp はリソースグループまたはノード内の最大スレッド数です。

FFTアルゴリズムの概略

$ \omega $ を

\[ \omega^k_n = e^{\frac{-2 \pi i}{n}k} = \cos \frac{2 \pi }{n}k - i \sin\frac{2 \pi }{n}k \]

と定義するとDFTは以下のような行列とベクトルの積により与えられます。

\[ \left[ \begin{array}{lllcl} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0 \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\ \vdots & \vdots & \vdots & & \vdots \\ \omega_n^0 & \omega_n^{(n-1)} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \\ \end{array} \right] \left[ \begin{array}{c} x_0 \\ x_1 \\ x_3 \\ \vdots \\ x_{n-1} \\ \end{array}\right] = \left[ \begin{array}{c} y_0 \\ y_1 \\ y_3 \\ \vdots \\ y_{n-1} \\ \end{array}\right] \]

$ \omega_{n} $ のべき乗の計算を無視するならば、DFTの浮動小数点演算数は $ 8n^{2} $ になります。1965年に Cooley と Tukey [1] が、DFTを計算するFFTアルゴリズムを開発しました。 $ n $ が2のべき乗である時、FFTアルゴリズムでは演算数が $ 5n \lg n $ と大幅に少なくなります。一般的な $ n $ では、演算数は $ O(n \lg n) $ です。当時CooleyとTukeyは知りませんでしたが、FFTのアルゴリズムはすでに Gauss [2] により1805年に提唱され、その死後1866年に発表されていました。 本アルゴリズムを理解するため、 $ n $ が2のべき乗である場合を考えます。 $ \omega $ をDFTの式に代入します。

\[ y_k = \sum^{n-1}_{j=0} x_j \omega_n^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1 \]

$ m = n/2 $ とし $ x_{j} $ を偶数成分 $ \underline{x}_{j}=x_{2j} $ と奇数成分 $ \overline{x}_{j}=x_{2j+1} $ に分割します。

\[ y_k = \sum^{m-1}_{j=0} \underline{x}_j \omega_n^{2jk} + \sum^{m-1}_{j=0} \overline{x}_j \omega_n^{(2j+1)k} \; \; \; \; \; \; \; k = 0, 1, \ldots, n-1 \]

$ \omega_{n}^{2jk} = \omega_{m}^{ jk} $ という周期性を使って $ y_{k} $ を変換します。

\[ y_k = \sum^{m-1}_{j=0} \underline{x}_j \omega_{m}^{jk} + \omega_n^{k} \sum^{m-1}_{j=0} \overline{x}_j \omega_{m}^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1 \]

$ y_{k} $ を前半の $ m $ 個と後半の $ m $ 個のエントリに分割します。

\begin{eqnarray} y_k & = & \sum^{m-1}_{j=0} \underline{x}_j \omega_{m}^{jk} + \omega_n^{k} \sum^{m-1}_{j=0} \overline{x}_j \omega_{m}^{jk} \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; k = 0, 1, \ldots, m-1 \nonumber \\ y_{k+m} & = & \sum^{m-1}_{j=0} \underline{x}_j \omega_{m}^{j(k+m)} + \omega_n^{k+m} \sum^{m-1}_{j=0} \overline{x}_j \omega_{m}^{j(k+m)} \; \; \; \; \; k = 0, 1, \ldots, m-1 \nonumber \end{eqnarray}

$ \omega_{m}^{jm} = 1 $ と $ \omega_{n}^{m} = -1 $ という周期性を使って $ y_{k+m} $ を変換します。

\begin{eqnarray} y_k & = & \sum^{m-1}_{j=0} \underline{x}_j \omega_{m}^{jk} + \omega_n^{k} \sum^{m-1}_{j=0} \overline{x}_j \omega_{m}^{jk} \; \; \; \; \; k = 0, 1, \ldots, m-1 \nonumber \\ y_{k+m} & = & \sum^{m-1}_{j=0} \underline{x}_j \omega_{m}^{jk} - \omega_n^{k} \sum^{m-1}_{j=0} \overline{x}_j \omega_{m}^{jk} \; \; \; \; \; k = 0, 1, \ldots, m-1 \nonumber \end{eqnarray}

サイズが $ n $ 、演算数が $ 8n^{2} $ のオリジナルDFTを、サイズが $ m=n/2 $ 、演算数が $ 2n^{2} $ の2つのDFTと置き換えました。 $ \lg n $ ステップ続けて、長さ $ n $ のオリジナルDFTを長さ1のDFTと置き換えることができます。 これを基数2のFFTと呼びます。
本アルゴリズムは2のべき乗だけでなく、一般の $ n $ まで拡張することが可能です。 上記の分割処理は、時間間引き形(DIT)Cooley-Tukeyアルゴリズムに対応しています。 合計で8回の分割があります。 MathKeisan FFTライブラリで用いられる分割は、 周波数間引き形(DIF)Stockhamアルゴリズムです。 本アルゴリズムから派生したアルゴリズムは Van Loan [4] に記載されています。 DIF Stockhamアルゴリズムの演算数はCooley-Tukeyと同じですが、演算を行う順序が違います。

多重1次元FFTとベクトル長

FFTの計算は、一連の基数のステップ $ n_{1}, n_{2}, \ldots, n_{r} $ として行われます。 $ n_{1}, n_{2}, \ldots, n_{r} $ は一般的に、 $ n = n_{1} n_{2} \cdots n_{r} $ を素因数分解して得られる小さな素数です。 計算中のベクトル長は基数 $ n_{1} $ ステップで $ n_{1} $ 、基数 $ n_{2} $ ステップで $ n_{1} n_{2} $ 、基数 $ p $ ステップでは以下のようになります。

\[ \prod_{j=1}^{p}n_j \]

$ n $ が小さいと、FFT計算時のベクトル長は短くなります。 アプリケーションによっては、同じ長さの1次元FFTがたくさん必要になります。 FFTの長さが $ n $ で、FFTの数が $ h $ だとすると、多重1次元FFTは以下のように記述することができます。

\[ y_{kl} = \sum_{j=0}^{n-1} x_{jl} \omega_n^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1 \; \; \; \; \; l = 0, 1, \ldots, h-1 \]

ここで $ y_{kl} $ と $ x_{jl} $ は複素数の2次元配列です。 多重1次元FFT計算時の基数 $ n_{p} $ ステップにおけるベクトル長は、以下のようになります。

\[ h \prod_{j=1}^{p}n_j \]

このようにベクトル長が長いため、1次元FFTより多重1次元FFTの方が性能上有利です。

2次元FFTおよび3次元FFT

複素数 $ x_{jl} $ の2次元フーリエ変換 $ y_{kp} $ は、以下の式で与えられます。

\begin{eqnarray} y_{kp} & = & \sum_{l=0}^{h-1} \sum_{j=0}^{n-1} x_{jl} \omega_h^{lp} \omega_n^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1 \; \; \; \; \; p = 0, 1, \ldots, h-1 \nonumber \\ & = & \sum_{l=0}^{h-1} \omega_h^{lp} \sum_{j=0}^{n-1} x_{jl} \omega_n^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1 \; \; \; \; \; p = 0, 1, \ldots, h-1 \nonumber \end{eqnarray}

ここで

\[ z_{kl} = \sum_{j=0}^{n-1} x_{jl} \omega_n^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1 \; \; \; \; \; l = 0, 1, \ldots, h-1 \]

と定義すると $ y_{kp} $ は以下の式になります。

\[ y_{kp} = \sum_{l=0}^{h-1} z_{kl} \omega_h^{lp} \; \; \; \; \; k = 0, 1, \ldots, n-1 \; \; \; \; \; p = 0, 1, \ldots, h-1 \]

こうして2次元フーリエ変換は、 $ z_{kl} $ と $ y_{kp} $ の2つの多重1次元FFTにより置き換えられました。 同様に、3次元フーリエ変換は3つの多重1次元FFTにより置き換えられます。 MathKeisanルーチンの内部では、2次元および3次元フーリエ変換は 多重1次元計算カーネルを呼び出して計算しています。 これをrow-column法といいます。

実数-複素数FFTおよび複素数-実数FFT

実数-複素数FFTは対称性を持っているので、複素数-複素数FFTに比べ、 演算数およびメモリ使用量において約2倍性能が良くなっています。 周期性があることがわかるように、DFTを三角関数 $ \sin $ と $ \cos $ を使って記述します。

\[ y_k=\sum_{j=0}^{n-1} x_j \cos \frac{2\pi }{n} jk - i \sum_{j=0}^{n-1} x_j \sin \frac{2\pi }{n} jk \; \; \; \; \; k = 0, 1, \ldots, n-1 \]

$ x_{j} $ が実数である場合、 $ y_{k} $ は $ k = n/2 $ について周期性を持つ複素共役となります。これにより、 $ y_{k} $ の値の半分しか格納しなくて良いことになります。値 $ y_{0} $ は、DC周波数と呼ばれます。 $ n $ が偶数の場合、値 $ y_{n/2} $ はナイキスト周波数と呼ばれます。恒等式 $ \sin 0 = 0 $ および $ \sin \pi j = 0 $ から、 $ y_{0} $ と $ y_{n/2} $ が実数であることが分かります。 入力ベクトルがサイズ $ n $ の実数である実数-複素数FFTの場合、速度とメモリの優位性を活かすには、 出力ベクトルはサイズ $ n/2+1 $ の複素数である必要があります。 入力ベクトルがサイズ $ n/2+1 $ の複素数である複素数-実数FFTの場合、出力ベクトルはサイズ $ n $ の実数である必要があります。サイズが $ n/2 $ ではなく $ n/2+1 $ である理由は、DC周波数とナイキスト周波数の冗長ゼロ部が格納されるからです。 LibSciTMやVECLIBのインタフェースの中には 空間節約を利用しているものもありますが、そうでないものもあります。 インタフェース情報は、 FFTのリファレンス(man page) を参照してください。

整合寸法の効果

多重1次元、2次元、3次元FFTサブルーチンは、 データを2次元または3次元の配列に格納しています。 配列の整合寸法が2の高次べき乗である場合、 バンクコンフリクトにより性能が劣化します。 配列の整合寸法は奇数に設定してバンクコンフリクトを防止するのが最善です。

FFTルーチン一覧

名称欄の ?  は、接頭辞欄のいずれかに置き換わります。 接頭辞はそれぞれ精度を表します。
凡例： S = REAL(kind=4)、 D = REAL(kind=8)、 C = COMPLEX(kind=4)、 Z = COMPLEX(kind=8)
  名称 接頭辞 説明
LibSci ?FFT CC ZZ 1次元複素数-複素数FFT
?FFT C Z 1次元複素数-複素数FFT
?FFT2 C Z 1次元複素数-複素数FFT
?FFT SC DZ 1次元実数-複素数FFT
?FFT CS ZD 1次元複素数-実数FFT
?FFT2 RC DZ 1次元実数-複素数FFT
?FFT2 CR ZD 1次元複素数-実数FFT
?FFT2D CC ZZ 2次元複素数-複素数FFT
?FFT2D C Z 2次元複素数-複素数FFT
?FFT2D SC DZ 2次元実数-複素数FFT
?FFT2D CS ZD 2次元複素数-実数FFT
?FFT3D CC ZZ 3次元複素数-複素数FFT
?FFT3D C Z 3次元複素数-複素数FFT
?FFT3D SC DZ 3次元実数-複素数FFT
?FFT3D CS ZD 3次元複素数-実数FFT
?FFTM CC ZZ 多重1次元複素数-複素数FFT
M?FFT C Z 多重1次元複素数-複素数FFT
?FFTMLT C Z 多重1次元複素数-複素数FFT
?FFTM SC DZ 多重 1次元実数-複素数FFT
?FFTM CS ZD 多重 1次元複素数-実数FFT
?FFTMLT R D 多重 1次元実数-複素数、複素数-実数FFT
?FTFAX F D 実数のFFT因数を生成、三角関数表を計算する
?FTFAX C Z 複素数のFFT因数を生成、三角関数表を計算する
XERFFT   FFTエラー処理を行う
VECLIB ?1DFFT C Z 1次元複素数-複素数FFT、複素数引数
?1DFFT S D 1次元複素数-複素数FFT、実数引数
?RC1FT C Z 1次元実数-複素数、複素数-実数FFT、複素数引数
?RC1FT S D 1次元実数-複素数、複素数-実数FFT、実数引数
?2DFFT C Z 2次元複素数-複素数FFT、複素数引数
?2DFFT S D 2次元複素数-複素数FFT、実数引数
?RC2FT C Z 2次元実数-複素数、複素数-実数FFT、複素数引数
?RC2FT S D 2次元実数-複素数、複素数-実数FFT、実数引数
?3DFFT C Z 3次元複素数-複素数FFT、複素数引数
?3DFFT S D 3次元複素数-複素数FFT、実数引数
?RC3FT C Z 3次元実数-複素数、複素数-実数FFT、複素数引数
?RC3FT S D 3次元実数-複素数、複素数-実数FFT、実数引数
?FFTS C Z 多重 1次元複素数-複素数FFT、複素数引数
?FFTS S D 多重 1次元複素数-複素数FFT、実数引数
?RCFTS C Z 多重 1次元実数-複素数、複素数-実数FFT、複素数引数
?RCFTS S D 多重 1次元実数-複素数、複素数-実数FFT、実数引数

参考文献

  1. J. Cooley, J. Tukey, "An algorithm for the machine calculation of complex Fourier series," Math. Comp. 19 (1965), 297-301.
  2. C.F. Gauss, "Theoria interpolationis methodo nova tracta," Carl Friedrich Gauss Werke, Band3, Koniglichen Gesellschaft de Wissenschafen, Gotingen, 1866.
  3. MLIB User's Guide, vol 2.
  4. C. Van Loan, Computational Frameworks for the Fast Fourier Transform , SIAM, Philadelphia, 1992.