So-net無料ブログ作成
  • ブログをはじめる
  • ログイン

fftwの使い方 [プログラミング]

自分でインストールしておきながら忘れてしまうのでメモしておこう。

FFTのフリーのライブラリとしてfftwというのがある。
http://www.fftw.org/
マルチプラットフォームで使い勝手が良くて最高のパフォーマンスでかつフリーという売りで有名。
MacにはnativeのvDSPと言うベクトルユニットを使うライブラリがあるが、ベクトルユニットがサポートしていないデータには当然使えない。
fftwはベクトルユニットを指定するコンパイル用スイッチがある。

ダウンロードはhttp://www.fftw.org/download.htmlから。07/10現在3.1.2がstable。3.2 alpha2も人柱向けに上がっている。
MacではFinkを使うことも出来る。

まず、ファイルを展開

% tar xvf fftw-3.1.2.tar


したあとのフォルダ

% cd fftw-3.1.2


に入り、いつもの

% ./configure
% make
% sudo make install


で/usr/localにインストールされる。
configureにはいっぱいスイッチがあり、いつものディレクトリ変更のスイッチの他に
1.計算精度を変える(デフォルトはdouble)
-enable-float
--enable-long-double(これはコンパイラ依存)
2.複数のthreadを使う
--enable-threads
3.ベクトルユニットを使う
--enable-sse, --enable-sse2, --enable-k7, --enable-altivec
等がある。
スイッチの簡単な説明は

% ./configure --help


でリストされる。
精度を変えたら別々のファイルができる。

% ls /usr/local/lib/libfft*
/usr/local/lib/libfftw3.a       /usr/local/lib/libfftw3f.a
/usr/local/lib/libfftw3.la*     /usr/local/lib/libfftw3f.la*


fftw3.*がdouble用でfftw3f.*がfloat用になっているので、リンクする時注意する。
ヘッダを読んでも良くわからないけどreferenceには関数の名前は
fftw_foo double用
fftwf_foo float用
fftwl_foo long double用
であるとなっている。
それ以外はファイルとしては同じ(上書き)。
例えばPowerBookでG4のaltivecベクトルユニットを使うfloat幅のライブラリを生成するには

% ./configure  --enable-float --enable-altivec


とすればよい。

使い方は
1.実データ領域と、Fourier領域のメモリをallocateする
2.planと呼ばれる構造体を作る。これにはサイズや正逆変換、アルゴリズムの指定などを設定する
3.fftを実行する
4.planを破棄する
となる。

特徴的なのはplanで設定するPlanning-rigor flagsである。内部アルゴリズムがいくつかあって速いのを選べるらしい。
flagには

FFTW_ESTIMATE
FFTW_MEASURE
FFTW_PATIENT
FFTW_EXHAUSTIVE


の四つがある。
適当なのを選ばせるFFTW_ESTIMATE、実際にちょっと計算してその結果で選択するFFTW_MEASURE、もうちょっとやってみてさらにましなのを探すFFTW_PATIENT、さらにがんばって一番ましなのを探すFFTW_EXHAUSTIVEとなっているとreferenceに書いてある。
入力配列を上書きしても良いというflagを設定してFFTW_ESTIMATE以外でplanを作ると、planが帰ってきた時、データ領域の一部が書き換えられている。試しにやってみたと言うことらしい。これ、知らないでデータを用意してからplanを作ってfftを実行するとでたらめな結果が出来てしまうので長く悩んでしまった。
planの中身は

void fftw_print_plan(const fftw_plan plan);


などの関数で書き出せるらしいが「nerd-readable」と書いてあってナードでなければ読めないようだ。

具体的にはチュートリアル
http://www.fftw.org/fftw3_doc/Complex-One_002dDimensional-DFTs.html#Complex-One_002dDimensional-DFTs
に一番一般的な1次元複素FFTのやり方が書いてある。

複素数のデータは

typedef double fftw_complex[2];


と言う型の配列、

fftw_complex	*in;


等とする。
in[0]に実部、in[1]に虚部が入るという簡単な構造になっている。2次元以上の配列はCで普通の並び方にすればいい(FORTRANのように縦が先、ということはない)ので、FFT以外の部分との相性も別段問題ない。2次元以上の場合packされていてもアラインメント等の都合で空いていてもOK。

実際にXcodeのprojectで使うには
1.「プロジェクト」メニューの「プロジェクトに追加...」を選ぶ
2.OpenPanelからインストールしたlibfftw3.aを指定する
3.「プロジェクト設定」を開いてヘッダ検索パスとライブラリ検索パスにfftw3.hとlibfftw3.aにパスが通るようにする。デフォルトでは/usr/local/lincludeと/usr/local/libになる
4.fftw3.hをソースファイルにインククードする
5.コンパイルする
でOK。

もし自分のマシン以外でも動く*.appを作りたかったら、アプリバンドルの中にlibfftw3.aを含めた方がいい。やったこと無いので良くわからないけど。

実際のコードは例えばコマンドラインで

#include < stdio.h>
#include < fftw3.h>
#include < math.h>

main()
{
	static const int	fftSize = 64;
	fftw_complex	*in, *out;
	fftw_plan	plan;
	int		i;

	in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fftSize);
	out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fftSize);
	plan = fftw_plan_dft_1d(fftSize, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

	fftw_print_plan(plan);
	printf("\n");

	for (i = 0 ; i < fftSize ; i ++) {
		in[i][0] = cos(M_PI *16.0 * i / fftSize) / fftSize;
		in[i][1] = 0.0;
	}

 	fftw_execute(plan); 
	fftw_destroy_plan(plan);
	fftw_free(in);
	for (i = 0 ; i < fftSize ; i ++)
		printf("[%4d]:(%lf,\t%lf)\n", i, out[i][0], out[i][1]);
	fftw_free(out);
}
で、
% cc fftwTest.c -lfftw3 -lm -o fftwTest
% ./fftwTest
とすれば、この場合、8番目と56番目(64-8番目)の実数部以外は0になっていると成功。 要素数はFFTなので当然2^nの場合が最も効率がいいが、任意の要素数で可能である。例えばソースのfftSizeを61等の素数にしてみてもちゃんと動く。その場合、fftw_print_plan()の出力が変わるがナードでないのでよくわからない。
nice!(0)  コメント(0)  トラックバック(0) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 0