FFTを固定小数点化する(1)
ネットで便利な信号処理プログラムを拾ってきたら、それが浮動小数点数に基づくものだったというのは良くある話です。
嘘です。
良くある話なんて甘っちょろいものじゃありません。全部が全部浮動小数点数で書かれていると言っても過言ではない21世紀。この甘ったれどもめ、と愚痴り*1ながら固定小数点化するのは我ら組み込み族の風物詩とも言える風景です。
てなわけで、浮動小数点演算の固定小数点数化について書いて見ます。あまり難しいのもなんなので、FFTにします*2。
FFTと言えば、大浦さんのページですので出発点を大浦さんの周波数間引きFFTプログラムにします。リンク先から持ってきたのがこれ。
/* http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/ftmn1_2.html#sec1_2_1 */ #include <math.h> void fft(int n, double theta, double ar[], double ai[]) { int m, mh, i, j, k; double wr, wi, xr, xi; for (m = n; (mh = m >> 1) >= 1; m = mh) { for (i = 0; i < mh; i++) { wr = cos(theta * i); wi = sin(theta * i); for (j = i; j < n; j += m) { k = j + mh; xr = ar[j] - ar[k]; xi = ai[j] - ai[k]; ar[j] += ar[k]; ai[j] += ai[k]; ar[k] = wr * xr - wi * xi; ai[k] = wr * xi + wi * xr; } } theta *= 2; } /* ---- unscrambler ---- */ i = 0; for (j = 1; j < n - 1; j++) { for (k = n >> 1; k > (i ^= k); k >>= 1); if (j < i) { xr = ar[j]; xi = ai[j]; ar[j] = ar[i]; ai[j] = ai[i]; ar[i] = xr; ai[i] = xi; } } }
このプログラムは、一番外側のループが少しトリッキーでわかりにくいので読みやすく修正します。原作者の大浦さんは、この部分を特定のコンパイラで効率が良くなる様に書いたのだと思われます。しかしながら、移植は読みやすいコードから始めるのが鉄則です。そこでわかりやすくなるよう手をいれました。
元のコード：
for (m = n; (mh = m >> 1) >= 1; m = mh) { for (i = 0; i < mh; i++) {
手をいれた後：
for (m = n; m > 1; m = m/2 ) { /* n, n/2, n/4, ... , 2 */ mh = m/2; /* n/2, n/4, n/8, ..., 1 */ for (i = 0; i < mh; i++) {
変数mと変数mhの変化の範囲が明確になりました。割り算がありますが、2で割る場合はたいていコンパイラが処理してくれるので問題ありません。もちろん、自分のコンパイラがどうするかは知っておかなければなりません。VisualDSP++4.5のコンパイラも、bfin-elf-gcc 3.1.2も大丈夫でした。
次は明日。
