余白の書きなぐり

aueweのブログ

高速フーリエ変換ライブラリ fftw のインストールと使用法

インストール

Linux(Ubuntu)への導入

インストール

$ sudo apt-get install libfftw3-3 libfftw3-dev libfftw3-doc

コンパイル方法

$ gcc hoge.c -lm -lfftw3
Windowsへの導入

fftwのwindows版インストーラのページから
32-bit version: fftw-3.3.3-dll32.zipをダウンロードした。
64bit環境なら64-bit version: fftw-3.3.3-dll64.zipの方が良いかも。

zipファイルを解凍すると、20種類ぐらいのファイルが生成される。その中に

  • fftw3.h (C:\path\to\fftw3.hに置いたとする)
  • libfftw3-3.dll (C:\path\to\libfftw3-3.dllに置いたとする)

の2つがある。この2つのファイルを適当な場所に置く。 コンパイル方法は cmd.exe を立ちあげて以下のコマンドを実行

gcc hoge.c -lm C:/path/to/libfftw3-3.dll

C:/path/to/libfftw3-3.dllの部分は
./相対パス/libfftw3-3.dllでもOK

サンプルコード

 (a_0, a_1, a_2, a_3) = (1,2,5,3)
という4成分からなる数列 {a} があったとしよう。 こいつに、N=4として
{\displaystyle b_n = \sum_{m=0}^{N-1} a_m \exp \left( -2 \pi i \frac{m}{N}n \right)}
というフーリエ変換をかませば、以下の数列{b}が得られる。
 b_0 = 1 (+1) + 2 (+1) + 5 (+1) + 3 (+1) =  11
 b_1 = 1 (+1) + 2 (-i) + 5 (-1) + 3 (+i) = -4  + i
 b_2 = 1 (+1) + 2 (-1) + 5 (+1) + 3 (-1) =  1
 b_3 = 1 (+1) + 2 (+i) + 5 (-1) + 3 (-i) = -4  - i

以上を計算するプログラムをサンプルとして用意した。

/***************************/
/* FourierTransformation.c */
/***************************/

// 要素数Nの数列 a[m] をフーリエ変換で b[n] にする。添字は 0,...,N-1
// b[n] = sum[m=0,N-1]  a[m]*exp(-2*PI*I* (m/N) *n)

#include <stdio.h>
#include <complex.h> // complex.h は fftw3.h より先に include する
#include <fftw3.h>   // windows環境では #include "C:/path/to/fftw3.h"
                     // あるいは        #include "./相対パス/fftw3.h"
 
int main( void ){
 
  int N=4;

  // a,b は double _Complex 型のC99標準複素配列と実質的に同じ
  // double _Complex a[4] としても動くけど計算速度が低下する可能性あり
  fftw_complex *a, *b;
  a = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  b = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

  // プランの生成
  // フーリエ逆変換つまり位相因子を exp(-k)ではなく exp(+k)とする場合は
  // FFTW_FORWARD の代わりに FFTW_BACKWARD とする
  fftw_plan plan;
  plan = fftw_plan_dft_1d( N, a, b, FFTW_FORWARD, FFTW_ESTIMATE );
 
  // フーリエ変換前の数列値を設定
  a[0] = 1.0 + 0.0*I;
  a[1] = 2.0 + 0.0*I;
  a[2] = 5.0 + 0.0*I;
  a[3] = 3.0 + 0.0*I;
 
  // フーリエ変換実行   b[n]に計算結果が入る
  fftw_execute(plan);
 
  // b[n]の値を表示
  int n;
  for( n=0; n<N; n++ ){
    printf("b_%d = %+lf %+lf*i\n", n, creal(b[n]), cimag(b[n]) );
  }

  // ここで a[m] の値を変えて再度 fftw_execute(plan) を実行すれば、
  // b[n] が再計算される。

  // 計算終了時、メモリ開放を忘れないように
  if(plan) fftw_destroy_plan(plan);
  fftw_free(a); fftw_free(b);
 
  return 0;
}

コンパイルと実行

$ gcc main.c -lm -lfftw3
$ ./a.out
b[0] = +11.000000 +0.000000*i
b[1] = -4.000000 +1.000000*i
b[2] = +1.000000 +0.000000*i
b[3] = -4.000000 -1.000000*i

参考:複素1次元離散フーリエ変換 Complex One-Dimensional DFTs-FFTW@wiki