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

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

【Rで高速フーリエ変換】2のべき乗の長さの時系列用FFTの実装についての解説

R言語 時系列解析

前回，2のべき乗の長さの時系列を，高速に離散フーリエ変換 (FFT)できるアルゴリズムのRスクリプトを作りました (高速じゃないけど)．
chaos-kiyono.hatenablog.com

　このFFTアルゴリズムは，Cooley-Tukey型と呼ばれるそうです (知らんけど)．今回はこのアルゴリズムについて少し解説します．

要は周波数について奇数番と偶数番に分けるだけ

　以下では，長さが2のべき乗の時系列\{x[1], x[2], \cdots, x[N]\}を離散フーリエ変換することを考えます．時系列解析の本では，データの番号を0からN-1で考えることが多いですが，Rでは，配列の番号が1からはじまるので，データの番号は1からNにしました．他の教科書とかインターネットの解説と番号がずれているので注意してください．

　このとき，離散フーリエ変換 X(k) は，

　\displaystyle X(k) = \sum_{i = 1}^{N} x[i]\, e^{-{\bf i} 2 \pi k (i-1)/N }

です．kは，0からN-1までです．サンプリング間隔が1のときは周波数は，f_k = k/Nです．

　Cooley-Tukey型FFTって，要は，kについて，奇数番と偶数番を分けるだけです．

N=4の例

　簡単な例として，N=4でやってみます (速くないけど)．kが偶数で，k=2のとき (eの肩の数値に注目)，

　\begin{eqnarray} X(k=2) &=& \sum_{i = 1}^{4} x[i] \, e^{-{\bf i} 2 \pi \cdot 2 (i-1)/N }\\ &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 2 \cdot 1/4} + x[3]\, e^{-{\bf i} 2 \pi \cdot 2 \cdot 2/4}+ x[4]\, e^{-{\bf i} 2 \pi \cdot 2 \cdot 3/4} \\ &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 1/2} + x[3]\, e^{-{\bf i} 2 \pi \cdot 4/4}+ x[4]\, e^{-{\bf i} 2 \pi \cdot 6/4} \\ &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 1/2} + x[3]+ x[4]\, e^{-{\bf i} 2 \pi \cdot 1/2} \\ &=& \left(x[1] + x[3] \right) + \left( x[2]+ x[4]\right) e^{-{\bf i} 2 \pi \cdot 1/2} \\ &=& \sum_{i = 1}^{2} \left(x[i]+x[i+2] \right) \, e^{-{\bf i} 2 \pi \cdot (i-1)/2} \end{eqnarray}

と変形できます．この変形の注目点は

  • 最初は\displaystyle \sum_{i = 1}^{4}で4項の和だった部分が，変形すると\displaystyle \sum_{i = 1}^{2}になり，半分の2項の和になった．
  • e^{-{\bf i} 2 \pi k (i-1)/N}の部分のk=2が半分の1になった (周波数が半分になった)．
  • e^{-{\bf i} 2 \pi k (i-1)/N}の部分のN=4が半分の2になった．

ということです．

　また，kが奇数で，k=1のとき，

　\begin{eqnarray} X (k=1) &=& \sum_{i = 1}^{4} x[i] \, e^{-{\bf i} 2 \pi \cdot 1 \cdot (i-1)/N }\\ &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 1 \cdot 1/4} + x[3]\, e^{-{\bf i} 2 \pi \cdot 1 \cdot 2/4}+ x[4]\, e^{-{\bf i} 2 \pi \cdot 1 \cdot 3/4} \\ &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 1/4} + x[3]\, e^{-{\bf i} \pi}+ x[4]\, e^{-{\bf i} 2 \pi \cdot 3/4} \\ &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 1/4} - x[3] - x[4]\, e^{-{\bf i} 2 \pi \cdot 1/4} \\ &=& \left(x[1] - x[3] \right) + \left( x[2]- x[4]\right) e^{-{\bf i} 2 \pi \cdot 1/4} \\ &=& \sum_{i = 1}^{2} \left\{\left(x[i]-x[i+2] \right) \, e^{-{\bf i} 2 \pi \cdot (i-1)/4}\right\} e^{-{\bf i} 2 \pi \cdot 0 \cdot (i-1)/2} \end{eqnarray}

と変形できます (最後の式は結果を知っているので無理矢理書きました)．やはり，\displaystyle \sum_{i = 1}^{4}で4項の和だった部分を，変形して\displaystyle \sum_{i = 1}^{2}にできました．

　一般にどんな関係式が成り立つか考えてみてください．

一般形

　Nが2の倍数のとき，離散フーリエ変換

　\displaystyle X(k) = \sum_{i = 1}^{N} x[i]\, e^{-{\bf i} 2 \pi k (i-1)/N }

について，以下の関係が成り立ちます．

　\begin{eqnarray} X(2k) &=& \sum_{i = 1}^{N/2} \left(x[i]+x[N/2+i] \right) e^{-{\bf i} 2 \pi k \cdot (i-1)/(N/2)} \\ X[2k+1] &=& \sum_{i = 1}^{N/2}\left\{ \left(x[i]-x[N/2+i] \right) e^{-{\bf i} 2 \pi (i-1)/N} \right\} e^{-{\bf i} 2 \pi k \cdot (i-1)/(N/2)} \end{eqnarray}

すっきり書くために，W_{N}=e^{-{\bf i} 2 \pi /N}とすれば，

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

となります．

　Nが2のべき乗の長さのとき，この分割を繰り返し行えるので，

　x[i]+x[N/2+i]

と，

　\left(x[i]-x[N/2+i] \right) W_{N}^{i-1}

を繰り返し計算すれば離散フーリエ変換できちゃいます．これが，Cooley-Tukey型FFTです．

Rスクリプトとの対応

　前回，RでCooley-Tukey型FFTのスクリプトを書いてみました．
chaos-kiyono.hatenablog.com

　このスクリプトの中の

  # 要素の半分をまとめて指定：1, 2, 3, ..., N/2
  ii <- 1:(N/2)
  # バタフライ演算の右上がり線 (/)に対応する計算
  x1 <- x[ii]+x[ii+N/2] # 足すだけ
  # バタフライ演算の右下がり線 (\)に対応する計算
  x2 <- (x[ii]-x[ii+N/2])*W[ii] # 引いて回転因子をかける

の部分が，

　x[i]+x[N/2+i]

と，

　\left(x[i]-x[N/2+i] \right) W_{N}^{i-1}

の計算です．

次回予告：バタフライダイアグラムの読み方

　FFTのポイントはこれだけですが，次回，バタフライ演算の図の読み方と，ビット反転について説明します．

　未来から帰ってきました．つづきの記事は，これです．
chaos-kiyono.hatenablog.com