一般に，バタフライ演算でデータの上書き(in-place 演算)を行う ためには，データの並べ替え(スクランブル)が必要になります．たとえば， 周波数間引きの FFT では，出力データの並べ替えまたは入力データと 三角関数表の並べ替えが必要になります．この並べ替えは，基数が 2 の 整数乗のときに，ビット反転順の並べ替えになります．ここでは主に， 少ない作業領域で行えるビット反転の並べ替えの方法について考えます．
一般にテーブルを用いないビット反転の並べ替えの方法は，ビット反転の カウンタを用いる方法をとります．すなわち，通常のカウンタ j に 1 を 加えたとき，それに対応するビット反転カウンタ i に上位のビットが 一ビット目とみなして 1 を加えることをします．なぜこのような面倒な ことをするかというと，一般の j に対してそのビット反転 i を直接計算する 方法は，たとえアセンブラを用いたとしても重い計算となることが知られて いるからです．次に，この方法によるビット反転並べ替えのルーチンを 示します．
| リスト1.2.5-1. ビット反転並べ替え |
|---|
/* ビット反転並べ替えを行う */
void bitrv(int n, double a[])
{
int i, j, k;
double tmp;
i = 0;
for (j = 1; j < n - 1; j++) {
for (k = n >> 1; k > (i ^= k); k >>= 1);
if (j < i) {
tmp = a[j];
a[j] = a[i];
a[i] = tmp;
}
}
}
|
このルーチンでは，ビット反転の添字の対 (i,j) の交換は (j,i) の交換と 重複するので， j < i となる対についてのみ交換しています． このルーチンは，排他的論理和を用いてビット反転のカウントを行って いますが，ビット演算を避けたい場合は，排他的論理和を加減算に 置き換えます．ビット演算を四則演算に置き換えると次のようになります．
| リスト1.2.5-2. ビット演算を用いないビット反転並べ替え |
|---|
/* ビット反転並べ替えを行う */
void bitrv(int n, double a[])
{
int i, j, k;
double tmp;
i = 0;
for (j = 1; j < n - 1; j++) {
k = n / 2;
while (k <= i) {
i = i - k;
k = k / 2;
}
i = i + k;
if (j < i) {
tmp = a[j];
a[j] = a[i];
a[i] = tmp;
}
}
}
|
次に，対称性を用いてビット反転の計算を減らすことを考えます． 例えば，k のビット反転を rev(k) とすると
2j のビット反転は rev(2j)
2j + 1 のビット反転は rev(2j) + N/2
2j + N/2 のビット反転は rev(2j) + 1
2j + N/2+1 のビット反転は rev(2j) + N/2+1
となるので，0 <= j < N/4 の範囲ですべてのビット反転の対ができる ことになります．この対称性を利用すれば，ビット反転のカウントの計算を 1/4 に減らすことができます[参考文献]． 次に，このルーチンを示します．
| リスト1.2.5-3. より効率的なビット反転並べ替え |
|---|
/* ビット反転並べ替えを行う */
void bitrv(int n, double a[])
{
int nh, nh1, i, j, k;
double tmp;
nh = n >> 1;
nh1 = nh + 1;
i = 0;
for (j = 0; j < nh; j += 2) {
if (j < i) {
tmp = a[j];
a[j] = a[i];
a[i] = tmp;
tmp = a[j + nh1];
a[j + nh1] = a[i + nh1];
a[i + nh1] = tmp;
}
tmp = a[j + nh];
a[j + nh] = a[i + 1];
a[i + 1] = tmp;
for (k = nh >> 1; k > (i ^= k); k >>= 1);
}
}
|
次に，ほんの少しの作業領域を用いて効率的なビット反転並べ替えを行う 方法を示します．この方法はよく知られた方法で，次のように計算します [参考文献]． まず，k のビット反転を rev(k) とすると
M = sqrt(N) が整数の場合:
j + rev(i) のビット反転は rev(j) + i
M = sqrt(N/2) が整数の場合:
j + rev(i) のビット反転は rev(j) + i
j + M + rev(i) のビット反転は rev(j) + M + i
となるので，0 <= j < M, 0 <= i < M の範囲ですべての ビット反転の対ができることを利用します．ビット反転対の作成は， ビット反転列を格納するための M の長さの配列を用いて行います． ビット反転対の交換は，重複を避けるため j < i となる対についてのみ 交換します．このとき，ループ内の条件分岐を排除することができます． このときのルーチンは次のようになります．
| リスト1.2.5-4. 作業領域を用いたビット反転並べ替え |
|---|
/* ビット反転並べ替えを行う．
ip[] は作業領域で，サイズは sqrt(n) または sqrt(n/2) の
割り切れる方．
*/
void bitrv(int n, double a[], int ip[])
{
int i, ij, j, ji, k, m;
double tmp;
/* ---- initialization ---- */
ip[0] = 0;
k = n;
m = 1;
while (2 * m < k) {
k = k / 2;
for (j = 0; j < m; j++) {
ip[m + j] = ip[j] + k;
}
m = m * 2;
}
/* ---- scramble ---- */
if (m == k) {
for (i = 1; i < m; i++) {
for (j = 0; j < i; j++) {
ji = j + ip[i];
ij = i + ip[j];
tmp = a[ji];
a[ji] = a[ij];
a[ij] = tmp;
}
}
} else {
for (i = 1; i < m; i++) {
for (j = 0; j < i; j++) {
ji = j + ip[i];
ij = i + ip[j];
tmp = a[ji];
a[ji] = a[ij];
a[ij] = tmp;
tmp = a[ji + m];
a[ji + m] = a[ij + m];
a[ij + m] = tmp;
}
}
}
}
|
このルーチンの作業領域は例えば N = 2048 のとき M = 32 で済みます． このルーチンは sqrt(N) が割り切れる場合とそうでない場合とに場合分けを していますが，さらに少しの工夫で一つのループにまとめることもできます． しかし，これは速度的にはあまり意味はないと思います．
最後に，in-place 型の混合基数 FFT のデータの並べ替えについて少し 考えます．この場合の並べ替えはビット反転の並べ替えにはならず， その場での並べ替えが困難になります．この並べ替えの方法は，まず添字の 順序を格納する整数型配列を用意して，並べ替えの順序を計算しておきます． 次に，このテーブルを参照して並べ替えを行います．この混合基数 FFT の データの並べ替えはビット反転の並べ替えよりも複雑になり，二つのデータの 交換では並べ替えができなくなります．このため，並べ替えたデータを 格納するための実数型配列が必要になります．