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

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

基本事項
- 回転因子W
- 離散フーリエ変換を周波数の偶数番と奇数番に分けると
N=16でFFTの計算式を確認
ここまでの計算で何がやりたかったのか
まとめ
過去の記事

基本事項

回転因子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)$ との対応を考えてみてください．