FFTのソフトウェアによる実装

[←]　[↑]　[→]

---

前々項で導いたCooley-TukeyのフローグラフからFFTのプログラムを書き下してみる。

まず、入力と各ステージのインデックスnがビットリバース順になっているので、これを自然な順番に振り直す。

各ステージの演算式は次のようになる。

このN=8の場合から一般のNの場合を類推すると、次の漸化式になる。

各ステージの変数は2次元配列x_m(N/2^m,2^m)になっているので、行数と列数がステージごとに変わるが要素数は常にNである。そこでプログラムで扱いやすいように

と1次元配列に置き換えると

さらにインデックス計算を簡単にするためにnを

に置き換えると、次のCooley-Tukey,DITの漸化式が得られる。ただし、元のnはビットリバーサルしていたことに注意。

結局、Cooley-Tukey,DITのFFTを計算するプログラムは次のとおりになる。

#include <complex>

static const int N = 1024; // FFT点数

static const double PI = 3.14159265358979; // 円周率

typedef complex<double> cpx_t; // 複素数タイプ

/*******************************************************************

 * サブルーチン - FFT Cooley-Tukey,DIT

 *******************************************************************/

void fft(cpx_t x[N],cpx_t X[N],cpx_t w[N],int bitrev[N]) {

  cpx_t xbuf;

/* 入力を並び替えてバッファに入れる */

  for(int n=0; n<N; n++) {

    X[bitrev[n]] = x[n];

  }

/* FFTの計算 */

  for(int a=1,b=N/2; a<N; a=a*2,b=b/2) {

    for(int k=0; k<a; k++) {

      for(int n=k; n<N; n=n+2*a) {

        xbuf = X[n+a] * w[b*k];

        X[n+a] = X[n] - xbuf;

        X[n]   = X[n] + xbuf;

      }

    }

  }

}

/*******************************************************************

 * メインルーチン

 *******************************************************************/

int main(int argc, char *argv[]) {

  static int bitrev[N];

  static cpx_t w[N];

  static cpx_t x[N];

  static cpx_t X[N];

/* ビットリバース順の事前計算 */

  bitrev[0] = 0;

  for(int a=1,b=N/2; a<N; a=a*2,b=b/2) {

    for(int k=0; k<a; k=k+1) {

      bitrev[k+a] = bitrev[k] + b;

    }

  }

 /* 係数の事前計算 */

  for(int k=0; k<N; k++) {

    w[k].real() =  cos(2*PI*k/N);

    w[k].imag() = -sin(2*PI*k/N);

  }

/* x(n)に値を設定 */

  for(int n=0; n<N; n++) {

    x[n] = 0;

    if(n==1) {x[n].real() = 1; x[n].imag() = 0;}

  }

/* x(n)をFFTしてX(k)を計算 */

fft(x,X,w,bitrev);

}

◇◇

同様にしてStockhamのアルゴリズムのFFTのプログラムも導いてみる。

このフローグラフは、前々項で導いた3ステージの式そのものであるから

このN=8の場合から一般のNの場合を類推すると、次の漸化式になる。

各ステージの変数は2次元配列x_m(N/2^m,2^m)になっているので、行数と列数がステージごとに変わるが要素数は常にNである。そこでプログラムで扱いやすいように

と1次元配列に置き換えると

結局、Stockham,DITのFFTを計算するプログラムは次のとおりになる。Cooley-Tukeyのようにビットリバース順への並べ替えは必要ないが、左辺と右辺の配列のインデックスが異なる(out-of-place)ので上書きできないことに注意。

/*******************************************************************

 * サブルーチン - FFT Stockham,DIT

 *******************************************************************/

void fft(cpx_t x[N],cpx_t X[N],cpx_t w[N],int bitrev[N]) {

  static cpx_t X1[N];

  cpx_t xbuf;

/* 入力をバッファに入れる */

  for(int n=0; n<N; n++) {

    X[n] = x[n];

  }

/* FFTの計算 */

  for(int a=1,b=N/2; a<N; a=a*2,b=b/2) {

    for(int n=0; n<b; n=n+1) {

      for(int k=0; k<a; k=k+1) {

        xbuf = x[a*n+N/2+k] * w[b*k];

        X1[2*a*n+k]   = X[a*n+k] + xbuf;

        X1[2*a*n+k+a] = X[a*n+k] - xbuf;

      }

    }

    for(int n=0; n<N; n=n+1) X[n]=X1[n];

  }

}

---

[←]　[↑]　[→]