Chapter 03

高速フーリエ変換

a3geek
a3geek
2021.05.02に更新

高速フーリエ変換の導出

 計算機で離散フーリエ変換の計算をするとき、一つの点kについての係数を計算するためにはN回の加算が行われる。よって、全ての点kについて係数を計算するときの計算量がO(N^2)となり、計算量が多くなる。また、積の計算は一般に計算機への負荷が高い処理であり、可能な限りこれを省略したい。そこで、フーリエ変換の計算をするときの計算量や計算機への負荷を削減するアルゴリズムを考える。
 計算量や計算機への負荷を削減するために、複素指数関数の特性である周期性と対称性を利用する。周期性とは、複素平面において角度が一周すると一周前と同じ値を取る特性である。対称性とは、複素平面上の単位円の中心に対して、点対称の位置では元の点と符号だけが異なる同じ値を取る特性である。
 離散フーリエ変換の複素指数関数部の表記を簡単化するために

\tag{34} W_N = e^{-j\frac{2\pi}{N}}

とする。図1N = 8のときのW_Nの複素平面上での値を示す。


図1 N = 8の複素平面

 図1より、W_Nが複素平面の単位円上に均等に配置され、周期性と対称性を有していることが分かる。式(34)を用いて離散フーリエ変換を書き直すと

\tag{35} F[k] = \sum_{n = 0}^{N - 1}f[n]W_N^{nk}

が得られる。ここで、f[n]が1周期内で取る値の数N2の累乗であると仮定する。式(34)の対称性を利用するために点kを偶数要素と奇数要素に分割すると、

\tag{36} \begin{aligned} F[2k] &= \sum_{n=0}^{N-1}f[n]W_N^{2nk} \\ &= \sum_{n=0}^{\frac{N}{2}-1}f[n]W_N^{2nk} + \sum_{n=0}^{\frac{N}{2}-1}f[n + \frac{N}{2}]W_N^{2(n + \frac{N}{2})k} \\ &= \sum_{n=0}^{\frac{N}{2}-1}f[n]W_{\frac{N}{2}}^{nk} + \sum_{n=0}^{\frac{N}{2}-1}f[n + \frac{N}{2}]W_{\frac{N}{2}}^{(n + \frac{N}{2})k} \quad (\because W_N^{2k} = W_{\frac{N}{2}}^k) \\ &= \sum_{n=0}^{\frac{N}{2}-1}(f[n] + f[n + \frac{N}{2}])W_{\frac{N}{2}}^{nk} \quad (\because W_{\frac{N}{2}}^{n+\frac{N}{2}} = W_{\frac{N}{2}}^{n}) \\ F[2k + 1] &= \sum_{n=0}^{N-1}f[n]W_N^{n(2k+1)} \\ &= \sum_{n=0}^{\frac{N}{2}-1}f[n]W_N^{n(2k+1)} + \sum_{n=0}^{\frac{N}{2}-1}f[n + \frac{N}{2}]W_N^{(n + \frac{N}{2})(2k+1)} \\ &= \sum_{n=0}^{\frac{N}{2}-1}(f[n] - f[n + \frac{N}{2}])W_{\frac{N}{2}}^{nk}W_N^n \end{aligned}

が得られる。式(36)より、点kにおける出力であるF[k]について、偶数要素は要素数\frac{N}{2}f[n] + f[n + \frac{N}{2}]による離散フーリエ変換に等しく、奇数要素は要素数\frac{N}{2}(f[n] - f[n + \frac{N}{2}])W_N^nによる離散フーリエ変換に等しいことが分かる。つまり、要素数Nの離散フーリエ変換は、要素数\frac{N}{2}の離散フーリエ変換を2回計算し加算することによって求められる。
 このとき、2つの離散フーリエ変換の計算量はそれぞれO(\frac{N}{2})であるため、元の離散フーリエ変換の計算量はO(\frac{N^2}{4} + \frac{N^2}{4}) = O(\frac{N^2}{2})となる。これは、素朴な離散フーリエ変換の計算量の半分である。さらに、分割した離散フーリエ変換は、再帰的に要素数が半分の離散フーリエ変換に分割していくことができる。Nは2の累乗であると仮定しているため、分割を\log_2 N回繰り返すことで要素数1の離散フーリエ変換になる。分割を繰り返す際に\frac{N}{2}回のW_Nの乗算を行うため、最終的な計算量はO(N \log_2N)となる。
 このような、要素を偶奇に分割することで計算量を削減し、高速に離散フーリエ変換を計算するアルゴリズムを(基数2の)Cooley-Tukeyのアルゴリズム(クーリー - テューキーのアルゴリズム)といい、高速化した離散フーリエ変換を高速フーリエ変換(Fast Fourier Transform, 以下「FFT」とする)という。また、式(34)W_Nを回転因子、または回転子という。なお、式(36)では周波数領域を偶奇で分割することを考えた。周波数領域に対して分割を行うことを周波数間引きといい、これに対して時間領域に対して分割を行うことを時間間引きという。
 時間間引きに関するFFTについても同様に考える。離散フーリエ変換における総和内の各要素を添字nの偶奇で分けると、

\tag{37} \begin{aligned} F[k] &= \sum_{n = 0}^{N - 1}f[n]e^{-j\frac{2\pi}{N}nk} \\ &= \sum_{n = 0}^{\frac{N}{2} - 1}f[2n]e^{-j\frac{2\pi}{N}2nk} + \sum_{n = 0}^{\frac{N}{2}-1}f[2n + 1]e^{-j\frac{2\pi}{N}(2n+1)k} \\ &= \sum_{n = 0}^{\frac{N}{2}-1}f[2n]e^{-j\frac{2\pi}{\frac{N}{2}}nk} + e^{-j\frac{2\pi}{N}k}\sum_{n = 0}^{\frac{N}{2}-1}f[2n+1]e^{-j\frac{2\pi}{\frac{N}{2}}nk} \end{aligned}

が得られる。式(37)に式(34)を代入すると、

\tag{38} F[k] = \sum_{n=0}^{\frac{N}{2}-1}f[2n]W_{\frac{N}{2}}^{nk} + W_N^k\sum_{n=0}^{\frac{N}{2}-1}f[2n+1]W_{\frac{N}{2}}^{nk}

が得られる。式(38)の右辺について、第一項は要素数\frac{N}{2}の離散フーリエ変換であり、第二項は要素数\frac{N}{2}の離散フーリエ変換にW_N^kをかけたものである。すなわち、要素数Nの離散フーリエ変換は、要素数\frac{N}{2}の2つの離散フーリエ変換に対して片方にW_N^kをかけて和を取った処理に分割できる。周波数間引きと同様に、2つの離散フーリエ変換の計算量はそれぞれO(\frac{N}{2})であるので、元の離散フーリエ変換の計算量はO(\frac{N^2}{4} + \frac{N^2}{4}) = O(\frac{N^2}{2})となる。これは、素朴な離散フーリエ変換の半分である。さらに、分割を再帰的に繰り返すことで要素数1の離散フーリエ変換になり、最終的な計算量はO(N \log_2N)となる。
 実際のFFTの計算を考える。要素数NN = 2のとき、式(35)より

\tag{39} \begin{aligned} F[0] &= f[0]W_2^0 + f[1]W_2^0 \\ F[1] &= f[0]W_2^0 + f[1]W_2^{(0 + \frac{2}{2})(0+1)} = f[0]W_2^0 + f[1]W_2^1 \end{aligned}

である。記述を簡単化するために行列で表現することを考える。F[k]=X_{2,k}, f[n]=x_nとすると、式(39)より

\tag{40} \begin{pmatrix}X_{2,0}\\X_{2,1}\end{pmatrix}= \begin{pmatrix} W_2^0&W_2^0\\ W_2^0&W_2^1\end{pmatrix} \begin{pmatrix}x_0\\x_1\end{pmatrix} = \begin{pmatrix} 1&1\\ 1&-1\end{pmatrix} \begin{pmatrix}x_0\\x_1\end{pmatrix}

が得られる。また、N = 4のときF[k]=X_{4,k},f[n]=x_nとし、W = W_4と略記すると、

\tag{41} \begin{pmatrix}X_{4,0}\\X_{4,1}\\X_{4,2}\\X_{4,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}x_0\\x_1\\x_2\\x_3\end{pmatrix}

である。このとき、周波数間引きのためにkについて偶奇で分けて行列を変形すると、式(41)より

\begin{aligned} {\begin{pmatrix}X_{4,0}\\X_{4,2}\\X_{4,1}\\X_{4,3}\end{pmatrix}} &={\begin{pmatrix} W^0&W^0&W^0&W^0\\ W^0&W^2&W^4&W^6\\ W^0&W^1&W^2&W^3\\ W^0&W^3&W^6&W^9 \end{pmatrix}}\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^0 & W^0 & -W^0 \\ W^0 & W^1 & -W^0 & -W^1 \\ W^0 & -W^1 & -W^0 & W^1\\ \end{pmatrix}}{\begin{pmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\end{pmatrix}} \\ &(\because W^{k+N} = W^{k}, W^{k+\frac{N}{2}} = -W^{k}) \\ &= {\begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 \\ 1 & W^1 & -1 & -W^1 \\ 1 & -W^1 & -1 & W^1\\ \end{pmatrix}}{\begin{pmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\end{pmatrix}} \end{aligned}

が得られる。ここで、N = 2における回転因子の行列を

\tag{43} W'_2= \begin{pmatrix} W_2^0&W_2^0\\ W_2^0&W_2^1\end{pmatrix} = \begin{pmatrix} 1&1\\ 1&-1\end{pmatrix}

とし、行列D_4

\tag{44} D_4 = \begin{pmatrix} 1&W^1\\ 1&-W^1\end{pmatrix}

と定義する。このとき、式(42)に式(43), (44)を代入すると、

\tag{45} \begin{pmatrix}X_{4,0}\\X_{4,2}\\X_{4,1}\\X_{4,3}\end{pmatrix}=\begin{pmatrix}W_2' & W_2' \\ D_4 & -D_4\end{pmatrix}\begin{pmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\end{pmatrix}

が得られる。式(45)より、N = 4の計算においてN = 2で計算した値を利用することができるため、計算量や計算負荷が削減できることが分かる。
 図2N = 2の、図3N = 4の計算を図式化したものを、図4N = 2N = 4に関する計算手順の内包関係をそれぞれ示す。


図2 N = 2の計算手順


図3 N = 4の計算手順


図4 N = 2N = 4に関する計算手順の内包関係

 図2、 図3より、FFTの計算を図式化すると蝶のような形の図になることから、この図をバタフライ図といい、この計算手順をバタフライ演算という。また、図4よりN = 4のバタフライ演算はN = 2のバタフライ演算を組み合わせることで実現できることが分かる。

Split-radix FFTの導出

 FFTのアルゴリズムの一つであるCooley-Tukeyのアルゴリズムによって、f[n]が1周期内で取る値の数N2の累乗である時、要素を偶奇に分けることで高速にフーリエ変換を計算することができる。言い換えると、要素を偶奇に分ける操作は、基数を2として要素を分割する操作である。そこで、要素数が4の累乗であると仮定し、基数を4にした場合を考える。
 離散フーリエ変換において、時間間引きで基数4の分割を行うと、

\tag{46} \begin{aligned} F[k] &= \sum_{n=0}^{\frac{N}{4}-1}f[n]W_N^{nk} + \sum_{n = 0}^{\frac{N}{4}-1}f[n + \frac{N}{4}]W_N^{(n + \frac{N}{4})k} \\ &+ \sum_{n = 0}^{\frac{N}{4}-1}f[n + \frac{N}{2}]W_N^{(n + \frac{N}{2})k} + \sum_{n = 0}^{\frac{N}{4}-1}f[n + \frac{3N}{4}]W_N^{(n + \frac{3N}{4})k} \end{aligned}

が得られる。ここで、回転因子の特性について

\tag{47} \begin{aligned} W_N^{\frac{N}{4}k} &= (\cos\frac{\pi}{2} - j\sin\frac{\pi}{2})^k = (-j)^k \\ W_N^{\frac{N}{2}k} &= (-1)^k \\ W _N^{\frac{3N}{4}k} &= j^k \end{aligned}

であるため、式(46)より

\tag{48} \begin{aligned} F[k] &= \sum_{n = 0}^{\frac{N}{4}-1}f[n]W_N^{nk} + W_N^{\frac{N}{4}k}\sum_{n = 0}^{\frac{N}{4}-1}f[n + \frac{N}{4}]W_N^{nk} \\ &+ W_N^{\frac{N}{2}k}\sum_{n = 0}^{\frac{N}{4}-1}f[n + \frac{N}{2}]W_N^{nk} + W_N^{\frac{3N}{4}k}\sum_{n = 0}^{\frac{N}{4}-1}f[n + \frac{3N}{4}]W_N^{nk} \\ &= \sum_{n = 0}^{\frac{N}{4}-1}f[n]W_N^{nk} + (-j)^k\sum_{n = 0}^{\frac{N}{4}-1}f[n + \frac{N}{4}]W_N^{nk} \\ &+ (-1)^k\sum_{n = 0}^{\frac{N}{4}-1}f[n + \frac{N}{2}]W_N^{nk} + j^k\sum_{n = 0}^{\frac{N}{4}-1}f[n + \frac{3N}{4}]W_N^{nk} \\ &= \sum_{n=0}^{\frac{N}{4}-1}(f[n] + (-j)^kf[n + \frac{N}{4}] + (-1)^kf[n + \frac{N}{2}] + j^kf[n + \frac{3N}{4}])W_N^{nk} \end{aligned}

が得られる。また、周波数間引きで基数4の分割を行うと、式(48)より

\tag{49} \begin{aligned} F[4k] &= \sum_{n=0}^{\frac{N}{4}-1}(f[n] + f[n + \frac{N}{4}] + f[n + \frac{N}{2}] + f[n + \frac{3N}{4}])W_{\frac{N}{4}}^{nk} \\ &(\because W_N^{4nk} = W_{\frac{N}{4}}^{nk}, k \in \mathbb{Z} \Rightarrow (-i)^{4k} = 1) \\ F[4k + 1] &= \sum_{n=0}^{\frac{N}{4}-1}(f[n] - jf[n + \frac{N}{4}] - f[n + \frac{N}{2}] + jf[n + \frac{3N}{4}])W_{\frac{N}{4}}^{nk}W_N^n \\ F[4k + 2] &= \sum_{n=0}^{\frac{N}{4}-1}(f[n] - f[n + \frac{N}{4}] - f[n + \frac{N}{2}] - f[n + \frac{3N}{4}])W_{\frac{N}{4}}^{nk}W_N^{2n} \\ F[4k + 3] &= \sum_{n=0}^{\frac{N}{4}-1}(f[n] + jf[n + \frac{N}{4}] - f[n + \frac{N}{2}] - jf[n + \frac{3N}{4}])W_{\frac{N}{4}}^{nk}W_N^{3n} \end{aligned}

が得られる。
 式(48)(49)より、基数2の分割と同様に、基数4の分割によって要素数Nの離散フーリエ変換は要素数\frac{N}{4}の4つの離散フーリエ変換に分割できることが分かる。このとき、4つの離散フーリエ変換の計算量はそれぞれO(\frac{N}{4})であるので、元の離散フーリエ変換の計算量はO(\frac{N^2}{16} + \frac{N^2}{16} + \frac{N^2}{16} + \frac{N^2}{16}) = O(\frac{N^2}{4})と基数2の分割と同様、\frac{1}{4}となる。分割を再帰的に\frac{1}{2}\log_2N回繰り返すことで、要素数1の離散フーリエ変換になる。
 図5、図6に、F[k]=X_{N,k},f[n]=x_nとしたときの、要素数Nにおける基数2と基数4についてx_k, x_{k + \frac{N}{4}}, x_{k + \frac{N}{2}}, x_{k + \frac{3N}{4}}のみを抽出したバタフライ図を示す。


図5 基数2のバタフライ図の一部


図6 基数4のバタフライ図の一部

 図5, 図6より、基数2と比べて基数4の方が計算量が少なくなることが分かる。まず、図5より基数2における第1段目のバタフライ演算に後置するW_N^k, W_N^{k + \frac{N}{4}}は、共通因子W_N^kを持っている。次に、図6より基数4では共通因子W_N^kを第2段目のバタフライ演算に後置することで、第2段目のW_Nの乗算が増えることになる代わりに第1段目のバタフライ演算に後置する乗算を削減できることが分かる。その結果、基数2ではW_Nの乗算が4箇所なのに対し、基数4では3箇所に減らすことができる。バタフライ演算全体にも同様にすると、乗算の計算量が\frac{3}{4}倍になる。このとき、-jによる乗算は、単に複素数の実部と虚部を入れ替える操作だけで実現できるため、計算量には影響しない。
 基数4のバタフライ演算でも、x_k, x_{k + \frac{N}{4}}から伸びる処理では乗算の計算量に変化はない。ここで、x_k, x_{k+\frac{N}{4}}から伸びる処理はX_{N, k}, X_{N, k + \frac{N}{2}}に対応している。これは、出力の添字が偶数の項に対しては乗算の削減が行われないことを示している。x_{k + \frac{N}{2}}, x_{k + \frac{3N}{4}}から伸びる処理はX_{N, k + \frac{N}{4}}, X_{N, k + \frac{3N}{4}}に対応している。これは、出力の添字が奇数の項に対しては乗算の削減が行われることを示している。
 そこで、偶数の添字の項に対しては基数2の分解を行い、奇数の添字の項に対しては基数4の分解を行うことを考える。式(36), (49)より偶数項と基数項で分割を行うと、

\tag{50} \begin{aligned} F[2k] &= \sum_{n=0}^{\frac{N}{2}-1}(f[n] + f[n + \frac{N}{2}])W_N^{2nk} \\ F[4k + 1] &= \sum_{n=0}^{\frac{N}{4}-1}((f[n] - f[n + \frac{N}{2}]) - j(f[n + \frac{N}{4}] - f[n + \frac{3N}{4}]))W_{N}^{4nk}W_N^n \\ F[4k + 3] &= \sum_{n=0}^{\frac{N}{4}-1}((f[n] - f[n + \frac{N}{2}]) + j(f[n + \frac{N}{4}] - f[n + \frac{3N}{4}]))W_{N}^{4nk}W_N^{3n} \end{aligned}

が得られる。式(50)のアルゴリズムをSplit-Radix(基数分割) FFTという。このアルゴリズムでは基数2のバタフライ演算と基数4のバタフライ演算が混在するため、構造が複雑になる欠点がある。しかし、他の既存のFFTと比較して、乗算と加算の回数が最小になる利点がある。