📂フーリエ解析

高速フーリエ変換アルゴリズム

概要¹

離散フーリエ変換^DFTは、数式的な定義に従って計算すると、$\mathcal{O}(N^{2})$の時間計算量を持ちますが、以下で説明するアルゴリズムに従って計算すると、時間計算量が$\mathcal{O}(N\log_{2}N)$に低減します。この高速フーリエ変換を使用して離散フーリエ変換を高速に実行することができます。これは高速フーリエ変換^{fast Fourier transform, FFT}と呼ばれています。

構築

2つの数字を掛けて、それを別の数字に加える操作を$1$回の演算^operationと呼びましょう。すると、$\sum\limits_{i=0}^{n-1}a_{n}b_{n}$の値を求めるには$n$回の演算が必要です。

$$ \begin{align*} \sum\limits_{n=0}^{0} a_{n}b_{b} &= a_{0}b_{0} = \overbrace{0 {\color{red} +} a_{0} {\color{red} \times} b_{0}}^{\color{red} 1 \text{ operation}} \\ \sum\limits_{n=0}^{1} a_{n}b_{b} &= a_{0}b_{0} + a_{1}b_{1} = \overbrace{\big( \overbrace{0 {\color{red} +} a_{0} {\color{red} \times} b_{0}}^{\color{red} 1 \text{ operation}} \big) {\color{#5882FA} + } a_{1} {\color{#5882FA} \times} b_{1}}^{\color{#5882FA}2 \text{ operations}} \\ \sum\limits_{n=0}^{2} a_{n}b_{b} &= a_{0}b_{0} + a_{1}b_{1} + a_{2}b_{2} = \overbrace{\bigg( \overbrace{\big( \overbrace{0 {\color{red} +} a_{0} {\color{red} \times} b_{0}}^{\color{red} 1 \text{ operation}} \big) {\color{#5882FA} + } a_{1} {\color{#5882FA} \times} b_{1}}^{\color{#5882FA}2 \text{ operations}} \bigg) {\color{#FE9A2E} + } a_{2} {\color{#FE9A2E} \times} b_{2}}^{\color{#FE9A2E}3 \text{ operations}} \\ \end{align*} $$

さて、離散フーリエ変換の定義を思い出してみましょう。

線型変換 $\mathcal{F}_{N} : \mathbb{C}^{N} \to \mathbb{C}^{N}$を離散フーリエ変換と呼びます。

$$ \mathcal{F}_{N}(\mathbf{a}) = \hat{\mathbf{a}} = \begin{bmatrix} \hat{a}_{0} \\ \hat{a}_{1} \\ \dots \\ \hat{a}_{N-1} \end{bmatrix} ,\quad \hat{a}_{m} = \sum_{n=0}^{N-1}e^{-i2\pi mn /N} a_{n}\quad (0\le m < N) \tag{1} $$

このとき、$\mathbf{a} = \begin{bmatrix} a_{0}& a_{1}& \dots& a_{N-1} \end{bmatrix}^{T}$です。

$\hat{a}_{m}$を計算するには$N$回の演算が必要で、$\hat{\mathbf{a}}$を計算するにはこれを$N$回繰り返す必要があるため、離散フーリエ変換を計算するには合計で$N^{2}$回の演算が必要です。つまり、$\mathcal{O}(N^{2})$の時間計算量を持っています。これは、コンピュータ計算の観点からフーリエ変換がかなりのコストを要することを意味します。

アルゴリズム

データの長さ$N$を合成数$N = N_{1}N_{2}$としましょう。そして、インデックス$m, n$を次のように定義します。

$$ m = m^{\prime}N_{1} + m^{\prime \prime},\quad n = n^{\prime}N_{2} + n^{\prime \prime} $$

すると、$0 \le m^{\prime}, n^{\prime \prime} \le N_{2}-1$および$0 \le m^{\prime \prime}, n^{\prime} \le N_{1}-1$です。$(1)$の指数部分を次のように表現できます。

$$ \begin{align*} e^{-i2\pi mn /N} &= e^{-i2\pi (m^{\prime}N_{1} + m^{\prime \prime})(n^{\prime}N_{2} + n^{\prime \prime})/(N_{1}N_{2})} \\ &= e^{-i2\pi (m^{\prime}n^{\prime}N_{1}N_{2} + m^{\prime}n^{\prime \prime}N_{1} + m^{\prime \prime}n^{\prime}N_{2} + m^{\prime \prime}n^{\prime \prime})/(N_{1}N_{2})} \\ &= e^{-i2\pi m^{\prime}n^{\prime}} e^{-i2\pi \big( (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime}/N_{1}) + (m^{\prime \prime}n^{\prime \prime}/N) \big)} \\ &= e^{-i2\pi \big( (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime}/N_{1}) + (m^{\prime \prime}n^{\prime \prime}/N) \big)} \\ \end{align*} $$

これを$(1)$に代入すると、

$$ \begin{align*} \hat{a}_{m} &= \sum_{n=0}^{N-1}e^{-i2\pi mn /N} a_{n} \\ &= \sum_{n^{\prime \prime}=0}^{N_{2}-1} \sum_{n^{\prime}=0}^{N_{1}-1} e^{-i2\pi \big( (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime}/N_{1}) + (m^{\prime \prime}n^{\prime \prime}/N) \big)}a_{n^{\prime}N_{2}+n^{\prime \prime}} \\ &= \sum_{n^{\prime \prime}=0}^{N_{2}-1} \left[ \sum_{n^{\prime}=0}^{N_{1}-1} e^{-i2\pi m^{\prime \prime}n^{\prime}/N_{1}}a_{n^{\prime}N_{2}+n^{\prime \prime}} \right] e^{-i2\pi [ (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime \prime}/N) ] } \end{align*} $$

上記の式に従うと、各括弧内の$\left[ \sum_{n^{\prime}=0}^{N_{1}-1} \right]$を計算するのに$N_{1}$回の演算、括弧の外の$\sum_{n^{\prime \prime}=0}^{N_{2}-1}$を計算するのに$N_{2}$回の演算が必要です。したがって、$\hat{a}_{m}$を計算するには合計で$(N_{1} + N_{2})$回の演算が必要です。$\hat{\mathbf{a}}$を得るにはこれを$N$回繰り返す必要があるため、合計で$N(N_{1} + N_{2})$のコストがかかり、$N^{2}$よりも減少することが確認できます。

括弧内を注意深く見ると、$N_{1}$が再び合成数の場合、同じロジックを適用できることがわかるでしょう。したがって、データの長さが$N = N_{1} N_{2} \cdots N_{k}$といった合成数の場合、時間計算量は次のように減少します。

$$ \mathcal{O}(N^{2}) \searrow \mathcal{O}\big( N(N_{1} + N_{2} + \cdots + N_{k}) \big) $$

ここで、$N$を$2$のべき乗$N = 2^{k}$と仮定してみましょう。すると、$\log_{2}N = k$であり、$N^{2} = 2^{k}$から$2^{k}(2k)$だけ減少するため、時間計算量は次のように減少します。

$$ \mathcal{O}(N^{2}) \searrow \mathcal{O}(2N \log_{2}N) $$

補足

これは1965年にCooleyとTukey²によって提案されたため、Cooley-Tukeyアルゴリズムとも呼ばれています。ただし、彼らが最初に発明したわけではありません。ガウスも同様のアルゴリズムを研究しましたが、正しく発表しなかったため、この事実は後に明らかになりました³。