高速フーリエ変換アルゴリズム
📂フーリエ解析 高速フーリエ変換アルゴリズム 概要 離散フーリエ変換DFT は、数式的な定義に従って計算すると、O ( N 2 ) \mathcal{O}(N^{2}) O ( N 2 ) の時間計算量を持ちますが、以下で説明するアルゴリズムに従って計算すると、時間計算量がO ( N log 2 N ) \mathcal{O}(N\log_{2}N) O ( N log 2 N ) に低減します。この高速フーリエ変換を使用して離散フーリエ変換を高速に実行することができます。これは高速フーリエ変換 fast Fourier transform, FFT と呼ばれています。
構築 2つの数字を掛けて、それを別の数字に加える操作を1 1 1 回の演算operation と呼びましょう。すると、∑ i = 0 n − 1 a n b n \sum\limits_{i=0}^{n-1}a_{n}b_{n} i = 0 ∑ n − 1 a n b n の値を求めるにはn n n 回の演算が必要です。
∑ n = 0 0 a n b b = a 0 b 0 = 0 + a 0 × b 0 ⏞ 1 operation ∑ n = 0 1 a n b b = a 0 b 0 + a 1 b 1 = ( 0 + a 0 × b 0 ⏞ 1 operation ) + a 1 × b 1 ⏞ 2 operations ∑ n = 0 2 a n b b = a 0 b 0 + a 1 b 1 + a 2 b 2 = ( ( 0 + a 0 × b 0 ⏞ 1 operation ) + a 1 × b 1 ⏞ 2 operations ) + a 2 × b 2 ⏞ 3 operations
\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*}
n = 0 ∑ 0 a n b b n = 0 ∑ 1 a n b b n = 0 ∑ 2 a n b b = a 0 b 0 = 0 + a 0 × b 0 1 operation = a 0 b 0 + a 1 b 1 = ( 0 + a 0 × b 0 1 operation ) + a 1 × b 1 2 operations = a 0 b 0 + a 1 b 1 + a 2 b 2 = ( ( 0 + a 0 × b 0 1 operation ) + a 1 × b 1 2 operations ) + a 2 × b 2 3 operations
さて、離散フーリエ変換 の定義を思い出してみましょう。
線型変換 F N : C N → C N \mathcal{F}_{N} : \mathbb{C}^{N} \to \mathbb{C}^{N} F N : C N → C N を離散フーリエ変換 と呼びます。
F N ( a ) = a ^ = [ a ^ 0 a ^ 1 … a ^ N − 1 ] , a ^ m = ∑ n = 0 N − 1 e − i 2 π m n / N a n ( 0 ≤ m < N ) (1)
\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}
F N ( a ) = a ^ = a ^ 0 a ^ 1 … a ^ N − 1 , a ^ m = n = 0 ∑ N − 1 e − i 2 πmn / N a n ( 0 ≤ m < N ) ( 1 )
このとき、a = [ a 0 a 1 … a N − 1 ] T \mathbf{a} = \begin{bmatrix} a_{0}& a_{1}& \dots& a_{N-1} \end{bmatrix}^{T} a = [ a 0 a 1 … a N − 1 ] T です。
a ^ m \hat{a}_{m} a ^ m を計算するにはN N N 回の演算が必要で、a ^ \hat{\mathbf{a}} a ^ を計算するにはこれをN N N 回繰り返す必要があるため、離散フーリエ変換を計算するには合計でN 2 N^{2} N 2 回の演算が必要です。つまり、O ( N 2 ) \mathcal{O}(N^{2}) O ( N 2 ) の時間計算量 を持っています。これは、コンピュータ計算の観点からフーリエ変換がかなりのコストを要することを意味します。
アルゴリズム データの長さN N N を合成数 N = N 1 N 2 N = N_{1}N_{2} N = N 1 N 2 としましょう。そして、インデックスm , n m, n m , n を次のように定義します。
m = m ′ N 1 + m ′ ′ , n = n ′ N 2 + n ′ ′
m = m^{\prime}N_{1} + m^{\prime \prime},\quad n = n^{\prime}N_{2} + n^{\prime \prime}
m = m ′ N 1 + m ′′ , n = n ′ N 2 + n ′′
すると、0 ≤ m ′ , n ′ ′ ≤ N 2 − 1 0 \le m^{\prime}, n^{\prime \prime} \le N_{2}-1 0 ≤ m ′ , n ′′ ≤ N 2 − 1 および0 ≤ m ′ ′ , n ′ ≤ N 1 − 1 0 \le m^{\prime \prime}, n^{\prime} \le N_{1}-1 0 ≤ m ′′ , n ′ ≤ N 1 − 1 です。( 1 ) (1) ( 1 ) の指数部分を次のように表現できます。
e − i 2 π m n / N = e − i 2 π ( m ′ N 1 + m ′ ′ ) ( n ′ N 2 + n ′ ′ ) / ( N 1 N 2 ) = e − i 2 π ( m ′ n ′ N 1 N 2 + m ′ n ′ ′ N 1 + m ′ ′ n ′ N 2 + m ′ ′ n ′ ′ ) / ( N 1 N 2 ) = e − i 2 π m ′ n ′ e − i 2 π ( ( m ′ n ′ ′ / N 2 ) + ( m ′ ′ n ′ / N 1 ) + ( m ′ ′ n ′ ′ / N ) ) = e − i 2 π ( ( m ′ n ′ ′ / N 2 ) + ( m ′ ′ n ′ / N 1 ) + ( m ′ ′ n ′ ′ / N ) )
\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*}
e − i 2 πmn / N = e − i 2 π ( m ′ N 1 + m ′′ ) ( n ′ N 2 + n ′′ ) / ( N 1 N 2 ) = e − i 2 π ( m ′ n ′ N 1 N 2 + m ′ n ′′ N 1 + m ′′ n ′ N 2 + m ′′ n ′′ ) / ( N 1 N 2 ) = e − i 2 π m ′ n ′ e − i 2 π ( ( m ′ n ′′ / N 2 ) + ( m ′′ n ′ / N 1 ) + ( m ′′ n ′′ / N ) ) = e − i 2 π ( ( m ′ n ′′ / N 2 ) + ( m ′′ n ′ / N 1 ) + ( m ′′ n ′′ / N ) )
これを( 1 ) (1) ( 1 ) に代入すると、
a ^ m = ∑ n = 0 N − 1 e − i 2 π m n / N a n = ∑ n ′ ′ = 0 N 2 − 1 ∑ n ′ = 0 N 1 − 1 e − i 2 π ( ( m ′ n ′ ′ / N 2 ) + ( m ′ ′ n ′ / N 1 ) + ( m ′ ′ n ′ ′ / N ) ) a n ′ N 2 + n ′ ′ = ∑ n ′ ′ = 0 N 2 − 1 [ ∑ n ′ = 0 N 1 − 1 e − i 2 π m ′ ′ n ′ / N 1 a n ′ N 2 + n ′ ′ ] e − i 2 π [ ( m ′ n ′ ′ / N 2 ) + ( m ′ ′ n ′ ′ / N ) ]
\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*}
a ^ m = n = 0 ∑ N − 1 e − i 2 πmn / N a n = n ′′ = 0 ∑ N 2 − 1 n ′ = 0 ∑ N 1 − 1 e − i 2 π ( ( m ′ n ′′ / N 2 ) + ( m ′′ n ′ / N 1 ) + ( m ′′ n ′′ / N ) ) a n ′ N 2 + n ′′ = n ′′ = 0 ∑ N 2 − 1 [ n ′ = 0 ∑ N 1 − 1 e − i 2 π m ′′ n ′ / N 1 a n ′ N 2 + n ′′ ] e − i 2 π [( m ′ n ′′ / N 2 ) + ( m ′′ n ′′ / N )]
上記の式に従うと、各括弧内の[ ∑ n ′ = 0 N 1 − 1 ] \left[ \sum_{n^{\prime}=0}^{N_{1}-1} \right] [ ∑ n ′ = 0 N 1 − 1 ] を計算するのにN 1 N_{1} N 1 回の演算、括弧の外の∑ n ′ ′ = 0 N 2 − 1 \sum_{n^{\prime \prime}=0}^{N_{2}-1} ∑ n ′′ = 0 N 2 − 1 を計算するのにN 2 N_{2} N 2 回の演算が必要です。したがって、a ^ m \hat{a}_{m} a ^ m を計算するには合計で( N 1 + N 2 ) (N_{1} + N_{2}) ( N 1 + N 2 ) 回の演算が必要です。a ^ \hat{\mathbf{a}} a ^ を得るにはこれをN N N 回繰り返す必要があるため、合計でN ( N 1 + N 2 ) N(N_{1} + N_{2}) N ( N 1 + N 2 ) のコストがかかり、N 2 N^{2} N 2 よりも減少することが確認できます。
括弧内を注意深く見ると、N 1 N_{1} N 1 が再び合成数の場合、同じロジックを適用できることがわかるでしょう。したがって、データの長さがN = N 1 N 2 ⋯ N k N = N_{1} N_{2} \cdots N_{k} N = N 1 N 2 ⋯ N k といった合成数の場合、時間計算量は次のように減少します。
O ( N 2 ) ↘ O ( N ( N 1 + N 2 + ⋯ + N k ) )
\mathcal{O}(N^{2}) \searrow \mathcal{O}\big( N(N_{1} + N_{2} + \cdots + N_{k}) \big)
O ( N 2 ) ↘ O ( N ( N 1 + N 2 + ⋯ + N k ) )
ここで、N N N を2 2 2 のべき乗N = 2 k N = 2^{k} N = 2 k と仮定してみましょう。すると、log 2 N = k \log_{2}N = k log 2 N = k であり、N 2 = 2 k N^{2} = 2^{k} N 2 = 2 k から2 k ( 2 k ) 2^{k}(2k) 2 k ( 2 k ) だけ減少するため、時間計算量は次のように減少します。
O ( N 2 ) ↘ O ( 2 N log 2 N )
\mathcal{O}(N^{2}) \searrow \mathcal{O}(2N \log_{2}N)
O ( N 2 ) ↘ O ( 2 N log 2 N )
補足 これは1965年にCooleyとTukey によって提案されたため、Cooley-Tukeyアルゴリズム とも呼ばれています。ただし、彼らが最初に発明したわけではありません。ガウスも同様のアルゴリズムを研究しましたが、正しく発表しなかったため、この事実は後に明らかになりました 。