円周率.jp (http://xn--w6q13e505b.jp/method/fft/implement.html)

FFT を素直に実装すると計算量の導出に書いたような処理になる。 ただしこれらのコードは並べ替えのために単純にコピーを行ったり同じ $\omega$ の計算を繰り返し行ったりするなど効率がよくない。 そのため一般的にはそのような無駄を取り除いた方式の実装を FFT の実装として考えることが多い。 ここでは広く知られている 2 種類の実装方式を紹介する。

プログラムを実装する際にどちらの方が良いかは細かな実装や ユーザの需要次第となるので決定的な言及はできないが、 大まかには小メモリ動作を目指すなら In-place 型を、 速さを求めるなら Self-sorting 型を使用する印象である。

様々に提案されている FFT の各種高速化手法については基本的に 実装方法がどちらの方式であっても適用することができる。 In-place 型 FFT での効率化実装については[JW06]が詳しいので、 本サイトでは Self-sorting 型 FFT を用いて説明している。

冒頭部分で述べたような、再帰呼び出しや同じ $\omega$ の計算を素直に省いた実装を Self-sorting 型 FFT、もしくは開発者の名をとって Stockham 型 FFT という。 Self-sorting の名の通りデータが常に並んでいるため計算がしやすい。 特徴としてデータ整列のために元のデータと同程度のデータ領域を必要とする。 領域 A → 領域 B → 領域 A →… と計算結果を交互に置いていくので Ping-pong FFT とも言われる。

時間間引き FFT

void SelfSortDITCore(int n, int m, Complex* x, Complex* y) {
  const int l = n / m / 2;
  for (int k = 0; k < m; ++k) {
    double wt = -M_PI * k / m;
    Complex w = Complex::polar(0, wt);
    for (int j = 0; j < l; ++j) {
      int xi0 = j * m + k, xi1 = xi0 + n / 2;
      int yi0 = j * 2 * m + k, yi1 = yi0 + m;
      Complex wx1 = x[xi1] * w;
      y[yi1] = x[xi0] - wx0;
      y[yi0] = x[xi0] + wx0;
    }
  }
}

void SelfSortDIT(int n, Complex* x, Complex* y) {
  for (int m = 1, i = 0; m < n; m *= 2, ++i) {
    if (i % 2) {
      SelfSortDITCore(n, m, y, x);
    } else {
      SelfSortDITCore(n, m, x, (m == n / 2) ? x : y);
    }
  }
}

周波数間引き FFT

void SelfSortDIFCore(int n, int m, Complex* x, Complex* y) {
  const int l = n / m / 2;
  for (int k = 0; k < m; ++k) {
    double wt = -M_PI * k / m;
    Complex w = Complex::polar(0, wt);
    for (int j = 0; j < l; ++j) {
      int xi0 = j * 2 * m + k, xi1 = xi0 + m;
      int yi0 = j * m + k, yi1 = yi0 + n / 2;
      Complex t0 = x[xi0] + x[xi1];
      Complex t1 = (x[xi0] - x[xi1]) * w;
      y[yi0] = t0;
      y[yi1] = t1;
    }
  }
}

void SelfSortDIF(int n, Complex* x, Complex* y) {
  int m = n / 2;
  if (static_cast<int>(log2(n)) % 2) {
    SelfSortDIFCore(n, m, x, x);
    m /= 2;
  }
  for (int i = 0; m >= 1; m /= 2, ++i) {
    if (i % 2) {
      SelfSortDIFCore(n, m, y, x);
    } else {
      SelfSortDIFCore(n, m, x, y);
    }
  }
}

計算量の導出に書かれたコードを見直すと、 計算の前後 (DIT では前、DIF では後ろ) にデータの並べ替えのみを行っている 部分があることが分かる。 これは事前にデータの移動先、もしくは元のデータ位置が分かっていれば 並べ替えの処理、ならびに並べ替えのデータ領域が必要ない事を示していて、 実際に (2 基底の FFT では) ビット逆順 (bit reversal order) と呼ばれる順序に並び替えることで実現することができる。 この方式に基づいてメモリ領域の無駄を省いた実装を in-place 型 FFT、 もしくは Cooley-Tukey 型 FFT という。 一般にこの方式では計算前後にデータの並べ替え処理を行う必要があるが、 畳み込み乗算への応用のみを考える場合は 並べ替えない方針を取ることもできる。

以下に bit reverse の順序に並べ替えるコードを書いておく。

void BitReverse(int n, Complex* x) {
  for (int i = 1, j = 0; i < n; ++i) {
    for (int b = n >> 1; b; b >>= 1) {
      j ^= b;
      if (j & b) break;
    }
    if (i < j)
      swap(x[i], x[j]);
  }
}

時間間引き FFT

void InPlaceDIT(int n, Complex* x) {
  BitReverse(n, x);
  for (int m = 1; m < n; m *= 2) {
    for (int k = 0; k < m; ++k) {
      double wt = -M_PI * k / m;
      Complex w = Complex::polar(0, wt);
      for (int j = 0; j < n; j += 2 * m) {
        Complex &x0 = x[k + j], &x1 = x[k + j + m];
        Complex x1w = x1 * w;
        x1 = x0 - x1w;
        x0 = x0 + x1w;
      }
    }
  }
}

周波数間引き FFT

void InPlaceDIF(int n, Complex* x) {
  for (int m = n / 2; m >= 1; m /= 2) {
    for (int k = 0; k < m; ++k) {
      Complex w = Complex::polar(0, wt);
      for (int j = 0; j < n; j += 2 * m) {
        Complex &x0 = x[k + j], &x1 = x[k + j + m];
        Complex x1w = (x0 - x1) * w;
        x0 = x0 + x1;
        x1 = x1w;
      }
    }
  }
  BitReverse(n, x);
}