FFTを固定小数点化する(4)

前回の話を元に、大浦さんのFFTを書き換えてみましょう。下のコードは前半をとりだしたものです。

/* based on http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/ftmn1_2.html#sec1_2_1 */

void fft(int n, double theta, double ar[], double ai[])
{
    int m, mh, i, j, k;
    double wr, wi, xr, xi;

    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++) {
            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;
    }

このプログラムを前回の話を元に書き換えます。ポイントは

• doubleをshortに置き換える
• 乗算後15bit右シフトする
• 絶対値1未満の値は32768倍して格納する

の、3点です。結果は以下の通り。

/* based on http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/ftmn1_2.html#sec1_2_1 */

void fft(int n, double theta, short ar[], short ai[])
{
    int m, mh, i, j, k;
    short wr, wi, xr, xi;

    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++) {
            wr = cos(theta * i) * 32768;                 // ここ
            wi = sin(theta * i) * 32768;                 // ここ
            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 ) >> 15;    // ここ
                ai[k] = ( wr * xi + wi * xr ) >> 15;    // ここ
            }
        }
        theta *= 2;
    }

一番内側のループはバタフライ演算を実行する部分です。この部分には積が2行でつかわれていますので、その両者で15bitのシフトをおこなっています。
複素係数であるwr/wiはshortに変更していますが、この変数への係数生成はcos/sin関数を残しています。次回はこれをテーブル参照に変更しましょう。