【Rで高速フーリエ変換】バタフライダイアグラムを読む - ケィオスの時系列解析メモランダム

私は肉まんの白い部分が好きです．ほのかに甘い，あの白い部分だけ食べたいです．

　ということで，今回は最近のFFTのお話
chaos-kiyono.hatenablog.com
の続きです．

　FFTの説明で，バタフライ演算ってのが登場します．バタフライ演算は下の図のようなバタフライダイアグラムで表すことができます．

「バタフライ」ってのは，FFTの計算を図で表すと蝶蝶の羽のように見えることに由来していると思います (知らんけど)．今回はバタフライダイアグラムが表す計算過程を説明します．

右上がりは足し，右下がりは引く

　時系列の長さ $N$ が4の時系列 $\{x[0], x[1], x[2], x[3]\}$ を考えます．前回はRでの実装を想定して，時系列の順番について1から番号をふりましたが，今回は0からです．ビット反転を説明するために，一般的な表現を使います．

　 $N=4$ のとき，バタフライダイアグラムに補足の記号 (プラスマイナスと回転因子)を加えると下図のようになります．

　上の図では，左から右へ計算が進んでいきます．赤 (右上がり線)と青 (右下がり線)で描いた線は，それぞれ，「足し算」と「引き算」を表します．

　前回，長さ $N$ が2の倍数であれば，時系列 $\{x[0], x[1], \cdots, x[N-1]\}$ の離散フーリエ変換

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

は，偶数番目の周波数の離散フーリエ変換

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

と，奇数番目の周波数の離散フーリエ変換

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

に分けられることを説明しました (詳しくは説明してないけど)．ここで， $W_{N}=e^{-{\bf i} 2 \pi /N}$ です．

　バタフライダイアグラムの上半分は偶数番目の周波数の離散フーリエ変換，下半分は奇数番目の周波数の離散フーリエ変換に対応します．

　ダイアグラムの見方を理解するために，以下ではバタフライダイアグラムのパーツに注目してみます．

右上がりの線 (足し算)

　下の図のように，右上がりの線 (赤)は，●でつないだ2本の線に対応する値を足すことを表しています．

　これは，偶数番目の周波数の離散フーリエ変換

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

に登場する

　 $\left(x[i]+x[N/2+i] \right)$

の計算を表しています．

右下がりの線 (引いて回転演算子をかける)

　下の図のように右下がりの線 (青)は，●でつないだ2本の線に対応する値を引くことを表しています．さらに，下の●の近くに $W_4^1$ と書いてあるので，それをかけます．

　これは，奇数番目の周波数の離散フーリエ変換

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

に登場する

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

の計算を表しています．

回転演算子

　回転演算 $W_N$ をかける必要があるのは，奇数番目の周波数の離散フーリエ変換です．つまり，バタフライダイアグラムでは下半分の部分です．かける必要がある回転演算子 $W_{N}^{i}$ は， $i=0 ,1 , \cdots, N/2-1$ の場合です．したがって， $N=4$ のときは，上から下に，回転演算子を2つ ( $N/2=2$ )，つまり，

　 $W_{4}^{0}, W_{4}^{1}$

をかけます．図に書き込めば，下のようになります．

ビット反転でフーリエ変換の結果を並び替え

　以上の説明で，バタフライダイアグラムが表している計算アルゴリズムを理解できたでしょうか．例えば，下の図のような， $N=8$ のバタフライダイアグラムが表す計算式をイメージできますか．

　とりあえず，バタフライダイアグラムの読み方は分かったとしても，右端に書いた離散フーリエ変換の結果 $X(k)$ の順番がめちゃくちゃです．この順番は，どうやって決められるのでしょうか．そこで，登場するのがビット反転 (ビットリバース)です．

ビット反転の計算

　ビット反転では，10進数を２進数に変換して，２進数の並びを反転してから10進数に戻します．ただし，時系列の長さが $N$ のときは， $\log_2 N$ 桁の２進数で表し，上位の桁がないときは，0を詰めます (わかった？)．

　例えば， $N=8$ のとき， $\log_2 8 = 3$ 桁の2進数を考えます． $N=8$ のとき，1をビット反転してみます．10進数の $1_{(10)}$ は，２進数で表しても $1_{(2)}$ ですが，3桁にしたいので，0を詰めて，

　 $001_{(2)}$

とします．数字の順番を逆にならべると，

　 $100_{(2)}$

になります．これを10進数にすると，

　 $100_{(2)} = 1\times 2^2 + 0 \times 2^1 + 0 \times 2^0 = 4$

になります．確かに，上の図では， $x[1]$ の右端に， $X(4)$ が書いてあります．

なぜ，ビット反転でうまくいくの？

　ビット反転で対応する周波数の順番が決まる理由を考えてみます． $N=8$ の例を表した下の図を見てください．この図を使って考えてみます．

　まず，バタフライダイアグラムの最初の段 (一番左)で，上半分にあれば，それは偶数番目の周波数です．だって，

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

の計算をしているからです． $N=8$ の例では，周波数の番号は $0, 2, 4, 6$ のどれかです．逆に，下半分は奇数で， $1, 3, 5, 7$ のどれかです．

　偶数ってことは，２進数で表したとき1の位は0です．奇数ってことは，２進数で表したとき，1の位は1です．この桁は，左端の $x[i]$ の順番 $i$ を２進数で表したときの，100の位に一致します．

　後は計算式を追ってやればビット反転との対応が見えてきます．

　例えば， $X(4)$ では， $4$ は偶数なので，最初は偶数のフーリエ変換の式

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

を使います．これで，2進数で表したときの1桁目は，0に確定です．この式では，周波数が $2k$ のフーリエ変換が，半分の $k$ でのフーリエ変換に書き換わっています．

　ということは， $X(4)$ の $4$ を2で割って， $X(2)$ のフーリエ変換の式を次の段階で使うことになります． $2$ は偶数なので，バタフライダイアグラムの真ん中の段は，上半分のパスを通ります．なので，2進数で表したときの10の位は，0になります．

　最後は， $X(2)$ の $2$ を2で割って， $X(1)$ のフーリエ変換の式を使うことになります．1は奇数なので，バタフライダイアグラムの一番右側の段は，下半分のパスを通ります．なので，2進数で表したときの100の位は，1になります．

　順番を整理すると，0→0→1です．これは，4を2進数に変換した $100_{2}$ を逆に並べたものと同じです．

　エビフライダイアグラムの左側の $x[i]$ の順番 $i$ は，上の図に書いた0と1を左から右へ並べた2進数になってます．一方，右側の $X(k)$ の順番 $k$ (フーリエ変換結果)では，左が1の位に対応するので，上の図に書いた0と1を右から左へ並べた2進数になってます．

　こんな説明で，理解できますでしょうか．

まとめ

　Rのコマンドとか，パッケージを使って，それっぽい結果出すだけでなく，使う道具にも興味や疑問をもって基礎知識を身につけてください．