2の基数のFFTの実装上の注意

昨日、一番基本的な2の基数のFFTのアルゴリズムをまとめてみた。式で書くと結構簡単に見えるけどそれを実装するとなるとそれなりに面倒なことがある。

実装上の問題

再帰呼び出しとループ

周波数間引きと時間間引きのどちらの場合も、入れ子はrecursive call（再帰呼び出し）で書くと簡単になるけど、どちらも一番内側の呼び出しはふたつのデータに対してスタックにつんでふたつで計算してスタックを戻して、ということをN/2回やることになるので関数呼び出しのオーバーヘッドがバカにならない。ふつうはrecursive callをループに展開して書く。そういうのってすぐ自分の足に穴があくので気をつけなければいけない。FFTの場合ちゃんと動作するものがいっぱいあるので答え合わせをしながらデバグできる。今回はループに展開するところまではやらないつもり。

ビット逆順

もっと問題なのはデータの入れ替え。

ひとつのステップで入力データを前半と後半にわけて、ふたつのFFTにするけど、それを偶数番目と奇数番目の結果としないといけない。これが入れ子になるのでわけがわからなくなる。

普通はどうするか、というとひとつのステップは入れ替えを無視してデータを上書きししてしまって、あとから並べ替えるということをする。

ひとつのステップで

としてこれを入れ子にしたとする。ここで

で

はk番目の結果をn番目の入力位置に上書きすることを表している。

Nが2のベキでこれを入れ子にして繰り返すと規則正しくぐちゃぐちゃに入れ替わる。この入れ替えのことをビット逆順と呼ぶ。場所を2進で表してその上下位を逆に読んだのと同じになる。古いDSPではこれをハードで実行するようになっているものがあった。今でもあるかどうかは知らない。

ビット逆順の話題はいろいろなサイトにあがっているので省略する。今回は2次元に対応するので同じNに対して何度もFFTをすることになる。こういう場合は単純にビット逆順のテーブルを持つのが効率もいいし簡単。