FFT

DFT (離散フーリエ変換) を高速に実行するアルゴリズムのことを FFT (高速フーリエ変換) と呼ぶ. 離散フーリエ変換とフーリエ変換は似て非なるもので, この言葉はいささか不用意であり, 本来 FDFT と名付けるべきだったが今更である. 離散フーリエ変換という言葉自体, 本当は「フーリエ級数展開」のバリエーションと見なせるものなので, これを「フーリエ変換」と言ってしまうのは乱暴である.

大雑把に言うと, 時系列データに対して DFT を施すと周波数成分が計算できる. よくオーディオ機器にスペクトル表示器が備わっているが, アレをイメージするとよい.

今日は, そのような本来の目的は一旦忘れて, 純粋にアルゴリズムの話を進めたい.

定義にはいくつか流儀があるが, おおむね次のようなものとなる.

X[n] = ΣW(k * n, N) x[k], n,k = 0...N-1.

x[k] は変換元となるデータ列, X[n] は変換結果のデータ列, W(m, N) は回転子と呼ばれるもので, 1 の N 乗根のうち 1 でないものの m 乗を表している. 通常,

W(m, N) = exp(-2πi * m / N)
= cos(-2π * m / N) + i sin(-2π * m / N)

指数 (角度) のマイナス符号は, あまり深い意味はないが, 変換結果 X[n] が exp(2πi * m / N) の係数となるようにするためには必要である. (もともと i と -i は同等の立場なので, こだわる意味はまったくない. 純粋に形式だけの問題.)

サイズ N = 2 の DFT

定義は

X[0] = Σ W(k * 0, 2) x[k], k = 0, 1,
X[1] = Σ W(k * 1, 2) x[k], k = 0, 1

となるので,

X[0] = Σ W(0, 2) x[k] = x[0] + x[1],
X[1] = Σ W(k, 2) x[k] = x[0] - x[1]

ここで, W(1, 2) = exp(-2πi / 2) = exp(-πi) = cos(-π) + i sin(-π) = -1 などを使った.

サイズ 2 の DFT の計算には加減算だけが登場することが分かる.

N が 2 のべきとなっているとき, 例えば N = 8 としよう. X[n], x[k] の添え字 n, k を 2 進数で n = pqr, k = hijと表そう. pqr = 4p + 2q + r, hij = 4h + 2i + j.

DFT の定義は,

X[pqr] = ΣW(hij * pqr, 8) x[hij]

となる.

バタフライとビットリバース

回転子の変形を試みる.

W(hij * pqr, 8)
=W((h00 + ij) * pqr, 8)
=W(h00 * pqr, 8) * W(ij * pqr, 8)
=W(h00 * (pq0 + r), 8) * W(ij * pqr, 8)
=W(hpq000, 8) * W(h00 * r, 8) * W(ij * pqr, 8)
=W(h * r, 2) * W(ij * pqr, 8)

元の DFT は,

X[pqr] = Σ[ij] W(ij * pqr, 8) Σ[h] W(h * r, 2) x[hij]

となるが, Σ[h] の部分は, サイズ 2 の DFT となっている. (Σ[x] は x についての和という意味.)

さらに変形する.

W(ij * pqr, 8)
= W(ij * (pq0 + r), 8)
= W(ij * pq0, 8) * W(ij * r, 8)
= W(ij * pq, 4) * W(ij * r, 8)

なので,

X[pqr]
= Σ[ij] W(ij * pqr, 8) Σ[h] W(h * r, 2) x[hij]
= Σ[ij] W(ij * pq, 4) * W(ij * r, 8) Σ[h] W(h * r, 2) x[hij]

このように変形すると, サイズ 8 の DFT は, サイズ 2 の DFT の結果 Σ[h] に W(ij * r, 8) を掛けたものに対するサイズ 4 の DFT となっていることが分かる.

そしてこのサイズ 4 の DFT もまた同様にサイズ 2 の DFT に分解できる. 最終的には

X[pqr]
= Σ[j] W(j * p, 2) * W(j * qr, 8) * Σ[i] W(i * q, 2) * W(i * r, 4) * Σ[h] W(h * r, 2) * x[hij]

となる. 各 DFT を分けて書くと, 次のようになる.

y[rij] = Σ[h] W(h * r, 2) * x[hij],
z[qrj] = Σ[i] W(i * q, 2) * W(i *  r, 4) * y[rij],
X[pqr] = Σ[j] W(j * p, 2) * W(j * qr, 8) * z[rqj]

計算する要素の順序がかなりややこしいこの演算は「バタフライ」演算などと呼ばれる. このように添え字の 2 進展開を利用して回転子をうまくまとめて処理することで, 演算回数を劇的に減らすことができる. DFT の定義に忠実に計算すると O(N^2) の計算量となるが, 上記の方法を採用すれば O(N*log(N)) の計算量で済む.

実際に処理するときには, いわゆる「インプレイス」型の処理を行うことが多いが, それには変換元データを変換後データで置き換えていかなければならないため,

y[rij] = Σ[h] W(h * r, 2) * x[hij]
z[rqj] = Σ[i] W(i * q, 2) * W(i *  r, 4) * y[rij]
X[rqp] = Σ[j] W(j * p, 2) * W(j * qr, 8) * z[rqj]

などとする必要がある. このとき X の添え字がビット逆転してしまっていることに注意. このため, X が得られた後, いわゆる「ビットリバース」処理を行うことになる.

応用

ここで説明した時系列データの添え字の上位桁から下位桁に向かって DFT を実行していくアルゴリズムは, 「周波数間引き型 FFT」と呼ばれている. 逆に下位桁から上位桁に向かって実行すると「時間間引き型 FFT」となるが, もはや説明するまでもない.

DFT とほぼ同様の処理で逆 DFT が考えられる. DFT が時系列データを周波数成分に変換する演算なら, 逆 DFT はもちろん周波数成分を時系列データに変換する演算となる.

FFT の目的は, 前述の通り一つは周波数成分を計算することであるが, もう一つ重要な目的として畳み込み, 英語では convolution と呼ばれる計算がある.

多項式 (a[0]x^0 + a[1]x^1 + a[2]x^2 + ... + a[N-1]x^(N-1)) と (b[0]x^0 + b[1]x^1 + b[2]x^2 + ... + b[N-1]x^(N-1)) の積を計算したいとき, 各係数同士の積和が必要となるので, O(N^2) の計算量となるが, それぞれの係数の DFT を求め, 同一周波数成分同士を乗算し, 結果を逆 DFT すると積和が計算できてしまっていると言う一見不思議な性質がある. この計算量は O(N*log(N)) で済む. もちろん FFT 自体がそれほど軽い演算ではないため, 項数が多くないときにはあまり有利にはならないが, 巨大な N についてはほとんど必須と言ってよい工夫となる.

畳み込みに FFT を利用する場合, 周波数領域での要素の並び順はまったく重要ではない (どうせ最後に逆 FFT する) ため, 前述のビットリバース処理は省略することができる.

FFT は複素数体上の演算であるが, 原理上, 任意の数体上の演算として拡張できる. 特に有限体上で処理する FFT (のような処理) のことを FMT, 高速剰余変換と呼び, 多桁の乗算にはこちらの方が重要となる.