2012年4月12日木曜日

FFTW


通信分野において出番の多い離散フーリエ変換。ディジタル処理ではそれを高速に計算するアルゴリズム(高速フーリエ変換, FFT)が使われ、みなさんもシミュレーションに組み込む機会があるかと思います。
1度自分で作りその仕組を理解することは非常に重要なので最初から使うことは勧めませんが、その後研究を進めていく上で計算速度やバグ取りの効率など考えると配布されているものを使ってもいいと思います。


・FFTは「FFTW」というものが配布されていますので、これからその導入の仕方を追っていきましょう。


【セットアップ】
1.まず、FFTWをパソコンにインストールする必要があります。私が導入した時点では「fftw-3.2.2」くらいのヴァージョンです。google先生に「FFTW」「mac」とか入れれば教えてくれるでしょう。


2.インストールが完了したら、ターミナルでホームフォルダから/usr/local/includeの場所に「fftw3.f」「fftw3.h」があるか確認しましょう。ホームフォルダから
ls /usr/local/include
とコマンドすれば見れます。確認できたら、前準備はおしまいです。




ここからは使い方を具体的なソースコードを例にしながら確認していきます。
実際あんまよくわかってないんですがわかっている事を書いていくので、参考程度に。

#include <complex.h>
#include <fftw3.h>
#define FFTSIZE 4
// Complex型を使うと便利です。fftw3.hをインクルードします。
// FFTSIZEはNとかなんでもいいです。FFTの大きさです。

int main(void)
{
int i;
complex *in, *out1, *out2;
fftw_plan p1,p2;
// fftw_complexという型はcomplexと思っておkです。
// fftw_planというのは見慣れないのですが、これにどういうFFTを行うのかの入れます。

in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * FFTSIZE);
out1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * FFTSIZE);
out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * FFTSIZE);
// complex型の配列を準備してるだけです。

in[0] = 1 + 1*I;
in[1] = 1 + 1*I;
in[2] = 0 ;
in[3] = 0 ;
// 適当に数値を設定しました。好きなように変えてください。

p1 = fftw_plan_dft_1d(FFTSIZE,in,out1,FFTW_FORWARD,FFTW_ESTIMATE);
p2 = fftw_plan_dft_1d(FFTSIZE,out1,out2,FFTW_BACKWARD,FFTW_ESTIMATE);
// 第一引数にFFTの大きさ、第二引数に入力、第三引数に出力、
// 第四引数はFFTなのかIFFTなのかどうか、第五引数は要調査。

fftw_execute(p1);
 // これがFFTを実行する関数
fftw_execute(p2);
for(i=0;i<FFTSIZE;i++){
out2[i] /= FFTSIZE;
}
 // FFTWは1/Nしてくれないので自分で正規化しましょう
  // 通信の場合はFFTとIFFTの両側で1/sqrt(N)とかですね。
  // sqrtを使う場合はmath.hのインクルードを忘れずに。


fftw_destroy_plan(p1);
fftw_destroy_plan(p2);
// planを破棄、まーいちお書きましょう。

for(i=0;i<FFTSIZE;i++){
printf("%f %f\n",creal(out1[i]),cimag(out1[i]));
}
putchar('\n');
for(i=0;i<FFTSIZE;i++){
printf("%f %f\n",creal(out2[i]),cimag(out2[i]));
}
// 結果表示

fftw_free(in); fftw_free(out1); fftw_free(out2);
return 0;
}

それでは、また(^o^)ノ

0 件のコメント :

コメントを投稿