FFTWについて

画像処理やってると FFT （高速フーリエ変換）をよく使います。名前に「高速」と付いてますが、それは DFT（離散フーリエ変換）を高速に実行出来るアルゴリズムだからです。 DFT を普通に実行するよりかなり高速に実行出来るので良いのですが、純粋な FFT アルゴリズムでは信号サイズが２の累乗（ 2, 4, 8, 16, ...）でないといけないと言う、かなり使いにくい仕様/制限があります。

そこで、この使いにくい仕様をどうにか簡単に克服しようと色々なアルゴリズムが提案されています。また、それらのアルゴリズムでは、簡単に FFT が出来るだけでは無く、 FFT を更に高速実行できるように工夫されています（普通に書かれた FFT より 5 倍程度早いです）。有名なアルゴリズムとしては FFTW, Apple vecLib, intel Math Kernel Library などがありますが、この中でも FFTW, Apple vecLib についてはフリーなので簡単に手に入れることが出来ます。 intel Math Kernel Library は intel の C コンパイラなどを買うとセットで付いてきます。これ高い。

インストール方法

まずは WEB サイトに行ってソースコードを取ってきます。バージョン2系列と3系列がありますが、自分のプログラムの中で FFT したいだけであればバージョン3系列のものを取ってきましょう。取ってきて展開したら、ここに書いてあるように

$ ./configure
$ make
$ sudo make install

で入ると思います。

簡単な使い方

#include <fftw3.h>
void fftw( TYPE reImg[line][pel],
	   TYPE imImg[line][pel],
	   int flag ){
  fftw_complex *in, *out;
  fftw_plan p;
  int i, j;

  in  = fftw_malloc(sizeof(fftw_complex) * line * pel);
  out = fftw_malloc(sizeof(fftw_complex) * line * pel);
  if( flag==0 ){
    p   = fftw_plan_dft_2d(line, pel, in, out, FFTW_FORWARD, FFTW_ESTIMATE); 
  }
  else {
    p   = fftw_plan_dft_2d(line, pel, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
  } 
  for( i=0; i<line; i++ ){
    for( j=0; j<pel; j++ ){
      in[j + pel * i][0]=reImg[i][j];
      in[j + pel * i][1]=imImg[i][j];
    }
  }
  fftw_execute(p);

  for( i=0; i<line; i++ ){
    for( j=0; j<pel; j++ ){
      reImg[i][j]=out[j + pel * i][0];
      imImg[i][j]=out[j + pel * i][1];
    }
  }
  fftw_free(in);
  fftw_free(out);
}

こんな感じに定義して

TYPE reImg[line][pel], imImg[line][pel];
...
...
...
fftw( reImg, imImg, 0 ); // forward
fftw( reImg, imImg, 1 ); // inverse

こんな感じで使う。

#include<stdio.h>
#include<math.h>
#define line 256
#define pel 256
#define COLOR 3
typedef double TYPE;
#include "fftw.h"
#include "out.h"

main( int argc, char *argv[] )
{
  char ch1[80], ch2[80];
  static TYPE reImg[line][pel], imImg[line][pel], tmp, outImg[line][pel][COLOR];
  static unsigned char inImg[line][pel][COLOR];
  FILE *fp1, *fp2;
  int i, j, k;

  fp1 = fopen( argv[1], "r" );
  if( fp1 == NULL ){
    printf("original file open error.?n");
  }
  fread( inImg, sizeof( unsigned char ), line*pel*COLOR, fp1 );
  fclose( fp1 );

  fp2 = fopen( argv[2], "w" );
  
  for( k=0; k<COLOR; k++ ){
    for( i=0; i<line; i++ ){
      for( j=0; j<pel; j++ ){
        reImg[i][j] = (double)inImg[i][j][k];
        imImg[i][j] = 0.0;
      }
    }
    fftw( reImg, imImg, 0 );

    for( i=0; i<line; i++ ){
      for( j=0; j<pel; j++ ){
        tmp = sqrt( reImg[i][j]*reImg[i][j]+imImg[i][j]*imImg[i][j] )/(100);
        if( (i<line/2) && (j<pel/2) ){
          outImg[i+(line/2)][j+(pel/2)][k]=tmp;
        }
        else if( (line/2<=i) && (pel/2<=j) ){
          outImg[i-(line/2)][j-(pel/2)][k]=tmp;
        }
        else if( (line/2<=i) && (j<pel/2) ){
          outImg[i-line/2][j+pel/2][k]=tmp;
        }
        else {
          outImg[i+line/2][j-pel/2][k]=tmp;
        }
      }
    }
  }
  x1out( outImg, fp2 );
}

コンパイル

$ cc -o fftw fftw.c -lfftw3 -lm -O2

実行テスト

$ ./fftw inputImg.raw outputImg.raw

これで、lena の周波数スペクトルが見れます。

ベンチマーク

単純に 512*512*3 の画像を FFT して IFFT するのを１セットにし、これを 100 回繰り返してみたところ

	普通のFFT	FFTW	vDSP
time[s]	88.71	35.13	54.16

となりました。繰り返し毎に fftw_malloc でメモリの動的確保と fftw_free での解放を行っているので、これを一回にすればさらに差が出ると思われます。