FFTを固定小数点化する(6)
前回で演算と係数の固定小数点化、およびshort型でのエミュレートが完了しました。しかし、そのままではFFTに自然信号をかけても十中八九うまく行きません。多くの場合、値が途中で飽和してしまいます。
これを理解するには、FFTアルゴリズムの中核となるバタフライ演算について理解しなければなりません。バタフライ演算は次のような演算です。
A = (a+b)
B = (a-b)*W
ここでa,b,c,dは複素数、Wは複素回転子。
模式化すると、交差的な信号の動きが蝶の形に似ていることからバタフライ演算の名が付きました。
さて、バタフライ演算の入力a,bの値の範囲が次のようであると考えましょう
a ≦ 1.0 b ≦ 1.0
つまり、いずれの値も複素平面上の単位円の中に納まるとします。その場合、以下の公式*1
|a+b| ≦ |a|+|b|
から、
a+b ≦ 2.0 a-b ≦ 2.0
となります。つまり、バタフライ演算を通すと、値の絶対値は最大2倍になります。こんな演算を繰り返せばあふれるに決まっています。そこで、対策として各バタフライの直前で入力値を1/2にしてやります。演算に余裕がある場合ならばバタフライ演算の直後に値を1/2にしてもかまいません。
さっそくこの1ビットシフトをバタフライ演算コードに適用しましょう。元のコードは以下のようになっています。
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; }
これは最内側のループです。これに変更を施します。
for (j = i; j < n; j += m) { k = j + mh; xr = (ar[j]-ar[k]) >> 1; // ここ xi = (ai[j]-ai[k]) >> 1; // ここ ar[j] = (ar[j]+ar[k]) >> 1; // ここ ai[j] = (ai[j]+ai[k]) >> 1; // ここ ar[k] = ( wr * xr - wi * xi ) >> 15; ai[k] = ( wr * xi + wi * xr ) >> 15; }
この補正によって、入力が自然信号ならば確実にあふれを回避できます*2。こうしてえられた複素FFTルーチンの全貌を以下に示します。
/* 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 /= 2) { mh = m/2; 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]) >> 1; xi = (ai[j]-ai[k]) >> 1; ar[j] = (ar[j]+ar[k]) >> 1; ai[j] = (ai[j]+ai[k]) >> 1; ar[k] = (wr * xr - wi * xi) >> 15; ai[k] = (wr * xi + wi * xr) >> 15; } } 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; } } }
