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

FFT

前回、残しておいたcos/sin関数を固定小数点化しましょう。FFTに限った話をするならば、これはテーブル参照にするのが一般的です。おそらくはFFTがベンチマークに使われた影響が大きいのですが、関数近似にして雑音の影響が大きくなるのはたまらんといった点もあるかと思います。
さて、テーブル参照にする場合、問題になるのはテーブルの大きさです。N点FFTの場合、cos/sin双方にN/2エントリのテーブルが必要になります*1。N点より小さなFFTにかんしては、N点FFTとテーブルを共用できるものの、設計時にNの最大値をいくつにとるか考えておかなければなりません。
あらかじめNの最大値を想定するのが不合理なら、FFT関数に引数としてテーブル・アドレスをあたえることになります。ここでは、

  • テーブルを外部から引数として与え
  • 小さなNのFFTとテーブルを共用できる

ようにしましょう。
この方針で作るFFTは、まずは関数ヘッダ部が変わります。以下、オリジナルとは前回の関数のことです。

/*オリジナルのヘッダ*/
void fft(int n, double theta, short ar[], short ai[])

/*変更後*/
void fft(int n, int theta, short ar[], short ai[], short cos[], short sin[] )

short型配列cos/sinを引数としてとる他、thetaが整数型に変更されています。
thetaは、Nの値を可変にしなければ無用な変数です。仮りにN点FFT用のテーブルを使って、M点のFFTを実行するとします。N>Mです。そのばあい、thetaにはN/Mの値を与えます。これによって不必要に細かいFFTテーブルを利用するときにも正しい答えを得ることができます。
cos/sinをshort型配列にしたことで、doubleからshortへの変換時のスケーリングも不要になります。

/* 変更前 */
            wr = cos(theta * i) * 32768;
            wi = sin(theta * i) * 32768;

/* 変更後 */
            wr = cos[theta * i];
            wi = sin[theta * i];

以上の変更を施したコードをいかに示します。

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

                // n : データ点数
                  // theta : m*2/n
                // ar[] : 入力データの実部(n点)
                // ai[] : 入力データ虚部(n点)
                // cos[] : 余弦テーブル(m点。cos(0..π))
                // sin[] : 正弦テーブル(m点。sin(0..π))
void fft(int n, int theta, short ar[], short ai[], short cos[], short sin[] ) // 変更
{
    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];                  // 変更
            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 ) >> 15;
                ai[k] = ( wr * xi + wi * xr ) >> 15;
            }
        }
        theta *= 2;
    }

さて、外部から与えるcos/sinテーブルですが、いくつか注意が必要です。
まずサイズはN/2エントリとし、引数が0からπまでに対応する値を格納します。つまり、半周期のテーブルとします。
また、値の範囲を[-1...+1]から32768倍した値にしなければなりません。ところが、+32768は16bit符号つき整数では表現できませんので、これは32767に修正しておきます。
以上でFFTの固定小数点化はおわりです。好きな波形を放りこんで楽しみましょう。と、いいたいところですがそうは問屋が卸しません。この関数に例に適当な波形を放りこんでも十中八九おかしな出力を生成します。あふれが生じるからです。
あふれを抑えるには、FFTの中核であるバタフライ演算について考えなければなりません。

*1:N/2エントリの複素テーブル