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

[] [] []


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

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

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

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

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

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

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

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

結局Cooley-Tukey,DITFFTを計算するプログラムは次のとおりになる。ただし複素数の四則演算が見やすいように疑似コードで書いてある。

/* FFTの疑似コード - Cooley-Tukey,DIT */

complex w[N],x[N],xbuf;

integer bitrev[N],n,k,a,b;

/* 係数の事前計算 */

for(n=0; n<N/2; n=n+1) w[n] = exp(-2πj*n/N);

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

bitrev[0] = 0;

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

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

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

  }

}

/* 入力の並び替え */

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

  if(n < bitrev[n]) {

    xbuf = x[n];

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

    x[bitrev[n]] = xbuf;

  }

}

/* FFTの計算 */

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

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

    for(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;

    }

  }

}

◇◇

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

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

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

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

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

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

/* FFTの疑似コード - Stockham,DIT */

complex w[N],x[N],x1[N],xbuf;

integer n,k,a,b;

/* 係数の事前計算 */

for(n=0; n<N; n=n+1) w[n] = exp(-2πj*n/N);

/* FFTの計算 */

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

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

    for(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(n=0; n<N; n=n+1) x[n]=x1[n];

}


[] [] []