以前フーリエ変換をやりましたね。今回は説明し損ねていたフーリエ変換の高速化の話です。離散フーリエ変換(DFT)というものはご存知でしょうか。以下のような配列 $f(j)$ を配列 $F(k)$ へ変換するアルゴリズムのことです。

\begin{align} F(k)=\sum_{j=0}^{N-1}f(j)e^{-2\pi ijk/N},\ k=0,1,\dots,N-1 \end{align}

フーリエ級数の係数を取り出すための積分演算が和で置き換わっていますね。DFT はデジタル信号処理にはなくてはならない非常に重要な技術です。しかし計算量に着目すると、1 つの $k$ に対して $F(k)$ を計算するためにはループ文を走らせなければなりません。したがってすべての $k$ に対して $F(k)$ を計算して変換を完了させるためには、2 重のループ文を実行し、ループ文の中身を $N^2$ 回計算する必要があります。

このループ回数を劇的に減らしたアルゴリズムが、これから説明する高速フーリエ変換(FFT)です。FFT の登場により信号処理分野は大きく発展しました。

FFT は $N$ の素因数分解を使います。説明を簡単にするために 2 つの素数 $p,q$ により $N=pq$ となっている前提で話を進めましょう。

配列 $f,F$ のレコード数を、$p\times q$ の 2 次元配列と考えてみましょう。$p$ と $q$ のどちらを行と考えるかにより 2 通りの表現方法があります。

\begin{align} j=qx_1+x_2,\ k=y_1+py_2 \end{align}

このような 2 次元配列として考えると $jk/N$ は以下のようになります。

\begin{align} \dfrac{jk}{N}&=\dfrac{pqx_1y_2+qx_1y_1+px_2y_2+x_2y_1}{pq}\\ &=x_1y_2+\dfrac{x_1y_1}{p}+\dfrac{x_2y_2}{q}+\dfrac{x_2y_1}{pq} \end{align}

これを指数関数の引数に代入すると

\begin{align} e^{-2\pi ijk/N}&=e^{-2\pi ix_1y_2}e^{-2\pi i x_1y_1/p}e^{-2\pi ix_2(y_2/q+y_1/pq)}\\ &=e^{-2\pi ix_1y_1/p}e^{-2\pi ix_2(y_2/q+y_1/pq)} \end{align}

$x_1y_2$ は整数であることと変数 $x_1,x_2$ の分離という点に注意して式を追ってみてください。変数分離によって DFT の総和計算も分離されます。

\begin{align} F(y_1+py_2)&=\sum_{x_1=0}^{p-1}\sum_{x_2=0}^{q-1}f(qx_1+x_2)e^{-2\pi ix_1y_1/p}e^{-2\pi ix_2(y_2/q+y_1/pq)}\\ &=\sum_{x_2=0}^{q-1}\Big(\sum_{x_1=0}^{p-1}f(qx_1+x_2)e^{-2\pi ix_1y_1/p}\Big)e^{-2\pi ix_2(y_2/q+y_1/pq)} \end{align}

変数分離によって計算量が減っていることに気が付くでしょうか。実際に数えてみましょう。内側のループ文はすべての $x_2,y_1$ に対して実行しなければならないため、計算回数は $pN$ です。外側のループ文に関しても同様にすべての $y_1,y_2$ に対して実行しなけれなならず、計算回数は $qN$ です。したがって合計で $(p+q)N$ 回の計算が必要になります。

DFT と FFT を比較すると $N$ の部分が $p+q$ になって減少していることが分ります。このことから更に $N$ の素因数分解が細かいほど高速化の効果は高い、すなわち $N$ が 2 の冪乗となる場合に最も高い高速化効果が得られることがわかります。

ご興味あれば $N$ が 2 の冪乗となる場合の FFT を勉強してみると良いでしょう。次回は FFT を用いた畳込み演算の高速化の話をします。(written by J.Nishigata)

その他9- 高速フーリエ変換、高速化の秘密を知ろう