%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Copyright (C) 1998 Keiichi Ishioka % % % % この著作権表示が記載してあるならば, どのような媒体にでも何も手を加え % % ずに複写して再配付することを許可します. これらのドキュメントとその使 % % 用許諾を他の言語に翻訳する際は, 以下の条件に従って下さい. % % % % 他の言語に翻訳する際に, 他の言語での口語表現への変更の範囲を越えて故 % % 意に意味を変更しないこと. % % % % 翻訳した使用許諾には, それが翻訳であること, また原文の使用許諾が翻訳 % % 全体に依然適用されることを明示すること. % % % % ハイパーテキストの場合, 同一サイト上で原文のコピーを保守すること. ま % % た, 翻訳したハイパーテキストのページからその原文へのリンクを提供する % % こと. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %履歴 98/10/27 石岡圭一 % \documentstyle[a4j]{jarticle} \begin{document} \begin{center} \Large\bf FFT -- 高速アルゴリズムの発見 -- \end{center} \begin{flushright} \bf 石岡圭一 \end{flushright} \section{はじめに} FFT(Fast Fourier Transform = 高速フーリエ変換)は, DFT(Discrete Fourier Transform = 離散フーリエ変換)を高速に実行する手法 です. ここでは, フーリエ変換の意味や性質については既知として, DFTを簡 単に導入したうえで, FFTそのものの構成や歴史に絞って解説していきます. \section{離散フーリエ変換} 周期$2\pi$もつ関数$f(x)$を以下のようにフーリエ級数展開することを考えま す: \begin{equation} f(x)=\sum_{k=-\infty}^{\infty}a_k e^{ikx}. \label{Fourier-Expansion} \end{equation} $f(x)$の値が一周期にわたって連続的に分かっていれば, 展開係数$a_k$の値 は以下のフーリエ積分で求めることができます: \begin{equation} a_k=\frac1{2\pi}\int^{2\pi}_0 f(x)e^{-ikx}dx. \label{Fourier-Integral} \end{equation} しかし, 関数$f(x)$が観測データ等である場合, 関数値は連続的には得られず, 離散的な$x$の値に対してしか関数値が与えられないことが普通です. 例えば, $N$個の等間隔の点$x_j$: \begin{equation} x_j=\frac{2\pi}N j, \quad (j=0,1,\cdots,N-1), \end{equation} での関数値$f_j$: \begin{equation} f_j=f(x_j), \quad (j=0,1,\cdots,N-1), \end{equation} のみが得られている場合を考えましょう. このとき, フーリエ展開 (\ref{Fourier-Expansion})における展開関数$e^{ikx}$は$N$個の点$x_j$でし か評価されませんので, $N$次元のベクトルとみなされます. $N$次元のベクトルですので, 当然一次独立な成分は$N$個しかありえません. 特に今の場合のように点$x_j$を等間隔にとっている場合には, \begin{equation} e^{i(k+N)x_j}=e^{i(k+N)\frac{2\pi}N j} =e^{2\pi i j}e^{ik\frac{2\pi}N j}=e^{ik\frac{2\pi}N j}=e^{ikx_j} \end{equation} となり, 波数$k$を$N$だけずらしたベクトルが波数$k$のベクトルと全く同じ ものになります. これは, 有限の点で関数を評価するために, ある波数の成分が別の波数の成分 と区別がつかなくなって同一視されてしまうことによります. この現象は, ある波数の成分に「別名」を与えてしまうということで, 「エリ アジング」と呼ばれています. 従って, 今の場合, フーリエ級数展開 (\ref{Fourier-Expansion})は無限個の和を考える代わりに, 独立な$N$個の 成分の和: \begin{equation} f_j=\sum_{k=0}^{N-1}a_ke^{ikx_j}, \quad (j=0,1,\cdots,N-1), \label{IDFT} \end{equation} のみを考えるべきであることが分かります. では, 展開(\ref{IDFT})において, 与えられた関数値$f_j$から展開係数$a_k$ を求めるにはどうすればよいでしょうか? すぐ思いつく方法は, 「(\ref{IDFT})の左辺が既知なのだから, (\ref{IDFT})を$N$元連立一次方程 式として解けば$a_k$が求まる」というものですが, 実はもっとうまい方法が あります. 等比級数の和の公式により, 整数$l$が$N$の倍数で ない場合には, \begin{equation} \sum_{j=0}^{N-1}e^{ilx_j} =\sum_{j=0}^{N-1}e^{il\frac{2\pi}{N}j} =\sum_{j=0}^{N-1}\left(e^{il\frac{2\pi}{N}}\right)^j =\frac{\left(e^{il\frac{2\pi}{N}}\right)^N-1}{e^{il\frac{2\pi}{N}}-1}=0, \end{equation} となることが分かりますし, 整数$l$が$N$の倍数の場合には \begin{equation} \sum_{j=0}^{N-1}e^{ilx_j} =\sum_{j=0}^{N-1}e^{il\frac{2\pi}{N}j} =\sum_{j=0}^{N-1}1^j=N, \end{equation} となることが分かります. 従って, (\ref{IDFT})の両辺に$e^{-ikx_j}\quad (k=0,1,\cdots,N-1)$を掛けて$j$について総和をとると, \begin{equation} (左辺)=\sum_{j=0}^{N-1}f_je^{-ikx_j}, \end{equation} \begin{equation} (右辺)=\sum_{j=0}^{N-1}\left(\sum_{k'=0}^{N-1}a_{k'}e^{ik'x_j}\right)e^{-ikx_j} =\sum_{k'=0}^{N-1}a_{k'}\left(\sum_{j=0}^{N-1}e^{i(k'-k)x_j}\right) =Na_k, \end{equation} が得られます. ここで, (\ref{IDFT})における$k$を$k'$に書換え, $k'-k$が$N$の倍数になるのは, $k'=k$となるときだけであることを用いまし た. 以上から, $a_k$を求めるための公式: \begin{equation} a_k=\frac1N\sum_{j=0}^{N-1}f_je^{-ikx_j}, \quad (k=0,1,\cdots,N-1), \label{DFT} \end{equation} が得られたことになります. $f(x)$が連続量として与えられていた場合のフー リエ積分の式(\ref{Fourier-Integral})と見くらべると, この式は, 積分を区 分求積法的な和に置き換えたものになっていることが分かります. 以上により, 関数の評価点が離散的に与えられた場合の展開係数$a_k$と関数 値$f_j$との間の関係式(\ref{IDFT}),(\ref{DFT})が導かれましたが, これ らの式は離散フーリエ変換(Discrete Fourier Transform = DFT)と呼ばれてい ます. \section{高速フーリエ変換} DFTの2つの式(\ref{DFT})と(\ref{IDFT})とは, 互いに変換と逆変換との関係 にありますが, 殆ど同じ形をしていることに気がつきます. 実際, (\ref{DFT})の両辺を$N$倍して複素共役をとると, \begin{equation} (Na_k)^*=\sum_{j=0}^{N-1}f_j^*e^{ikx_j} =\sum_{j=0}^{N-1}f_j^*e^{ik\frac{2\pi}N j} =\sum_{j=0}^{N-1}f_j^*e^{ij\frac{2\pi}N k} =\sum_{j=0}^{N-1}f_j^*e^{ijx_k}, \end{equation} となり, (\ref{IDFT})と全く同じ形に帰着されます. 従って, DFTの計算につ いて考える場合, (\ref{IDFT})式の形のみを考えるだけで十分であることが分 かります. というわけで, 以下(\ref{IDFT})式の計算について考えていくわけ ですが, 便宜上, $x_j=2\pi j/N$であることを使って, 以下のような形に書換 えておきます: \begin{equation} F(j)=\sum_{k=0}^{N-1}W_N^{jk}A(k), \quad (j=0,1,\cdots,N-1). \label{WDFT} \end{equation} ここに, $F(j)=f_j$, $A(k)=a_k$と書き直し, $W_N=e^{\frac{2\pi i}N}$と定 めています. さて, (\ref{WDFT})式の計算について考えていきましょう. この計算は, $N\times N$行列$(W_N^{jk})$に$N$次元ベクトル$a_k$を掛けることに他なり ませんから, 必要な計算量は, 積が$N^2$回と和が$N(N-1)$回で, ともに ${\cal O}(N^2)$回ということになります. 従って, 例えば$N=10000$の場合を考えると, DFTの計算のために $N^2=1$億回オーダーの演算が必要になります. 今では, パソコンの性能が向 上して, 最速のものになると1GFlops(1秒間に10億回の浮動小数点演算を実行 できる)に迫っていますから, この程度の計算ならパソコンでやれないことも ありません. しかし, 電子計算機の草創期ではそんな余裕なことを言っていら れませんでした. 1960年代になってやっとスーパーコンピューターと呼ばれる 高速計算機が現れてきましたが, それでもたかだか1MFlops(1秒間に100万回の 浮動小数点演算を実行できる)程度の能力しかありませんでした. というわけ で, 当時は最速の計算機を使っても, 上記のDFTの計算が100秒以上かかってし まうという状況であり, $N$が大きい場合のDFTは気軽には実行できませんでし た. そこで, 「DFTを計算するための何かうまい手法はないか?」というのが時代の 要請であり, 「必要は発明の母」の諺どおりに高速フーリエ変換(Fast Fourier Transform = FFT)が生まれることになりました. ここでは, やや冗長にはなりますが, FFTの「発見者」である Cooley and Tukey (1965)の論文に沿ってその概要を見ていきましょう. まず, 変換の長さ$N$が$N=N_1P$と2つの整数$N_1,P$の積になっている場 合を考えます. このとき, (\ref{WDFT})の添字$j,k$それぞれを2つの添字の組 $(j_p,j_1)$, $(k_1,k_p)$によって, \begin{equation} j=N_1j_p+j_1,\quad (j_1=0,1,\cdots,N_1-1; \ j_p=0,1,\cdots,P-1), \end{equation} \begin{equation} k=Pk_1+k_p,\quad (k_1=0,1,\cdots,N_1-1;\ k_p=0,1,\cdots,P-1), \end{equation} のように表示することにします. これは, Cなどのプログラミング言語におい て, $N_1\times P$および$P\times N_1$の2次元配列を使っていると考え ると分かりやすいかもしれません. この表示を使うと, (\ref{WDFT})は以下のように書換えられます: \begin{equation} F(j_p,j_1) = \sum_{k_1=0}^{N_1-1}\sum_{k_p=0}^{P-1} W_N^{(N_1j_p+j_1)(Pk_1+k_p)}A(k_1,k_p). \label{WDFT2} \end{equation} ここで, \begin{eqnarray*} W_N^{(N_1j_p+j_1)(Pk_1+k_p)} &=& W_N^{N_1Pj_pk_1+Pj_1k_1+(N_1j_p+j_1)k_p}\\ &=& W_N^{Nj_pk_1}W_N^{Pj_1k_1}W_N^{(N_1j_p+j_1)k_p}\\ &=& (W_N^N)^{j_pk_1}W_N^{Pj_1k_1}W_N^{(N_1j_p+j_1)k_p}\\ &=& W_N^{Pj_1k_1}W_N^{(N_1j_p+j_1)k_p} \end{eqnarray*} であることに注意すると, (\ref{WDFT2})は, \begin{equation} F(j_p,j_1) = \sum_{k_p=0}^{P-1}W_N^{(N_1j_p+j_1)k_p} \sum_{k_1=0}^{N_1-1}W_N^{Pj_1k_1}A(k_1,k_p), \label{WDFT3} \end{equation} すなわち, \begin{equation} F(j_p,j_1) = \sum_{k_p=0}^{P-1}W_N^{(N_1j_p+j_1)k_p}A_1(j_1,k_p), \quad (j_1=0,1,\cdots,N_1-1;\ j_p=0,1,\cdots,P-1), \label{WDFTD1} \end{equation} \begin{equation} A_1(j_1,k_p) = \sum_{k_1=0}^{N_1-1}W_N^{Pj_1k_1}A(k_1,k_p), \quad (j_1=0,1,\cdots,N_1-1;\ k_p=0,1,\cdots,P-1), \label{WDFTD2} \end{equation} と2段階に分解されます. さて, このように分解すると, 演算回数がどう変化 するか数えてみましょう. まず, (\ref{WDFTD2})には積が$N_1^2×P$回, 和が$N_1(N_1-1)×P$回必要で, (\ref{WDFTD2})については, 積が$P^2×N_1$回, 和が$P(P-1)×N_1$回 となりますから, 合計で, 積が$N_1P(N_1+P)=N(N_1+P)$回, 和が$N_1P(N_1+P-2)=N(N_1+P-2)$回必要ということになります. 従って, ${\cal O}(N^2)$の計算量が${\cal O}(N(N_1+P))$になったことに なります. 最初に$N=10000$の場合の計算量を例にとりましたが, この場合, $N_1=P=100$ととってあげると, $N^2=1$億回の計算量が$N(N_1+P)=200$万 回と, $1/50$に減じられることになります. これでも十分にすばらしい改良なのですが, さらに, $P=N_2Q$と$P$が2つの整 数の積に分解できる場合には, $j_p,k_p$それぞれを2つの添字の組 $(j_q,j_2)$, $(k_2,k_q)$によって, \begin{equation} j_p=N_2j_q+j_2,\quad (j_2=0,1,\cdots,N_2-1;\ j_q=0,1,\cdots,Q-1), \end{equation} \begin{equation} k_p=Qk_2+k_q,\quad (k_2=0,1,\cdots,N_2-1;\ k_q=0,1,\cdots,Q-1), \end{equation} と表すことにより, (\ref{WDFTD1})の計算は以下のように変形できます(以下, 添字の範囲につい ては省略しています): \begin{eqnarray} F(j_q,j_2,j_1) &=& \sum_{k_q=0}^{Q-1}\sum_{k_2=0}^{N_2-1} W_N^{(N_1(N_2j_q+j_2)+j_1)(Qk_2+k_q)}A_1(j_1,k_2,k_q)\\ &=& \sum_{k_q=0}^{Q-1}W_N^{(N_1(N_2j_q+j_2)+j_1)k_q} \sum_{k_2=0}^{N_2-1}W_N^{Q(N_1(N_2j_q+j_2)+j_1)k_2}A_1(j_1,k_2,k_q)\\ &=& \sum_{k_q=0}^{Q-1}W_N^{(N_1(N_2j_q+j_2)+j_1)k_q}A_2(j_1,j_2,k_q). \label{WDFT4} \end{eqnarray} ここに, $W_N^{QN_1N_2}=W_N^N=1$であることに注意して, \begin{equation} A_2(j_1,j_2,k_q)=\sum_{k_2=0}^{N_2-1}W_N^{Q(N_1j_2+j_1)k_2} A_1(j_1,k_2,k_q) \end{equation} と定めました. このように分解すると, 計算量は${\cal O}(N(N_1+N_2+Q))$に下げられること になります. さらに$Q=N_3R$と表されるときに(\ref{WDFT4})$を分解すると \cdots$と,この分解を繰り返していくことを考えると, 結局, $N=N_1N_2\cdots N_m$のとき, $F(j)$は以下のような漸化式的手順で計算できることになります. \begin{enumerate} \item まず, 以下のように初期値$A_1$を$A$から定めます: \begin{equation} \begin{array}{c} Q_1=N_1;\quad P_1=N/Q_1,\\ \displaystyle A_1(j_1,k_p) = \sum_{k_1=0}^{N_1-1}W_N^{P_1j_1k_1}A(k_1,k_p),\\ \quad (j_1=0,1,\cdots,N_1-1;\ k_p=0,1,\cdots,P_1-1). \end{array} \label{FFT1} \end{equation} \item あとは$l=2,3,\cdots,m$について, 漸化式的に順に$A_l$を求めていきます: \begin{equation} \begin{array}{c} Q_l=N_1N_2\cdots N_l;\quad P_l=N/Q_l,\\ \displaystyle A_l(j_1,j_2,\cdots,j_{l-1},j_l,k_p)=\sum_{k_l=0}^{N_l-1} W_N^{P_l(Q_{l-1}j_l+Q_{l-2}j_{l-1}+\cdots+Q_1j_2+j_1)k_l} A_{l-1}(j_1,j_2,\cdots,j_{l-1},k_l,k_p)\\ \quad (j_r=0,1,\cdots,N_r-1\ (r=1,2,\cdots,l); \ k_p=0,1,\cdots,P_l-1). \end{array} \label{FFT2} \end{equation} \item 以上を繰り返していくと最終的に$A_m$が求まるので, それから$F$を添字の並 べ替えによって以下のように求めます: \begin{equation} \begin{array}{c} \displaystyle F(j_m,j_{m-1},\cdots,j_2,j_1)=A_m(j_1,j_2,\cdots,j_{m-1},j_m),\\ \quad (j_r=0,1,\cdots,N_r-1\ (r=1,2,\cdots,m)). \end{array} \label{FFT3} \end{equation} \end{enumerate} これで, 入力$A$から出力$F$を得ることができました. 要される計算量は, 数えてみればすぐ分かりますように, ${\cal O}(N(N_1+N_2+\cdots +N_m))$にまで減らされることになります. 特に, $N$がある整数$r$の羃によって, $N=r^m$と書ける場合には, $N_1=N_2=\cdots =N_m=r$ととると, ${\cal O}(Nrm)={\cal O}(Nr\log_r N)$ と表せます. 最初の$N=10000$の場合について, $r=10$ととると, 計算量は, ${\cal O}(Nr\log_r N)={\cal O}(10^4\times 10\times 4)={\cal O}(400000)$, となり, 結局DFTを直接計算するのに必要な1億回の計算量が, 40万回となり, $1/250$にまで減じられることになります \footnote{もっとも, この場合$10000=(2\times 5)^4=2^4\times 5^4$ とより細かく分解した方が, 計算量が${\cal O}(10^4\times (2\times4 +5\times 4))={\cal O}(280000)$と, さらに少なくなります.}. 以上のアルゴリズムによって, ${\cal O}(N^2)$の計算量が必要だったDFTが ${\cal O}(N\log N)$の計算量でできることが分かります. このアルゴリズムは Cooley and Tukey (1965)によって発表され\footnote{Cooley and Tukey (1965)では, 具体的な表式としては $N=2^m$の場合のみが示されますので, $N$が上記のように一般の$N_r$で分 解される場合についての表式は, その若干の拡張である Bergland (1967)によっ ています.}, Fast Fourier Transform = FFT と呼ばれています. このアルゴリズムを使うことによる計算量の減少は劇的で すが, アルゴリズム自体は非常に簡明なので1965年まで発見されなかったとい うことを不思議に思われる方もいらっしゃるかもしれません. この疑問はもっ ともで, 最初に紹介した部分で「発見者」とかぎ括弧を付けたのもそのためで す. 実際, Cooley and Tukey (1965)以前にも同様のアルゴリズムを独立に発 見し, 実際, プログラミングまで行っていた人もいたのですが, FFTを実際に 世に広めたのはこの5ページの簡潔な論文でしたので, 今ではFFTの発見者の栄 誉は Cooley and Tukey (1965)に与えられています. このあたりのFFTの歴史 については, Brigham ``The Fast Fourier Transform'' (1974) (邦訳: 宮川・今井訳「高速フーリエ変換」(1978)) の第1章の末尾に詳しくまとめられ ています. \section{行列表示によるFFTの理解} 前節で述べた内容でFFTの基本的な構成法は尽されていて, これだでもFFTを計 算するプログラムが作れるはずですが, FFTの手法に慣れないうちは「確かに 数えてみると計算量が劇的に減るのは分かるけど, 添字が入り乱れていて手品 みたい.」という印象を持つ方も多いと思います. そのような場合は, 全体の 計算の流れを行列計算の形に表示しておくと理解が容易になります. $N$が2の 羃の場合についてはいろいろな教科書に載っていますので, ここでは $N=6=3\times 2$の場合について前節のFFTアルゴリズムを行列の形で書き下し てみましょう. まず, DFTの定義式(\ref{WDFT})を$N=6$の場合について書いておきます: \begin{equation} \left[\begin{array}{l} F(0)\\F(1)\\F(2)\\F(3)\\F(4)\\F(5) \end{array}\right] = \left[\begin{array}{llllll} \omega^0 &\omega^0 &\omega^0 &\omega^0 &\omega^0 &\omega^0 \\ \omega^0 &\omega^1 &\omega^2 &\omega^3 &\omega^4 &\omega^5 \\ \omega^0 &\omega^2 &\omega^4 &\omega^6 &\omega^8 &\omega^{10} \\ \omega^0 &\omega^3 &\omega^6 &\omega^9 &\omega^{12} &\omega^{15} \\ \omega^0 &\omega^4 &\omega^8 &\omega^{12} &\omega^{16} &\omega^{20} \\ \omega^0 &\omega^5 &\omega^{10} &\omega^{15} &\omega^{20} &\omega^{25} \\ \end{array}\right] \left[\begin{array}{l} A(0)\\A(1)\\A(2)\\A(3)\\A(4)\\A(5) \end{array}\right]. \label{DFT-M} \end{equation} ここに, $\omega=W_6$としています. さて, 前節で述べたFFTアルゴリズムを$N_1=3,N_2=2$として書き下しみます: \begin{equation} \left[\begin{array}{l} F(0)\\F(1)\\F(2)\\F(3)\\F(4)\\F(5) \end{array}\right] = \left[\begin{array}{llllll} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{array}\right] \left[\begin{array}{llllll} \omega^0 & \omega^0 & 0 & 0 & 0 & 0 \\ \omega^0 & \omega^3 & 0 & 0 & 0 & 0 \\ 0 & 0 & \omega^0 & \omega^1 & 0 & 0 \\ 0 & 0 & \omega^0 & \omega^4 & 0 & 0 \\ 0 & 0 & 0 & 0 & \omega^0 & \omega^2 \\ 0 & 0 & 0 & 0 & \omega^0 & \omega^5 \end{array}\right] \left[\begin{array}{llllll} \omega^0 & 0 & \omega^0 & 0 & \omega^0 & 0 \\ 0 & \omega^0 & 0 & \omega^0 & 0 & \omega^0 \\ \omega^0 & 0 & \omega^2 & 0 & \omega^4 & 0 \\ 0 & \omega^0 & 0 & \omega^2 & 0 & \omega^4 \\ \omega^0 & 0 & \omega^4 & 0 & \omega^8 & 0 \\ 0 & \omega^0 & 0 & \omega^4 & 0 & \omega^8 \end{array}\right] \left[\begin{array}{l} A(0)\\A(1)\\A(2)\\A(3)\\A(4)\\A(5) \end{array}\right]. \label{FFT-M} \end{equation} 両者がちゃんと同じ演算になっていることは, (\ref{FFT-M})の右辺の3つの正 方行列の積を計算し, $\omega^6=\omega^0=1$であることに注意すれば簡単に 確認できます. FFTのアルゴリズムによって計算量が減少するのも(\ref{FFT-M})を見れば一目 瞭然で, (\ref{DFT-M})の右辺の行列が(\ref{FFT-M})では非零成分の少ない疎 行列の積に分解されているため, 演算を右から行うと計算回数が少なくて済む ということに他ならないことが分かります. また, このような行列表示を使うと, FFTアルゴリズムの別の形式についての 見通しも良くなります. 例えば, (\ref{DFT-M})の右辺の行列は対称行列であ るため, 転置をとっても同じ形になります. そこで, (\ref{FFT-M})の右辺に 現れている3つの正方行列の積全体の転置をとっても良いことになり, そうす ると$A$のベクトルに演算される行列の順番が変えられることになります. こうして Cooley and Tukeyのアルゴリズムの``転置''に相当するアルゴリズ ムが導かれますし, また, (\ref{FFT-M})では最後に並べ替えの行列が掛けら れていますが, この並べ替えの操作を途中の行列に埋め込んだ ``self-sorting''型のアルゴリズムも提案されています. これらのFFTのいろ いろな形態については, Temperton (1983)に詳細に整理されています. \section{Faster Fourier Transform?} FFTの基本的アルゴリズムがCooley and Tukey (1965)によって発表されて以来, さらに計算を効率化するための研究が数多く行われました. FFTが計算量${\cal O}(N^2)$を${\cal O}(N\log N)$に減少させたようには劇 的ではありませんが, 場合によっては2倍程度の効率化につながりますので, 以下, その効率化の概略を見ていきましょう. \subsection{$N$が素数の場合のDFT} \label{Winograd} FFTは変換の項数$N$が$N=N_1N_2\cdots N_m$と分解できる場合のDFTの高速化 を図るものであり, $N$が素数の場合にはこの手は使えません. 例えば, $N=5$の場合のDFT: \begin{equation} \left[\begin{array}{l} F(0)\\F(1)\\F(2)\\F(3)\\F(4) \end{array}\right] = \left[\begin{array}{llllll} \omega^0 & \omega^0 & \omega^0 & \omega^0 & \omega^0 \\ \omega^0 & \omega^1 & \omega^2 & \omega^3 & \omega^4 \\ \omega^0 & \omega^2 & \omega^4 & \omega^6 & \omega^8 \\ \omega^0 & \omega^3 & \omega^6 & \omega^9 & \omega^{12} \\ \omega^0 & \omega^4 & \omega^8 & \omega^{12} & \omega^{16} \\ \end{array}\right] \left[\begin{array}{l} A(0)\\A(1)\\A(2)\\A(3)\\A(4) \end{array}\right], \label{DFT-M5} \end{equation} を考えてみましょう. ここで, $\omega=W_5$と略記しています. この場合, $\omega^0=1$なので, この成分との乗算が必要ない ことを考慮しても, この計算には乗算16回, 加算20回はどうしても必要なよう に見えます. しかし, 以下のように計算過程を工夫すると, 少ない演算回数で この計算が可能になります. まず, $\theta=2\pi/5$として5つの定数: \begin{equation} c_1=\frac{\cos\theta+\cos 2\theta}2-1,\ c_2=\frac{\cos\theta-\cos 2\theta}2,\ c_3=i(\sin\theta+\sin2\theta),\ c_4=i\sin2\theta,\ c_5=i(\sin\theta-\sin 2\theta), \end{equation} を前もって計算しておくと, $F(0),\cdots,F(4)$は, \begin{equation} \begin{array}{lllll} s_1=A(1)+A(4),& s_2=A(1)-A(4),& s_3=A(3)+A(2),& s_4=A(3)-A(2),& \\ t_1=s_1+s_3,& t_2=s_1-s_3,& t_3=s_2+s_4,& F(0)=t_1+A(0),& \\ u_1=c_1 t_1,& u_2=c_2 t_2,& u_3=c_3 s_2,& u_4=c_4 t_3,& u_5=c_5 s_4,\\ v_1=F(0)+u_1,& v_2=v_1+u_2,& v_3=v_1-u_2,& v_4=u_3-u_4,& v_5=u_4+u_5,\\ F(1)=v_2+v_4, & F(2)=v_3+v_5,& F(3)=v_3-v_5, & F(4)=v_2-v_4 & \end{array} \label{DFT5} \end{equation} の手順で計算できます. この手順に要する演算回数は, 乗算5回, 加減算17回 と, DFTの定義式どおりの計算よりもかなり少なくなっています. しかも, (\ref{DFT5})に現れる乗算はすべて, (実数または純虚数)$\times$複 素数の形になっているため, 計算機の内部のように複素数を2つの実数の組と して表した場合の乗算回数は, $5\times 2=10$回になりますが, (\ref{DFT-M5})の方はすべて 複素数$\times$複素数の形ですから, 複素数を2つの実数の組として表した場 合の乗算回数は$16\times 4=64$回であり, 改善が大きいことが分かります. このような効率化は, 最初Rader (1968)によって提案されましたが, Winograd (1978)によって様々な$N$に対して系統的に導く方法が整理されまし た. また, この手法はFFTそのものの高速化にも応用できます. FFTのアルゴリズムの中の漸化式(\ref{FFT2})は以下のように2段階に分離でき ます: \begin{eqnarray} A'_{l-1}(j_1,j_2,\cdots,j_{l-1},k_l,k_p) &=& W_N^{P_l(Q_{l-2}j_{l-1}+\cdots+Q_1j_2+j_1)k_l} A_{l-1}(j_1,j_2,\cdots,j_{l-1},k_l,k_p), \label{FFT-ROT}\\ A_l(j_1,j_2,\cdots,j_{l-1},j_l,k_p) &=& \sum_{k_l=0}^{N_l-1} W_{N_l}^{j_lk_l} A'_{l-1}(j_1,j_2,\cdots,j_{l-1},k_l,k_p). \label{FFT-Nl} \end{eqnarray} ここで, $W_N^{P_lQ_{l-1}}=W_N^{N/N_l}=W_{N_l}$であることを用いました. このようにすると, (\ref{FFT-Nl})の方は長さ$N_l$のDFTですから, このDFTの 計算を\ref{Winograd}のように最適化しておけば全体のFFTの高速化につなが ります. これは(\ref{FFT2})の計算には$N\cdot N_l$回の乗算が必要なのに対 して, (\ref{FFT-ROT})は$N$回, (\ref{FFT-Nl})は$(N/N_l)\cdot N_l^2$回の 乗算が必要で, 定義式どおりでは2段階に分けると(\ref{FFT-ROT})の分だけ乗 算回数が増えてしまいますが, \ref{Winograd}の最適化を行っておくと, 後者 の演算における乗算回数がかなり減らせるため効率化につながるわけです. また, このように変形した際の(\ref{FFT-ROT})の右辺の係数は 回転因子(= twiddle factor)と呼ばれています. \subsection{PFA = Prime Factor Algorithm} 通常のFFTのアルゴリズムでは回転因子の存在は不可避なのですが, 回転因子を伴わないアルゴリズムも存在します. 変換の長さ$N$が$N=N_1N_2$と2つの整数の積となっている場合を考えます. ただし, FFTの場合と異なり, $N_1$と$N_2$とは互いに素であるという条件を 付けておきます. まず, $N_1$と$N_2$とは互いに素であるという条件より, \begin{equation} p_1N_2=p_2N_1+1, \quad q_2N_1=q_1N_2+1, \quad (0