Bluestein アルゴリズムは何のためにあるのですか？

任意長 N の DFT を 2 のべき乗の FFT に帰着するための手法。kn = (k² + n² − (k − n)²) / 2 という分解を使い、DFT を chirp 乗算と通常の畳み込みの組み合わせに書き直す。畳み込みを 2N − 1 以上の 2 のべき乗 FFT で実行できるため、N が素数や大きな素因数を含む場合でも O(N log N) で計算できる。

Rader アルゴリズムとは何ですか？

N が素数のとき乗法群 Z/NZ* が巡回群であることを使って、長さ N の DFT を長さ N−1 の循環畳み込みに変換する手法。Bluestein と比べて N−1 がさらに合成数のときに効率が良いが、N−1 の素因数構造に依存する。

実数 FFT の出力はなぜ N/2 + 1 ビンで済むのですか？

実数入力 x[n] の DFT X[k] は共役対称性 X[N − k] = conj(X[k]) を満たすため、後半 N/2 ビンは前半から復元できる。独立な値は X[0] から X[N/2] までの N/2 + 1 個 (両端は実数で残りは複素) であり、メモリ・演算量とも約半分に削減できる。

packing 法とは何ですか？

二つの長さ N 実数列を一つの長さ N 複素数列 z[n] = x[n] + j y[n] に詰めてから 1 回の複素 FFT を実行し、共役対称性から X[k] と Y[k] を分離する手法。実数 FFT 一回分の複素 FFT で二本の信号を変換できるためバッチ処理で有用である。

Bluestein の数値精度は通常の FFT より悪いですか？

Chirp 因子 e^(iπn²/N) を高精度に評価する必要があるため、素朴な実装では Cooley-Tukey より相対誤差がやや大きくなる。n² の評価で先に N を法とした整数演算を行い、その後に三角関数を評価することで誤差を抑えられる。sangi は事前計算した chirp テーブルを使い、内部 FFT サイズも 2N 以上の 2 のべきにパディングして循環畳み込みの汚染を回避する。

Bluestein 任意長 FFT と実数 FFT

概要

Cooley-Tukey FFT は $N$ が小さい素因数の積であれば $O(N \log N)$ で計算できる。しかし $N$ が大きな素数や扱いにくい素因数を含む場合には混合基数の枠組みから外れる。本記事は任意長 $N$ に対応する Bluestein 法・Rader 法・素因数アルゴリズム (PFA) と、実数入力で計算量・メモリを半減する RealFFT の仕組みを解説する。

sangi の API は RealFFT クラスおよび FFT クラスとして公開されている。

Bluestein (chirp z 変換)

1968 年に Bluestein が発表した手法。 DFT の指数 $kn$ を次のように分解する:

$$kn = \frac{k^2 + n^2 - (k - n)^2}{2}$$

これを DFT 定義に代入すると、$X[k]$ は chirp 乗算と一次元畳み込みの組み合わせに書き直せる:

$$X[k] = \omega_N^{-k^2/2} \sum_{n=0}^{N-1} \bigl( x[n] \, \omega_N^{-n^2/2} \bigr) \, \omega_N^{(k-n)^2/2}$$

$a[n] = x[n] \, \omega_N^{-n^2/2}$, $b[m] = \omega_N^{m^2/2}$ と置けば、 DFT は $a$ と $b$ の (線形) 畳み込みを長さ $\geq 2N - 1$ の循環畳み込みにパディングして計算するだけで得られる。

循環畳み込みは 2 のべきの FFT で $O(N \log N)$ に計算できるため、 Bluestein 全体の計算量は $O(N \log N)$ となる。$N$ がいかなる整数であっても適用可能であることが最大の特長である。

$\omega_N^{m^2/2} = e^{-i\pi m^2 / N}$ の指数は時間とともに二次的に変化するため chirp 信号 (周波数が時間に比例して変化する波形) と呼ばれ、本手法は chirp z 変換とも呼ばれる。

関連記事: Bluestein FFT (チャープ z 変換)

Rader アルゴリズム

1968 年に Rader が発表した、$N$ が素数のときに有効な手法。 $\mathbb{Z}/N\mathbb{Z}^*$ が位数 $N-1$ の巡回群であることに着目し、原始根 $g$ を用いて DFT のインデックスを並べ替えると、長さ $N$ の DFT は長さ $N-1$ の循環畳み込みに変換される。

$N-1$ がさらに小さい素因数の積 (5-smooth など) なら混合基数 FFT で畳み込みが効率的に評価できる。 Bluestein と比べるとパディングが不要で内部 FFT サイズが小さく、$N-1$ の素因数分解が良性のときに有利である。

素因数アルゴリズム (PFA / Good-Thomas)

Good (1958) と Thomas (1963) による手法。 $N = N_1 N_2$ かつ $\gcd(N_1, N_2) = 1$ のとき、中国剰余定理を使って 2 次元の独立な軸に並べ替えると、 twiddle 因子が不要な単純な 2 次元 DFT に分解できる:

$$X[k_1, k_2] = \sum_{n_1=0}^{N_1-1} \sum_{n_2=0}^{N_2-1} x[n_1, n_2] \, \omega_{N_1}^{k_1 n_1} \, \omega_{N_2}^{k_2 n_2}$$

Cooley-Tukey と異なり段間の twiddle 乗算がないため、小さな素因数同士の合成数で乗算回数を削減できる。ただし入出力のインデックス再配置 (CRT 写像) が必要で、$\gcd = 1$ の条件が制限となる。現代の SIMD 実装では twiddle のコストが小さいため、PFA は特定のサイズで限定的に使われる。

実数入力の共役対称性

実数列 $x[n] \in \mathbb{R}$ の DFT $X[k]$ は次の共役対称性を満たす:

$$X[N - k] = \overline{X[k]}, \qquad k = 1, \dots, N-1$$

つまり負周波数ビンは正周波数ビンの複素共役であり、独立な情報は半分しかない。 $X[0]$ と (偶数 $N$ のとき) $X[N/2]$ は実数となり、残り $N/2 - 1$ 個が複素数として独立である。合計の自由度は $N$ 個の実数となり、入力の自由度と一致する。

RealFFT — N → N/2+1 圧縮

共役対称性を陽に使うと、長さ $N$ 実数 FFT の出力は長さ $N/2 + 1$ の複素数配列で表現できる。計算量・メモリとも約半分に削減される。sangi の RealFFT はこの圧縮表現を直接返す。

実数列を複素 FFT に詰める方法

実数入力を長さ $N/2$ の複素列に変換し、長さ $N/2$ の複素 FFT を 1 回呼び出して後処理で長さ $N$ の結果に変換する手法が標準的である。

入力 $x[0..N-1]$ を偶数・奇数で分け、$z[k] = x[2k] + i \, x[2k+1]$ ($k = 0, \dots, N/2 - 1$) を構成する。
長さ $N/2$ の複素 FFT $Z = \mathrm{FFT}(z)$ を計算する。
後処理で $X[k]$ を $Z[k]$ と $Z[N/2 - k]$ の線形結合として復元する: $$X[k] = \tfrac{1}{2}(Z[k] + \overline{Z[N/2-k]}) - \tfrac{i}{2} \omega_N^k (Z[k] - \overline{Z[N/2-k]})$$

$X[0]$ と $X[N/2]$ は実数として個別に処理される。この手法により実数 FFT は同サイズの複素 FFT の約半分の計算量で実現される。

逆変換

逆変換は同じ packing を逆向きに辿る。共役対称性を満たさない入力を渡すと結果が実数にならないため、sangi の RealFFT::ifft は引数サイズが output_size / 2 + 1 と一致することを検証する。

packing 法 (二信号同時 FFT)

二本の実数列 $x[n]$ と $y[n]$ を $z[n] = x[n] + i \, y[n]$ にまとめ、一度の複素 FFT $Z[k] = \mathrm{FFT}(z)$ を実行する。共役対称性から両者のスペクトルを分離できる:

$$X[k] = \tfrac{1}{2}\bigl(Z[k] + \overline{Z[N-k]}\bigr), \qquad Y[k] = \tfrac{1}{2i}\bigl(Z[k] - \overline{Z[N-k]}\bigr)$$

計算は実数 FFT 2 本分よりわずかに重いが、複素 FFT カーネルが 1 つで済むため、ステレオ音声処理など 2 チャンネル同期処理で有用である。 sangi の simd_fft::stereoFFT はインターリーブされたステレオ WAV に対しこの分離を 1 回の FFT で実行する。

Hermitian FFT と DCT/DST

実数 FFT の出力は Hermitian な複素数列となる。 DCT (離散コサイン変換) と DST (離散サイン変換) はそれぞれ偶対称・奇対称延長の DFT として書け、長さ $N$ の入力を長さ $2N$ や $4N$ の実数 FFT に帰着できる。

DCT-II: $X[k] = \sum_{n=0}^{N-1} x[n] \cos\!\bigl(\frac{\pi(2n+1)k}{2N}\bigr)$。 Makhoul (1980) の構成では長さ $N$ の実数 FFT 1 回 + $O(N)$ 後処理で計算できる。
DST-II: サイン基底による変換。DCT との対応で同様に長さ $N$ 実数 FFT に帰着可能。
DCT-III / DST-III: 逆変換 (IDCT / IDST)。同じ枠組みで対称。

sangi の dct / idct は $N$ が 5-smooth のときに Makhoul 経路を取り $O(N \log N)$、それ以外では直接定義通り $O(N^2)$ で計算する。 dst / idst は現在のところ直接定義通りの計算である。

関連記事: 重要なフーリエ変換ペア / フーリエ余弦変換 / フーリエ正弦変換

Bluestein の数値精度

Bluestein は chirp 因子 $\omega_N^{n^2/2} = e^{-i\pi n^2 / N}$ を多用する。 $n^2$ を直接 double で評価すると桁落ちが発生し、相対誤差が Cooley-Tukey 比で 1〜2 桁悪化することがある。以下の対策がある。

整数モジュロ評価: $n^2 \bmod 2N$ を整数演算で先に求め、その後で三角関数を呼ぶ。 $\pi n^2 / N$ の位相が $[0, 2\pi)$ に収まり、引数縮約の誤差が抑えられる。
chirp テーブルの事前計算: $b[m] = \omega_N^{m^2/2}$ を全長分テーブル化し、繰り返し計算を避ける。
内部 FFT サイズの padding: 循環畳み込みでの巻き込みを避けるため、$M \geq 2N - 1$ の 2 のべき乗を選ぶ。

これらを実装した Bluestein は通常の Cooley-Tukey FFT と相対誤差で同オーダーになる (Rabiner et al. 1969, Smith 1995 のベンチマーク参照)。

sangi の動作

sangi の FFT::fft(data) は次のように分岐する:

FFT::is_valid_size(N) が true ($N$ が 5-smooth) — Cooley-Tukey 混合基数
それ以外 — Bluestein を呼び出し、内部で長さ $M = $ next_valid_size(2N - 1) の 5-smooth FFT を 3 回実行

Bluestein 経路では作業領域 ($M$ の 2 倍程度) が確保されるため、ピークメモリは Cooley-Tukey より大きい。頻繁に呼び出すサイズが固定ならば FFTEngine を構築して twiddle/chirp テーブルをキャッシュすると良い。

実数入力には RealFFT<Real> を使うのが望ましい。内部で長さ $N/2$ の複素 FFT を 1 回 + 後処理で構成されるため、複素 FFT を直接呼ぶより約 2 倍速い。

参考文献

Bluestein, L. I. (1970). "A linear filtering approach to the computation of discrete Fourier transform". IEEE Trans. Audio Electroacoust., 18, 451–455.
Rader, C. M. (1968). "Discrete Fourier transforms when the number of data samples is prime". Proc. IEEE, 56, 1107–1108.
Good, I. J. (1958). "The Interaction Algorithm and Practical Fourier Analysis". J. Royal Stat. Soc. B, 20, 361–372.
Makhoul, J. (1980). "A fast cosine transform in one and two dimensions". IEEE Trans. ASSP, 28, 27–34.
Sorensen, H. V., Jones, D. L., Heideman, M. T., and Burrus, C. S. (1987). "Real-valued fast Fourier transform algorithms". IEEE Trans. ASSP, 35, 849–863.