FFTW

FFTW†

Visual C++ でのインストール†

1. http://www.fftw.org/install/windows.html から fftw-3.2.2.pl1-dll32.zip をダウンロードし， 適当な場所に解凍する．
2. コマンドプロンプトでVC++の環境変数をセットする（VC++のbinフォルダのvcvars32.batをコマンドラインで実行するか， スタートメニューから Visual Studio 2005 コマンド プロンプト などを起動）． 1で解凍したフォルダに移動し，

   lib /def:libfftw3-3.def
   lib /def:libfftw3f-3.def
   lib /def:libfftw3l-3.def

   を実行する．
3. プロジェクトのインクルードフォルダに

   fftw3.h

   ライブラリフォルダに

   libfftw3-3.lib
   libfftw3f-3.lib
   libfftw3l-3.lib

   実行フォルダに

   libfftw3-3.dll
   libfftw3f-3.dll
   libfftw3l-3.dll

   を入れておく． それぞれ，

• libfftw3-3 ：単精度
• libfftw3f-3 ：倍精度
• libfftw3l-3 ：ロング精度 である．

実数データの多次元DFT†

実数データの多次元DFTには以下のプランナを用いる．

fftw_plan fftw_plan_dft_r2c_2d(int n0, int n1, double *in, fftw_complex *out, unsigned flags);
fftw_plan fftw_plan_dft_r2c_3d(int n0, int n1, int n2, double *in, fftw_complex *out, unsigned flags);
fftw_plan fftw_plan_dft_r2c(int rank, const int *n, double *in, fftw_complex *out, unsigned flags);

r2cはreal to complexのことであり，complex to real の場合は関数名のr2cをc2rに置き換えればよい．
2次元の場合，入力inにはn0×n1の2次元データがrow-major orderで格納された配列を， 出力outにはn0×(n1/2+1)の大きさのfftw_complex型の配列のアドレスを渡す． プランナに渡すflagについては，Planner Flagsを参照．

出力のサイズが異なるのは，DFT変換後のデータがエルミート冗長性を持ち， out[i]がout[n-i]の共役になっており，メモリの節約のためにサイズを小さくしている． 描画などでフルデータが必要ならば

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14

    for(int i = 0; i < w; ++i){
        for(int j = 0; j < h; ++j){
            double re, im;
            if(j >= h/2+1){
                int j0 = h-j-1;
                re =  out[j0+i*h1][0];
                im = -out[j0+i*h1][1];
            }
            else{
                re = out[j+i*h1][0];
                im = out[j+i*h1][1];
            }
        }
    }

のように変換する．

2次元の画像データからの変換例として，

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31

 void FFTW_c2r(Array2 &reImg, Array2 &imImg, int w, int h)
 {
    fftw_complex *in;
    double *out;
    fftw_plan p;
    in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*w*(h/2+1));
    out = (double*)fftw_malloc(sizeof(double)*w*h);
    
        p = fftw_plan_dft_c2r_2d(w, h, in, out, FFTW_ESTIMATE);
    // 
    for(int i = 0; i < w; ++i){
        for(int j = 0; j < h/2+1; ++j){
            in[j+(h/2+1)*i][0] = reImg[i][j];
            in[j+(h/2+1)*i][1] = imImg[i][j];
        }
    }
    
    fftw_execute(p);
    
    for(int i = 0; i < w; ++i){
        for(int j = 0; j < h; ++j){
            reImg[i][j] = out[j+h*i];
            imImg[i][j] = 0.0;
        }
    }
    
    fftw_free(in);
    fftw_free(out);
    fftw_destroy_plan(p);
 }

2D FFTで得られたデータは四隅が低い周波数，中心が高い周波数の値を示している． 一般的な中心が低い周波数の状態にしたい場合は，第一象限と第三象限，第二象限と第四象限を入れ替えればよい．

TIPS†

• 逆フーリエ変換したデータはスケーリングされていないので，データ数で割ること

  reImg[i][j]/(w*h);

• Planner Flagの FFTW_PRESERVE_INPUT は入力データの保存を保証する(c2r,hc2r以外ではデフォルト)が， 多次元c2r変換でについてはこれは実装されておらず，これが指定されるとプランナはNULLを返す．