ケィオスの時系列解析メモランダム

時系列解析,生体情報学,数学・物理などの解説です.

【Rで高速フーリエ変換】N=16でCooley-Tukey型FFTを考えてみる

FFTを理解してもらいたい,というか,自分がわかりやすく説明できるようになりたいので,しつこくFFTの話をやっていきます.今回は,時系列の長さN=16でやっていきます.

基本事項

回転因子W

 フーリエ変換の式の表現をシンプルにするために,回転因子W_Nを定義します.

 \displaystyle 
 W_N=e^{-{\bf i} 2 \pi /N}

オイラーの式で,

 \begin{eqnarray} e^{-{\bf i} 2 \pi /N} &=& \cos\left(- \frac{2 \pi}{N} \right) + {\bf i} \sin\left(- \frac{2 \pi}{N} \right) \\
 &=& \cos\left(\frac{2 \pi}{N}\right) - {\bf i} \sin\left( \frac{2 \pi}{N} \right)
\end{eqnarray}

となるので,

 \begin{eqnarray} W_N^{k} = e^{-{\bf i} 2 \pi k/N} &=&  \cos\left(\frac{2 \pi k}{N}\right) - {\bf i} \sin\left( \frac{2 \pi k}{N} \right)
\end{eqnarray}

です.ということで,

 \displaystyle 
 W_N^{0}= 1 \displaystyle 
 W_N^{N}= 1 \displaystyle 
 W_N^{N/2}= -1

となります.

離散フーリエ変換を周波数の偶数番と奇数番に分けると

 時系列を\{x[0], x[1], x[2], x[3], \cdots, x[N-1]\}とすれば,離散フーリエ変換は,

 \displaystyle X(k) = \sum_{i = 0}^{N-1} x[i] \, W_N^{ki}

です.時系列の長さが2の倍数のとき,前に説明したように,kについて,偶数番 (k = 0, 2, \cdots, N/2)と奇数番 (k = 1, 3, \cdots, N/2-1)に分けて以下の式が成り立ちます.

\kappa = 0, 1, 2, \cdots, N/2-1として,
偶数番 k=2\kappaでは,

 \begin{eqnarray}  X(2\kappa) &=& \sum_{i = 0}^{N/2-1} \left(x[i]+x[N/2+i] \right) W_{N/2}^{\kappa i}
 \end{eqnarray}

奇数番 k=2\kappa+1では,

 \begin{eqnarray}
X(2\kappa+1) &=& \sum_{i = 0}^{N/2-1}\left\{ \left(x[i]-x[N/2+i] \right) W_{N}^{i} \right\} W_{N/2}^{\kappa i}
 \end{eqnarray}

この式の導出については,詳しく説明してあるページがいくつもあるので,検索して調べてください.

N=16でFFTの計算式を確認

 N=16=2^4で,時系列を\{x[0], x[1], x[2], x[3], \cdots, x[15]\}として,FFTの計算過程を見ていきます.

1回目の分割

 周波数の偶数番 k=0, 2, \cdots, 14では,

 \begin{eqnarray}  X(0) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{0 i} \\
 X(2) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{i} \\
 X(4) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{2 i} \\
 X(6) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{3 i} \\
 X(8) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{4 i} \\
 X(10) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{5 i} \\
 X(12) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{6 i} \\
 X(14) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{7 i} \\
 \end{eqnarray}

です.ここで,\left(x[i]+x[8+i] \right)の計算が共通なので,

 x_8^{(0)}[i] = x[i]+x[8+i]

として,i=0, 1, \cdots, 7について足し算すれば,全体で足し算する回数が減ります.xに添え字がたくさんつきますが,計算の確認のため記号を付けておきます.上付き添え字の{}^{(0)}は偶数番周波数の置き換え (足し算),上付き添え字の{}^{(1)}は奇数番周波数の置き換え (引き算&回転因子)を表すことにします.

 周波数の奇数番 k=1, 3, \cdots, 15では,

 \begin{eqnarray}  X(1) &=& \sum_{i = 0}^{7}\left\{ \left(x[i]-x[8+i] \right) W_{16}^{i} \right\} W_{8}^{0 i} \\
 X(3) &=& \sum_{i = 0}^{7}\left\{ \left(x[i]-x[8+i] \right) W_{16}^{i} \right\} W_{8}^{1 i} \\
 X(5) &=& \sum_{i = 0}^{7}\left\{ \left(x[i]-x[8+i] \right) W_{16}^{i} \right\} W_{8}^{2 i} \\
&\vdots&\\
 X(15) &=& \sum_{i = 0}^{7}\left\{ \left(x[i]-x[8+i] \right) W_{16}^{i} \right\} W_{8}^{7 i} \\
 \end{eqnarray}

です.ここでは,\left(x[i]-x[8+i] \right) W_{16}^{i}の計算が共通なので,

 x_8^{(1)}[i] = \left(x[i]-x[8+i] \right) W_{16}^{i}

とします.

 上半分に偶数番目の周波数,下半分に奇数番目の周波数をならべて書くと,

 \begin{eqnarray}  X(0) &=& \sum_{i = 0}^{7} x_8^{(0)}[i]\, W_{8}^{0 i} \\
 X(2) &=& \sum_{i = 0}^{7} x_8^{(0)}[i]\,  W_{8}^{1 i} \\
 X(4) &=& \sum_{i = 0}^{7} x_8^{(0)}[i] \, W_{8}^{2 i} \\
 X(6) &=& \sum_{i = 0}^{7} x_8^{(0)}[i] \, W_{8}^{3 i} \\
&\vdots& \\
 X(14) &=& \sum_{i = 0}^{7} x_8^{(0)}[i] \, W_{8}^{7 i} \\
X(1) &=& \sum_{i = 0}^{7} x_8^{(1)}[i]\, W_{8}^{0 i} \\
 X(3) &=& \sum_{i = 0}^{7} x_8^{(1)}[i]\,  W_{8}^{1i} \\
 X(5) &=& \sum_{i = 0}^{7} x_8^{(1)}[i] \, W_{8}^{2 i} \\
 X(7) &=& \sum_{i = 0}^{7} x_8^{(1)}[i] \, W_{8}^{3 i} \\
&\vdots& \\
 X(15) &=& \sum_{i = 0}^{7} x_8^{(1)}[i] \, W_{8}^{7 i} \\
 \end{eqnarray}

となります.ここで,右辺について気づくことは,上半分と下半分のそれぞれで,N=8のときのフーリエ変換になっていると言うことです.このフーリエ変換

 \begin{eqnarray} X_8^{(0)} (k) &=& \sum_{i = 0}^{7} x_8^{(0)}[i]\, W_{8}^{k i}\\
X_8^{(1)} (k) &=& \sum_{i = 0}^{7} x_8^{(1)}[i]\, W_{8}^{k i} \end{eqnarray}

と表し,さらに,偶数番,奇数番で分割すると,

 \begin{eqnarray}  X(0) &=&X_8^{(0)}(0) = \sum_{i = 0}^{3} \left(x_8^{(0)}[i]+x_8^{(0)}[4+i] \right)\, W_{4}^{0 i} \\
 X(2) &=& X_8^{(0)}(1) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(0)}[i]-x_8^{(0)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{0 i}\\
 X(4) &=& X_8^{(0)}(2) =  \sum_{i = 0}^{3} \left(x_8^{(0)}[i]+x_8^{(0)}[4+i] \right)\, W_{4}^{1 i} \\
 X(6) &=& X_8^{(0)}(3) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(0)}[i]-x_8^{(0)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{1 i}\\
&\vdots&\\
 X(14) &=& X_8^{(0)}(7) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(0)}[i]-x_8^{(0)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{3 i}\\
X(1) &=&X_8^{(1)}(0) = \sum_{i = 0}^{3} \left(x_8^{(1)}[i]+x_8^{(1)}[4+i] \right)\, W_{4}^{0 i} \\
 X(3) &=& X_8^{(1)}(1) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(1)}[i]-x_8^{(1)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{0 i}\\
 X(5) &=& X_8^{(1)}(2) =  \sum_{i = 0}^{3} \left(x_8^{(1)}[i]+x_8^{(1)}[4+i] \right)\, W_{4}^{1 i} \\
 X(7) &=& X_8^{(1)}(3) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(1)}[i]-x_8^{(1)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{1 i}\\
&\vdots&\\
 X(15) &=& X_8^{(1)}(7) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(1)}[i]-x_8^{(1)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{3 i}
 \end{eqnarray}

となります.

2回目の分割

 上の式にも共通の部分があることに気づくので,

 \begin{eqnarray}  x_4^{(00)}[i] &=& x_8^{(0)}[i]+x_8^{(0)}[4+i] \\
x_4^{(01)}[i] &=& \left(x_8^{(0)}[i]-x_8^{(0)}[4+i] \right) W_{8}^{i}\\
x_4^{(10)}[i] &=& x_8^{(1)}[i]+x_8^{(1)}[4+i] \\
x_4^{(11)}[i] &=& \left(x_8^{(1)}[i]-x_8^{(1)}[4+i] \right) W_{8}^{i}
 \end{eqnarray}

と置きます.N=4フーリエ変換を,

 \begin{eqnarray}  X_4^{(00)} (k) &=& \sum_{i = 0}^{3} x_4^{(00)}[i]\, W_{4}^{k i}\\
X_4^{(01)} (k) &=& \sum_{i = 0}^{3} x_4^{(01)} [i]\, W_{4}^{k i}\\
X_4^{(10)} (k) &=& \sum_{i = 0}^{3} x_4^{(10)}[i]\, W_{4}^{k i}\\
X_4^{(11)} (k) &=& \sum_{i = 0}^{3} x_4^{(11)} [i]\, W_{4}^{k i}
 \end{eqnarray}

と表し,またまた,偶数番,奇数番で分割すると,

 \begin{eqnarray}  X(0) &=&X_8^{(0)}(0) = X_4^{(00)}(0)  = \sum_{i = 0}^{1} \left(x_4^{(00)}[i]+x_4^{(00)}[2+i] \right)\, W_{2}^{0 i} \\
 X(2) &=& X_8^{(0)}(1) = X_4^{(01)}(0)  = \sum_{i = 0}^{1} \left(x_4^{(01)}[i]+x_4^{(01)}[2+i] \right) W_{2}^{0 i}\\
 X(4) &=& X_8^{(0)}(2) = X_4^{(00)}(1)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(00)}[i]-x_4^{(00)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{0 i}\\
 X(6) &=& X_8^{(0)}(3) = X_4^{(01)}(1)  = \sum_{i = 0}^{1} \left\{ \left(x_4^{(01)}[i]-x_4^{(01)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{0 i}\\
 X(8) &=& X_8^{(0)}(4) = X_4^{(00)}(2)  = \sum_{i = 0}^{1} \left(x_4^{(00)}[i]+x_4^{(00)}[2+i] \right)\, W_{2}^{1 i} \\
 X(10) &=&X_8^{(0)}(5) = X_4^{(01)}(2)  = \sum_{i = 0}^{1} \left(x_4^{(01)}[i]+x_4^{(01)}[2+i] \right)\, W_{2}^{1 i} \\
 X(12) &=& X_8^{(0)}(6) = X_4^{(00)}(3)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(00)}[i]-x_4^{(00)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{1 i}\\
 X(14) &=& X_8^{(0)}(7) = X_4^{(01)}(3)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(01)}[i]-x_4^{(01)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{1 i}\\
X(1) &=&X_8^{(1)}(0) = X_4^{(10)}(0)  = \sum_{i = 0}^{1} \left(x_4^{(10)}[i]+x_4^{(10)}[2+i] \right)\, W_{2}^{0 i} \\
 X(3) &=& X_8^{(1)}(1) = X_4^{(11)}(0)  = \sum_{i = 0}^{1} \left(x_4^{(11)}[i]+x_4^{(11)}[2+i] \right) W_{2}^{0 i}\\
 X(5) &=& X_8^{(1)}(2) = X_4^{(10)}(1)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(10)}[i]-x_4^{(10)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{0 i}\\
 X(7) &=& X_8^{(1)}(3) = X_4^{(11)}(1)  = \sum_{i = 0}^{1} \left\{ \left(x_4^{(11)}[i]-x_4^{(11)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{0 i}\\
 X(9) &=& X_8^{(1)}(4) = X_4^{(10)}(2)  = \sum_{i = 0}^{1} \left(x_4^{(10)}[i]+x_4^{(10)}[2+i] \right)\, W_{2}^{1 i} \\
 X(11) &=&X_8^{(1)}(5) = X_4^{(11)}(2)  = \sum_{i = 0}^{1} \left(x_4^{(11)}[i]+x_4^{(11)}[2+i] \right)\, W_{2}^{1 i} \\
 X(13) &=& X_8^{(1)}(6) = X_4^{(10)}(3)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(10)}[i]-x_4^{(10)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{1 i}\\
 X(15) &=& X_8^{(1)}(7) = X_4^{(11)}(3)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(11)}[i]-x_4^{(11)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{1 i}\\
 \end{eqnarray}

となります.

3回目と4回目の分割

 ちょっと込み入ってきたので,自信がなくなってきましたが,ここでも,共通部分を,

 \begin{eqnarray}  x_2^{(000)}[i] &=& \left(x_4^{(00)}[i]+x_4^{(00)}[2+i] \right)\\
x_2^{(001)}[i] &=& \left(x_4^{(00)}[i]+x_4^{(00)}[2+i] \right) W_4^i\\
x_2^{(010)}[i] &=& \left(x_4^{(01)}[i]-x_4^{(01)}[2+i] \right) \\
x_2^{(011)}[i] &=& \left(x_4^{(01)}[i]-x_4^{(01)}[2+i] \right) W_4^i\\
x_2^{(100)}[i] &=& \left(x_4^{(10)}[i]+x_4^{(10)}[2+i] \right)\\
x_2^{(101)}[i] &=& \left(x_4^{(10)}[i]+x_4^{(10)}[2+i] \right) W_4^i\\
x_2^{(110)}[i] &=& \left(x_4^{(11)}[i]-x_4^{(11)}[2+i] \right) \\
x_2^{(111)}[i] &=& \left(x_4^{(11)}[i]-x_4^{(11)}[2+i] \right) W_4^i
 \end{eqnarray}

と置きます.さらにさらに,N=2フーリエ変換

 \begin{eqnarray}  X_2^{(000)} (k) &=& \sum_{i = 0}^{1} x_2^{(000)}[i]\, W_{2}^{k i}\\
X_2^{(001)} (k) &=& \sum_{i = 0}^{1} x_2^{(001)} [i]\, W_{2}^{k i} \\
X_2^{(010)} (k) &=& \sum_{i = 0}^{1} x_2^{(010)}[i]\, W_{2}^{k i}\\
&\vdots&\\
X_2^{(111)} (k) &=& \sum_{i = 0}^{1} x_2^{(111)} [i]\, W_{2}^{k i}
 \end{eqnarray}

と表せば,

 \begin{eqnarray}  X(0) &=& X_8^{(0)}(0) = X_4^{(00)}(0) = X_2^{(000)} (0)  \\
&=&  x_2^{(000)}[0]\, W_{2}^{0} + x_2^{(000)}[1]\, W_{2}^{0} = x_2^{(000)}[0]+ x_2^{(000)}[1]\\
 X(2) &=& X_8^{(0)}(1) = X_4^{(01)}(0)  =  X_2^{(010)} (0)\\
&=&  x_2^{(010)}[0]\, W_{2}^{0} + x_2^{(010)}[1]\, W_{2}^{0} = x_2^{(010)}[0]+ x_2^{(010)}[1]\\
&\vdots& \\
 \end{eqnarray}

のようになります.この変形をして,さらに,これまでと同様のルールでxに添え字をつけて記述すれば,

 \begin{eqnarray}  X(0) &=& x_2^{(000)}[0]+ x_2^{(000)}[1] = x_1^{(0000)}[0]\\
X(1) &=& x_2^{(100)}[0]+ x_2^{(100)}[1] = x_1^{(1000)}[0]\\
X(2) &=& x_2^{(010)}[0]+ x_2^{(010)}[1] = x_1^{(0100)}[0]\\
 X(3) &=& x_2^{(110)}[0]+ x_2^{(110)}[1] = x_1^{(1100)}[0]\\
X(4) &=& x_2^{(001)}[0]+ x_2^{(001)}[1] = x_1^{(0010)}[0]\\
X(5) &=& x_2^{(101)}[0]+ x_2^{(101)}[1] = x_1^{(1010)}[0]\\
X(6) &=& x_2^{(011)}[0]+ x_2^{(011)}[1] = x_1^{(0110)}[0]\\
X(7) &=& x_2^{(111)}[0]+ x_2^{(111)}[1] = x_1^{(1110)}[0]\\
X(8) &=& x_2^{(000)}[0]- x_2^{(000)}[1] = x_1^{(0001)}[0]\\
X(9) &=& x_2^{(100)}[0]- x_2^{(100)}[1]  = x_1^{(1001)}[0] \\
X(10) &=& x_2^{(010)}[0]- x_2^{(010)}[1] = x_1^{(0101)}[0]  \\
X(11) &=& x_2^{(110)}[0]- x_2^{(110)}[1]   = x_1^{(1101)}[0]\\
X(12) &=& x_2^{(001)}[0]- x_2^{(001)}[1] = x_1^{(0011)}[0]  \\
X(13) &=& x_2^{(101)}[0]- x_2^{(101)}[1]  = x_1^{(1011)}[0] \\
X(14) &=& x_2^{(011)}[0]- x_2^{(011)}[1] = x_1^{(0111)}[0]  \\
X(15) &=& x_2^{(111)}[0]- x_2^{(111)}[1] = x_1^{(1111)}[0]
 \end{eqnarray}

となりました.最後は,順番に並べて書きました.

ここまでの計算で何がやりたかったのか

 Cooley-Tukey型FFTで何をやっているかと言えば,周波数の偶数番と奇数番を分割を繰り返しているだけです.この分割というのは,それぞれ,

  • 足し算を計算する x[i]+x[N/2+i]
  • 引き算して回転因子をかける  \left(x[i]-x[N/2+i] \right) W_{N}^{i}

に対応します.そして,最後まで,この2パターンの計算を繰り返しました.

 上の変形では最終的に,

 \begin{eqnarray}  X(0) &=& x_1^{(0000)}[0]\\
X(1) &=& x_1^{(1000)}[0]\\
 X(2) &=& x_1^{(0100)}[0]\\
 X(3) &=& x_1^{(1100)}[0]\\
X(4) &=& x_1^{(0010)}[0]\\
&\vdots&\\
X(14) &=& x_1^{(0111)}[0]  \\
X(15) &=& x_1^{(1111)}[0]
 \end{eqnarray}

という表現がえられました.

 xの上付き添え字の部分(0000), (0100), \cdotsは,周波数の偶数番と奇数番の分割で,どちら側だったかという順番を表しています.

 例えば,X(0)に登場する(0000)は,1回目偶数(0),2回目偶数(0),3回目偶数(0),4回目偶数(0)だったということです.

 X(5)に登場する(1010)は,1回目奇数(1),2回目偶数(0),3回目奇数(1),4回目偶数(0)だったということです.

 そして,X(5)に登場する(1010)は,5を2進数にする計算と同じことをしています.ただし,0, 1の値の並べ方が右から左ではなく,左から右になります.これが,なぜ,ビット反転をするとうまくいくのかの理由です.

 今回のお話は,前回の記事で紹介したバタフライダイアグラムを,式で表してみたということです.下の図のようなバタフライダイアグラムの構造と今回の結果の関係が見えてくるでしょうか.図の中に書いた0, 1と,x^{(\cdots)}の上付き添え字の部分(\cdots)との対応を考えてみてください.

N=16のバタフライダイアグラム.マイナスの記号や回転因子は省略してあります.その代わり,周波数の偶数番と奇数番に対応する演算の部分に,それぞれ,0と1を書きました.

まとめ

 Cooley-Tukey型FFTについては,いろいろな説明の仕方があります.例えば,行列を使った説明があります.私は,行列を使った説明はわかりやすいと思います.

しかし,私は個性的な説明を試みてみたかったので,今回は,周波数の分割をすべて書いてみるというアプローチをしました.数式が嫌いな人が多いので,わかりにくくなったかもしれません.私としては,数式からイメージがわいてくるので,知識の整理には役立ちました.