高速フーリエ変換その 1 & 添字の和の畳み込み

まえがき

この記事は投稿者(NokonoKotlin)の個人サイトの記事から Qiita 用に移植 & 加筆したものです。 (この文言は記事が剽窃でないことを周知するためのものであり、個人サイトの宣伝ではないことをあらかじめご了承ください。)

離散フーリエ変換

長さ $N$ の数列 $A$ から以下の数列 $F_{A}$ への変換を離散フーリエ変換と言います。

• $F_{A}[x] = \sum_{0\leq p < N} (A[p]\cdot e^{\frac{-2\pi i}{N}\cdot px})$

ここで ${\omega_N} := e^{\frac{-2\pi i}{N}}$ として、今後は離散フーリエ変換を以下で表記します。

• $F_{A}[x] = \sum_{0\leq p < N} (A[p]\cdot {\omega}^{px}_N)$

この変換をナイーブに実装すると要素 $1$ つあたり $O(N)$ 時間かかってしまうので、高速化を検討してみましょう。

その前に $A$ のサイズを二冪のサイズに整形し、$N = 2^{\log_2{N}+0.999...}$ とします。余分に確保した領域は $0$ で初期化しておけば良いです。

離散フーリエ変換の右辺を行列の積で表します。特に、$A$ にかける係数の行列を $M_N$ とします。例えば $N = 8$ の場合の $M_N \cdot A $ は以下のようになります。

以降、$N = 8$ の例をベースに解説していきます。

高速フーリエ変換 (FFT) の天才的な式変形

$A$ の要素のうち、偶数番目の要素を前半に、奇数
番目の要素を後半に移動し、これを新たに $A$ とします。計算結果 $F_A$ を保持するために係数の行列 $M_8$ の偶数番目の列を前半、奇数番目の列を後半に移動し、これを新たに $M_8$ とします。

計算結果はそのままで、右辺の行列の中身を整理することができました。

ここで ${\omega_k} = {\omega}^2_{2k}$ であることと、${\omega}^{0}_{k} = {\omega}^{k}_{k} = e^{-2\pi i} = e^{2\pi i} = 1$ であることに注目し、上記の行列を ${\omega}_{4} = {\omega}^2_{8}$ で括ることにします。

${\omega}_{4} = {\omega}^2_{8}$ で括ると $M_8$ は以上のようになります。ただしただ括っただけではなく、意図して少し変な括り方をしました。気になるなら、${\omega}^{k}_{k} = 1 $ であることを踏まえて展開すると、元の行列に一致させることが可能であることを確認できます。

特に $\omega_4$ に注目すると、行列全体が以下の行列が登場する $4$ つのブロックに分割されていることを確認できます。

これは $N = 4$ の場合の離散フーリエ変換の係数の行列 $M_4$ です。ここで、$M_8$ から $M_4$ を取り出して分解してみましょう。$M_8$ の左半分は単に $M_4$ が縦に二つ並んでおり、右半分は $M_4$ の各成分に ${\omega}_{8}$ が規則的にくっついています。この規則から、$M_8 \cdot A$ を以下のように分解することができます。

さて、右 $2$ つの行列の積から計算しましょう。これはサイズが $\frac{N}{2} = 4$ の FFT を $2$ 回行うことで計算できます。

列 $A[0],A[2],A[4],A[6]$ の FFT の結果を $f_A[0],f_A[1],f_A[2],f_A[3]$ とし, 列 $A[1],A[3],A[5],A[7]$ の FFT の結果を $f_A[4],f_A[5],f_A[6],f_A[7]$ とします。すると $ M_N\cdot A$ は以下に変形できます。

$N \times N$ 行列の掛け算が発生しますが、中身がスカスカなので愚直に積を取る必要はありません。$F_A[x]$ は以下の式で計算できます。

• $F_A[x] = f_A[x \% \frac{N}{2}] + f_A[( x \% \frac{N}{2} )+ \frac{N}{2}]\cdot {\omega}^{x}_{N}$
• $x\%y$ は $x$ を $y$ で割ったあまり

この計算は $O(N)$ 時間でできるので、あとの計算量解析はサイズ $\frac{N}{2}$ の FFT を計算する部分 ($f_A$ の計算) に依存します。サイズ $N$ の FFT の計算量を $T(N)$ とすると、$T(N) = O(N) + T(\frac{N}{2})$ です。再帰の深さが $O(\log{N})$ なので、$T(N) = O(N)\cdot O(log{N}) = O(N\log{N})$ です。よって FFT は $O(N\log{N})$ 時間で処理できます。

おさらい

• 長さ $N$ の列 $A$ を FFT するために列の長さを二冪に整形し、新たにできた領域を $0$ で埋める。これを新たに $A$ とし、長さを $N$ とする。
• フーリエ変換を行列の積 $M_N\cdot A$ の形で表現する。
• $A$ のインデックスが偶数のものを $A$ の前半に、奇数のものを後半に移動し、行列の積の結果が変わらないように行列の列も移動しておく。
• 行列を $\omega_{\frac{N}{2}}$ で括り、$M_N$ から $M_{\frac{N}{2}}$ を抽出して分解する。
• 右二つの行列の積は、$A$ の前半と後半それぞれにサイズ $\frac{N}{2}$ の FFT を行ったものを結合したものであり、それを $f_A$ とする。
• 残った行列は中身がスカスカなので、$F_A $ の計算は以下の式で $O(N)$ 時間で行える。
  • $F_A[x] = f_A[x \% \frac{N}{2}] + f_A[( x \% \frac{N}{2} )+ \frac{N}{2}]\cdot {\omega}^{x}_{N}$

逆変換

$\omega$ の指数部分の符号を逆にして上記の計算をしてから、全ての要素を $N$ で割るだけです。ただし、再帰的な実装をする場合に要素を割る操作が重複しないように注意してください。再帰的な実装で逆変換を行う場合、指数部分の符号を逆にした計算を先に完了させ、全ての再帰関数呼び出しが完了してから最上位のスコープで全ての要素を $N$ で割ってください。

畳み込み (Convolution)

$2$ つの列 $A,B$ に対して以下の長さ $|A| + |B| - 1$ の列 $C$ を計算することを、添字の和の畳み込みと言います。

• $C_t := \sum_{i+j = t}(A_i \cdot B_j)$
• $A,B$ のインデックス範囲外の部分は $0$ としておく。

結論から言うと、$A,B$ それぞれのサイズを $|A|+|B|-1$ にリサイズ (新たな領域は $0$ とする) して、それを離散フーリエ変換したものをそれぞれ $F_A,F_B$ とします。 $F_C[t] := F_A[t]\cdot F_B[t]$ なる列 $F_C$ を離散フーリエ逆変換すると、$A,B$ の畳み込み $C$ が得られます(不思議ですね〜)。