Cooley-Tukey FFT — 混合基数アルゴリズム

概要

離散フーリエ変換 (DFT) を $O(N \log N)$ で計算するアルゴリズム群を総称して FFT (高速フーリエ変換) と呼ぶ。 その中で最も広く使われるのが 1965 年に Cooley と Tukey が再発見した分割統治法 (Cooley-Tukey FFT) である。 sangi の FFT クラスはこの系統のアルゴリズムを核とし、 $N$ の合成性に応じて radix-2 SIMD バタフライ・混合基数バタフライ・Bluestein を自動的に切り替える。

API の使い方については FFT API リファレンス を参照。

DFT の定義と計算量

長さ $N$ の複素数列 $x[0], \dots, x[N-1]$ に対する離散フーリエ変換は次で定義される。

$$X[k] = \sum_{n=0}^{N-1} x[n] \, \omega_N^{kn}, \qquad \omega_N = e^{-2\pi i / N}, \qquad k = 0, \dots, N-1$$

定義通り計算すると各 $k$ で $N$ 回の複素乗算と $N-1$ 回の複素加算を要し、 全体では $\Theta(N^2)$ の演算量となる。$N = 1024$ で約 100 万回、$N = 2^{20}$ で約 $10^{12}$ 回に達し、 $N$ が増えると実用に耐えない。

Cooley-Tukey は $N = N_1 N_2$ という分解を用い、長さ $N$ の DFT を 長さ $N_1$ の DFT $N_2$ 個と長さ $N_2$ の DFT $N_1$ 個に置き換える。 これを再帰的に適用すると深さ $\log N$、各段の作業量 $O(N)$ となり、 全体の計算量が $O(N \log N)$ に削減される。

radix-2 Cooley-Tukey

$N = 2^m$ のとき、入力を偶数インデックスと奇数インデックスに分割する:

$$X[k] = \sum_{n=0}^{N/2-1} x[2n] \, \omega_N^{2nk} + \omega_N^{k} \sum_{n=0}^{N/2-1} x[2n+1] \, \omega_N^{2nk}$$

$\omega_N^{2} = \omega_{N/2}$ であるから、これは長さ $N/2$ の DFT 二つの線形結合となる:

$$X[k] = E[k] + \omega_N^{k} \, O[k], \qquad X[k+N/2] = E[k] - \omega_N^{k} \, O[k]$$

$E, O$ はそれぞれ偶数項・奇数項の長さ $N/2$ の DFT。 この再帰関係はバタフライ演算と呼ばれる対称的な複素加減算で表現される。

$$\begin{pmatrix} X[k] \\ X[k+N/2] \end{pmatrix} = \begin{pmatrix} 1 & \omega_N^k \\ 1 & -\omega_N^k \end{pmatrix} \begin{pmatrix} E[k] \\ O[k] \end{pmatrix}$$

1 段あたり $N$ 回の複素加減算と $N/2$ 回の複素乗算、$\log_2 N$ 段で 合計約 $\frac{N}{2} \log_2 N$ 回の複素乗算となる。

DIT と DIF

radix-2 Cooley-Tukey には双対な 2 つの定式化がある。

• DIT (Decimation-In-Time): 時間領域インデックス $n$ を偶奇で分割する。入力をビット反転順に並べ替えてからバタフライを進め、出力は自然順となる。
• DIF (Decimation-In-Frequency): 周波数領域インデックス $k$ を $k$ と $k+N/2$ に分割する。自然順入力にバタフライを進め、出力をビット反転で並べ替える。

算術演算の回数は両者で同一であるが、メモリアクセスパターン・twiddle 因子の出現順・SIMD 適合性が異なる。 実装上は順変換に DIT、逆変換に DIF を選ぶことで、共通の twiddle テーブルを順方向で読めるようにする手法もある。

ビット反転置換

DIT では入力配列をビット反転順に並べ替えてから本体ループに入る。 $N = 2^m$ のとき、インデックス $n = b_{m-1} b_{m-2} \cdots b_1 b_0$ (2 進表記) を $\tilde{n} = b_0 b_1 \cdots b_{m-1}$ に置換する。

$$N = 8 \text{ のときの例: } 0,1,2,3,4,5,6,7 \;\;\to\;\; 0,4,2,6,1,5,3,7$$

素朴な実装ではビット単位のループで $O(N \log N)$ となるが、 Gold-Rader アルゴリズムや段階的なビット交換で $O(N)$ に圧縮できる。 また、DIT-DIF を順逆で組み合わせるとビット反転置換が両端で打ち消し合い、 畳み込みなど往復で使う場面では明示的な置換を省ける場合がある。

radix-4 と radix-8

$N = 4^m$ のとき 4 分割を行うと、1 段あたりのバタフライは 4 点の複素線形変換となる:

$$\begin{pmatrix} X_0 \\ X_1 \\ X_2 \\ X_3 \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -j & -1 & j \\ 1 & -1 & 1 & -1 \\ 1 & j & -1 & -j \end{pmatrix} \begin{pmatrix} E_0 \\ \omega^{k} E_1 \\ \omega^{2k} E_2 \\ \omega^{3k} E_3 \end{pmatrix}$$

係数 $\pm 1, \pm j$ による乗算はビット反転と符号反転のみで実装でき、 実乗算回数が radix-2 比で約 25% 削減される。 radix-8 ではさらに $\omega_8 = e^{-i\pi/4}$ が現れ、$\pm 1, \pm j, \pm(1\pm j)/\sqrt{2}$ で構成される。

段数は $\log_4 N = \tfrac{1}{2} \log_2 N$、$\log_8 N = \tfrac{1}{3} \log_2 N$ となるため、 メモリアクセスとループオーバーヘッドが減る点でも有利である。

split-radix

Yavne (1968) および Duhamel, Hollmann (1984) による手法。 偶数項を radix-2、奇数項を radix-4 で分解する非対称な split を用いる:

$$X[k] = E[k] + \omega_N^k Z_1[k] + \omega_N^{3k} Z_3[k]$$

$E$ は長さ $N/2$ の DFT、$Z_1, Z_3$ は長さ $N/4$ の DFT。 実数乗算回数が $\frac{4}{3} N \log_2 N - \tfrac{38}{9} N + \cdots$ となり、 pure radix-2 ($2 N \log_2 N$) より約 33% 少ない。 2004 年に Johnson と Frigo が同じ係数を保つよりコンパクトな形式を発表した。

計算量の小ささに対して実装が複雑で、ループ構造が radix-2 SIMD と整合しにくい。 sangi は AVX2 で float なら 4 複素数 (8-wide)、double なら 2 複素数 (4-wide) を同時に扱う radix-2 SIMD バタフライを採用し、 理論的乗算回数より実機スループットを優先している。

混合基数 (radix-2,3,5)

$N$ が 2 の冪でなくても、$N$ が小さい素因数の積であれば Cooley-Tukey の枠組みが使える。 sangi は $N = 2^a \cdot 3^b \cdot 5^c$ (5-smooth) を直接処理する混合基数経路を持つ。

• radix-3 バタフライ: $\omega_3 = e^{-2\pi i / 3} = -\tfrac{1}{2} - \tfrac{\sqrt{3}}{2} i$ を用いる 3 点 DFT。実乗算 4 回・実加算 12 回。
• radix-5 バタフライ: $\omega_5$ を用いる 5 点 DFT。実乗算 10 回・実加算 34 回 (Winograd の最小乗算形式)。

混合基数の段の順序は、まず大きな素因数から先に処理して短い段ほど SIMD で並列化しやすくする戦略 (Singleton 1967, Bergland 1968) がある。 sangi は 2 のべきの段を AVX2 で集約し、3 と 5 の段はスカラーで処理する単純化を取る。

FFT::is_valid_size(n) は $n$ が 5-smooth かを判定し、 FFT::next_valid_size(n) は混合基数経路で扱える最小の $\geq n$ を返す。

twiddle 因子の精度

各バタフライで掛ける $\omega_N^k = \cos(2\pi k / N) - i \sin(2\pi k / N)$ を twiddle 因子と呼ぶ。 実装上は次の 3 通りがある。

• 直接計算: std::cos, std::sin を毎回呼ぶ。最も正確だが遅い。
• 事前テーブル: $\{\omega_N^0, \dots, \omega_N^{N-1}\}$ をプログラム開始時に計算してキャッシュ。$O(N)$ メモリで $O(1)$ アクセス。
• 漸化式: $w_{k+1} = w_k \cdot w_1$ を保持し逐次更新。$O(1)$ メモリだが誤差が累積する。

漸化式の相対誤差は段ごとに $O(\varepsilon)$ ずつ蓄積し、長さ $N = 2^m$ の変換では $O(m \varepsilon)$ に達する ($\varepsilon$ は丸め単位)。大きな $N$ では Tang (1989) の半角公式や Buneman の再帰式で誤差を抑える工夫が必要となる。 sangi は事前テーブル方式を基本とし、テーブル境界の値のみ cos/sin で補正する設計である。

in-place と stride 変換

radix-2 バタフライは 2 要素を読み 2 要素を書く一対の操作であり、入力配列上で in-place に進められる。 sangi の FFT::fft(std::span<T>) は入力スパンを直接書き換える in-place 変換である。

多次元 FFT (FFT-multidim) では行方向・列方向それぞれの軸に対して 1D FFT を適用するが、 列方向では要素のストライドが $N_x$ となり連続アクセスが崩れる。 実装上は次のいずれかを選ぶ:

• stride 対応バタフライ: ポインタ算術にストライドを直接組み込む。コードは単純だがキャッシュ効率に劣る。
• 転置 + 連続 FFT: 列方向を一度連続メモリに転置してから 1D FFT を適用する。ブロック転置でキャッシュをフレンドリーに保てる。

sangi の fft2d / fft3d は転置経路を採用し、各 1D FFT は連続メモリ上で SIMD radix-2 を走らせる。

sangi の dispatch

FFT<T, Real>::fft(data) はサイズ $N = $ data.size() に応じて以下のように振り分ける。

逆変換は順変換に対し $1/N$ 正規化を施す。 FFTEngine クラスを使うと FFTDivMode で正規化規約 (None / Forward / Inverse / Both) を選べる。

Bluestein 経路・任意長アルゴリズムの詳細は FFT-bluestein-real、 実数 FFT の packing は同記事の後半、NTT の詳細は FFT-ntt-bigint、 多次元と畳み込みは FFT-multidim を参照。

参考文献

• Cooley, J. W. and Tukey, J. W. (1965). "An Algorithm for the Machine Calculation of Complex Fourier Series". Mathematics of Computation, 19, 297–301.
• Singleton, R. C. (1969). "An Algorithm for Computing the Mixed Radix Fast Fourier Transform". IEEE Trans. Audio Electroacoust., 17, 93–103.
• Duhamel, P. and Hollmann, H. (1984). "Split-radix FFT algorithm". Electronics Letters, 20, 14–16.
• Johnson, S. G. and Frigo, M. (2007). "A Modified Split-Radix FFT With Fewer Arithmetic Operations". IEEE Trans. Signal Process., 55, 111–119.
• Van Loan, C. (1992). Computational Frameworks for the Fast Fourier Transform. SIAM.