高速フーリエ変換を自前で作りたい (3) - 離散フーリエ変換・FFTの実装方法

前回に引き続き、今回は離散フーリエ変換と、具体的な高速フーリエ変換の実装方法についてまとめてみます。離散フーリエ変換は、一定のサンプリングレートで取られた離散データに対してフーリエ変換を行う方法です。高速フーリエ変換はこれを高速で処理できるように実装したものです。

離散フーリエ変換

前回導出したフーリエ変換の定義式は次のものでした。

F(\omega) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(x)e^{-i\omega x}dx

これは、ある周波数に対しての波形を為す関数を得る式でした。離散フーリエ変換の場合、関数が連続ではないため、データ数が限られてしまいます。そこで、これを次のように和の形で表してあげます。データの数列 $f_{n}$ 、データ数 $N$ として、

F(t) = \sum_{n=0}^{N-1} f_{n}e^{\frac{-2n\pi it}{N}} (0 \leq t \leq N)

$\frac{2n\pi}{N}$ と $x$ が置き換わっています。フーリエ変換は非周期関数を無限大の周期を持つ関数としてフーリエ級数展開を行う処理でしたね。これを限られたデータ数Nにおいて実現するためにこのような処理をしています。

これをそのまま実装した場合、すべての $t \in N$ について、 $N$ 回の足し算を行うことになるので、計算量は $O(N^{2})$ となります。

私の場合Web AudioのAudio Workletを用いた処理にFFTを使いたいと考えているのですが、Audio Workletでは一度に128個の信号を取得できます。この時、計算量は $128^{2} = 16384$ となります。

この計算量を抑える工夫をしてみるのが、高速フーリエ変換です。今回はWikipediaにデカデカと載っているクーリー・テューキー型FFTアルゴリズムを取り上げます。

複素数を用いることで、1のn乗根を求めることができます。

\epsilon_{n}^{k} = e^{\frac{-2ik\pi}{N}} = cos(\frac{2k\pi}{n}) + isin(\frac{2k\pi}{n})

この一つ $W_{N} = e^{\frac{-2i\pi}{N}}$ 用いると、離散フーリエ変換の式は次のように表せます。

F(t) = \sum_{n=0}^{N-1} f_{n}W_{N}^{tn}

ここで、N=4として、フーリエ級数を行列を用いて表してみます。

\begin{pmatrix} X_{0}\\ X_{1}\\ X_{2}\\ X_{3}\\ \end{pmatrix} = \begin{pmatrix} W^{0} & W^{0} & W^{0} & W^{0}\\ W^{0} & W^{1} & W^{2} & W^{3}\\ W^{0} & W^{2} & W^{4} & W^{6}\\ W^{0} & W^{3} & W^{6} & W^{9}\\ \end{pmatrix} \begin{pmatrix} f_{0}\\ f_{1}\\ f_{2}\\ f_{3}\\ \end{pmatrix}

行を添え字の偶奇によって分け、

= \begin{pmatrix} W^{0} & W^{0} & W^{0} & W^{0}\\ W^{0} & W^{2} & W^{1} & W^{3}\\ W^{0} & W^{4} & W^{2} & W^{6}\\ W^{0} & W^{6} & W^{3} & W^{9}\\ \end{pmatrix} \begin{pmatrix} f_{0}\\ f_{2}\\ f_{1}\\ f_{3}\\ \end{pmatrix}

こうすると、N=2の結果を用いて演算できます。

\begin{pmatrix} X_{0}\\ X_{1}\\ X_{2}\\ X_{3}\\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & W^{0} & 0\\ 0 & 1 & 0 & W^{1}\\ 1 & 0 & W^{2} & 0\\ 0 & 1 & 0 & W^{3}\\ \end{pmatrix} \begin{pmatrix} W_{2}^{0} & W_{2}^{0} & 0 & 0\\ W_{2}^{0} & W_{2}^{1} & 0 & 0\\ 0 & 0 & W_{2}^{0} & W_{2}^{0}\\ 0 & 0 & W_{2}^{0} & W_{2}^{0}\\ \end{pmatrix} \begin{pmatrix} f_{0}\\ f_{2}\\ f_{1}\\ f_{3}\\ \end{pmatrix}

これを繰り返し行います。計算量は $N$ 個のデータに対して $log(N)$ 回繰り返すので、 $O(Nlog(N))$ となります。

このアルゴリズムは他のアルゴリズムと組み合わせることもできます。

フーリエ変換を離散データに対して行う離散フーリエ変換と、その具体的な実装方法である高速フーリエ変換についてまとめました。

3回に渡って高速フーリエ変換を理解するための理屈や知識をまとめてみました。書いていてとても勉強になったので、こういう記事はこれからも書いていきたいです。