ページ

2011年4月16日土曜日

大浦氏のFFTの使い方

さくらの花が、開いてきましたよん。

これまでC/C++でのFFTには、FFTWもしくは自前のコードを使っていた。FFTWは非常に高速なFFTライブラリとして有名なのだけど、ライブラリをコンパイルしてインストールした計算機でないと使えないので、人にプログラムを渡すときに説明するのが面倒。自前のFFTはあまり速くない(SlowなFast Fourier Transformation)。というわけで、FFTWと同等かそれ以上に高速といわれる大浦氏のFFTを使ってみた。ありがたいことに、こののライブラリは商用/非商用を問わず自由に利用して良いと記されている。FFTWと比べると、関数の使い方にやや癖があるのでメモ。

【 使用するファイル 】
配布されている関数群は、基数や作業用配列使用の有無ごとにいくつかのファイルに分けられている。今回は作業用配列を用いたSplit-radixの"fftsg.c"を使うことにする。組み込み系などを除けば、たぶん多くの環境でこれが一番高速だろう。

【 複素DFTの使い方 】
void cdft(int n, int isgn, double *a, int *ip, double *w);
nFFT点数×2の値を入力する。
sign-1または1。一般的な定義の場合、-1:離散フーリエ変換 1:離散逆フーリエ変換。
*a入出力配列であり。添字が偶数番目の位置に実部、奇数番目の位置に虚部を置く(だから配列の長さがFFT点数の2倍になる)。FFTの結果は配列aにそのまま上書き出力される(in-place演算)。
*ip作業用配列。長さは2+sqrt(n)以上が必要で、最初にcdft()を呼ぶときにはip[0]=0を代入しておく。
*w三角関数表用配列を渡す。長さはn/2。ip[0]=0が渡された場合に初期化される。
作業用配列ipの値には少し注意が必要で、プログラム内で最初にcdft()を呼ぶときと、FFT点数を変更する場合には、ip[0]=0を予め代入してから渡す。これは作業用配列ipと三角関数表wを初期化するためのフラグになっている。

【 使用例 】
サンプルはシンプルに。440.0[Hz]の正弦波信号をFFTして、IFFTでもとに戻す。
元信号と、FFT+IFFTで戻した信号。当然ながらピタリ一致。

パワースペクトル(一部)。

ソースコードはこちら。いろいろなプログラムにincludeして使いたいので、ヘッダファイル"fftsg.h"を作った。

ヘッダファイル"fftsg.h"

使用例

1 件のコメント:

やもり さんのコメント...

Mac用のFFTソフトに大浦さんのコードを組み込み際に参考にさせてもらいました。ありがとうございました。

http://hp.vector.co.jp/authors/VA009278/sparrow/index.html