「Complex One-Dimensional DFTs」の編集履歴（バックアップ）一覧はこちら
追加された行は緑色になります。
削除された行は赤色になります。
* 複素１次元離散フーリエ変換 Complex One-Dimensional DFTs
プランとは、偶然の結果を達成するための最善の方法をいろいろ思い悩むことである。
(Ambrose Bierce, The Enlarged Devil’s Dictionary.)
----
FFTWで、サイズNのデータの一次元DFTを計算する方法はいたって簡単で、おおよそ以下のようなコードになります。
#include <fftw3.h>
...
{
fftw_complex *in, *out;
fftw_plan p;
...
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
...
fftw_execute(p); /* 必要なだけ繰り返す */
...
fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);
}
コンパイル時には、fftw3ライブラリをリンクしなければなりません。つまり例えばUnix系では -lfftw3 -lmをリンカオプションとしてつけます。
まず初めに、入出力の配列のメモリを確保します。
メモリは好きなように確保して構いませんが、fftw_mallocを使うことを推奨します。
fftw_mallocはmallocのように動作しますが、違いとしてはSSEやAltivecなどのSIMD命令が利用出きるように配列が確保されます。(3.1.1節 SIMD alignment and fftw mallocを参照)
データはfftw_complex型の変数の配列になっており、fftw_complexはデフォルトでdouble二つ分の配列(double[2])です。それぞれ実部(in[i][0])、虚部(in[i][1])に対応します。
次のステップは、プランの生成です。プランはFFTWがFFTを計算する際に必要な全ての情報を含んだものです。以下の関数がプランを作ります。
fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
最初の引数nは計算しようとしている変換のサイズです。nは任意の正の整数をとりえますが、サイズが小さい素数のみを因数に持つ場合、最も効率よく計算できます。(しかしどんな素数であれ,O(nlogn)アルゴリズムを採用しています)
計算時間の次元が、データ数nに対して、nlogn
次の二つの引数は、入出力の配列へのポインタです。これらのポインタは同じものを入れることも可能で、この場合は入力配列が書き換わります。
四番目の引数signは、正変換か逆変換かを表す、FFTW_FORWARD(-1)かFFTW_BACKWARD(+1)のどちらかを入れます。技術的にはこれは、変換の指数部になります。
最後の引数flagは、通常FFTW_MEASUREかFFTW_ESTIMATEをとります。FFTW_MEASUREは実際にいくつかのFFTWを実行して実行時間を計り、サイズnの変換を計算するのに一番良い方法を選択します。このプロセスは少し時間がかかり、(通常、2、3秒)これは変換のサイズnや計算機環境に依存します。FFTW_ESTIMATEは、実際に実行することなく、おそらく準最適だと思われる方法を選択します。 つまり、同じサイズの変換を複数回行い、かつ初期化の時間が重要でない場合はFFTW_MEASUREを使い、そうでなければFFTW_ESTIMATEを利用すれば良いことになります。
inやoutといった入出力データの配列は、FFTW_MEASUREのプランニングを行う際に上書きされるので、FFTW_MEWSUREでのプランニングを行う場合は、入力がユーザーに初期化される前にプランニングを行う必要があります。
一旦プランが生成されれば、特定のint/out配列に対して実際の変換をfftw_execute(plan)で計算する際に、プランは何度でも再利用できます。
ちなみにfftw_executeは,
void fftw_execute(const fftw_plan plan);
です。
もし、同じサイズの違う配列を変換したい場合は、fftw_plan_dft_1dを使って新しいプランを作れます。するとFFTWは、可能であれば、自動的に前回のプランの情報を再利用します。(一方で、グルインターフェースでは、既に生成済みのプランを他の配列に用いることができますが、十分気をつける必要があります。4章の[[FFTW Reference>>FFTW Reference]]を参照のこと)
プランを用いて計算が終了したら、fftw_destroy_plan(plan)を呼び出すことによって、メモリを解放します。
void fftw_destroy_plan(fftw_plan plan);
fftw_mallocで確保された配列は、通常のfreeでなく、fftw_freeで開放してください。(ましてや、deleteなどしてはいけません)
離散フーリエ変換DFTの結果は、配列outに順番に格納され、周波数がゼロの直流成分(DC)がout[0]に入ります。もし、入出力配列が異なる配列なら(in!=out)、入力配列は変更されません。
そうでなければ、入力の配列は変換によって上書きされます。
ユーザーは&color(red){FFTWがノーマライズしていないDFTを計算する}ことに注意する必要があります。このようにして、逆変換のあと正変換(または、正変換のあと逆変換)した場合、元の配列がn倍された結果が得られます。
DFTの定義については、4.7のWhat FFTW Really Computesを参照してください。
もし、gccなどの最近のC99スタンダードをサポートするCコンパイラが利用可能な場合、以下のように、complex.hをfftw3.hの前にインクルードすれば、fftw_complexは倍精度の複素数型になり、通常の演算で扱うことができます。
#include <complex.h>
#include <fftw3.h>
その他の場合では、FFTWはC99と互換性のある独自の型fftw_complexを定義します。4.1.1 Complex numbersも参照のこと(C++のテンプレートクラスも、型をキャストすることによって使えます)
single(float)と、long-doubleの精度のバージョンもインストールされていて、それらを使うためには、fftw_というプレフィクスをfftw_fやfftwl_に置き換え、リンク時に -lfftw3fや -lfftw3lを用いてください。ヘッダーファイルは同じもの(fftw3.h)で構いません。
FFTW_MEASUREとFFTW_ESTIMATE以外にも多数のフラグがあります。例えば、FFTW_PATIENTは、プラン生成に相当程度長い時間がかかる代わりに、実に速くプランが実行できます(4章の[[FFTW Reference>>FFTW Reference]]を参照)
プランを今後の為に保存することもできて、それについては3.3の[[Words of Wisdom-Saving>>Words of Wisdom-Saving]]を参照してください。
* 複素１次元離散フーリエ変換 Complex One-Dimensional DFTs
プランとは、偶然の結果を達成するための最善の方法をいろいろ思い悩むことである。
(Ambrose Bierce, The Enlarged Devil’s Dictionary.)
----
FFTWで、サイズNのデータの一次元DFTを計算する方法はいたって簡単で、おおよそ以下のようなコードになります。
#include <fftw3.h>
...
{
fftw_complex *in, *out;
fftw_plan p;
...
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
...
fftw_execute(p); /* 必要なだけ繰り返す */
...
fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);
}
コンパイル時には、fftw3ライブラリをリンクしなければなりません。つまり例えばUnix系では -lfftw3 -lmをリンカオプションとしてつけます。
まず初めに、入出力の配列のメモリを確保します。
メモリは好きなように確保して構いませんが、fftw_mallocを使うことを推奨します。
fftw_mallocはmallocのように動作しますが、違いとしてはSSEやAltivecなどのSIMD命令が利用出きるように配列が確保されます。(3.1.1節 SIMD alignment and fftw mallocを参照)
データはfftw_complex型の変数の配列になっており、fftw_complexはデフォルトでdouble二つ分の配列(double[2])です。それぞれ実部(in[i][0])、虚部(in[i][1])に対応します。
次のステップは、プランの生成です。プランはFFTWがFFTを計算する際に必要な全ての情報を含んだものです。以下の関数がプランを作ります。
fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
最初の引数nは計算しようとしている変換のサイズです。nは任意の正の整数をとりえますが、サイズが小さい素数のみを因数に持つ場合、最も効率よく計算できます。(しかしどんな素数であれ,O(nlogn)アルゴリズムを採用しています)
計算時間の次元が、データ数nに対して、nlogn
次の二つの引数は、入出力の配列へのポインタです。これらのポインタは同じものを入れることも可能で、この場合は入力配列が書き換わります。
四番目の引数signは、正変換か逆変換かを表す、FFTW_FORWARD(-1)かFFTW_BACKWARD(+1)のどちらかを入れます。技術的にはこれは、変換の指数部になります。
最後の引数flagは、通常FFTW_MEASUREかFFTW_ESTIMATEをとります。FFTW_MEASUREは実際にいくつかのFFTWを実行して実行時間を計り、サイズnの変換を計算するのに一番良い方法を選択します。このプロセスは少し時間がかかり、(通常、2、3秒)これは変換のサイズnや計算機環境に依存します。FFTW_ESTIMATEは、実際に実行することなく、おそらく準最適だと思われる方法を選択します。 つまり、同じサイズの変換を複数回行い、かつ初期化の時間が重要でない場合はFFTW_MEASUREを使い、そうでなければFFTW_ESTIMATEを利用すれば良いことになります。
inやoutといった入出力データの配列は、FFTW_MEASUREのプランニングを行う際に上書きされるので、FFTW_MEWSUREでのプランニングを行う場合は、入力がユーザーに初期化される前にプランニングを行う必要があります。
一旦プランが生成されれば、特定のint/out配列に対して実際の変換をfftw_execute(plan)で計算する際に、プランは何度でも再利用できます。
ちなみにfftw_executeは,
void fftw_execute(const fftw_plan plan);
です。
もし、同じサイズの違う配列を変換したい場合は、fftw_plan_dft_1dを使って新しいプランを作れます。するとFFTWは、可能であれば、自動的に前回のプランの情報を再利用します。(一方で、グルインターフェースでは、既に生成済みのプランを他の配列に用いることができますが、十分気をつける必要があります。4章の[[FFTW Reference>>FFTW Reference]]を参照のこと)
プランを用いて計算が終了したら、fftw_destroy_plan(plan)を呼び出すことによって、メモリを解放します。
void fftw_destroy_plan(fftw_plan plan);
fftw_mallocで確保された配列は、通常のfreeでなく、fftw_freeで開放してください。(ましてや、deleteなどしてはいけません)
離散フーリエ変換DFTの結果は、配列outに順番に格納され、周波数がゼロの直流成分(DC)がout[0]に入ります。もし、入出力配列が異なる配列なら(in!=out)、入力配列は変更されません。
そうでなければ、入力の配列は変換によって上書きされます。
ユーザーは&color(red){FFTWがノーマライズしていないDFTを計算する}ことに注意する必要があります。このようにして、逆変換のあと正変換(または、正変換のあと逆変換)した場合、元の配列がn倍された結果が得られます。
DFTの定義については、4.7のWhat FFTW Really Computesを参照してください。
もし、gccなどの最近のC99スタンダードをサポートするCコンパイラが利用可能な場合、以下のように、complex.hをfftw3.hの前にインクルードすれば、fftw_complexは倍精度の複素数型になり、通常の演算で扱うことができます。
#include <complex.h>
#include <fftw3.h>
その他の場合では、FFTWはC99と互換性のある独自の型fftw_complexを定義します。4.1.1 Complex numbersも参照のこと(C++のテンプレートクラスも、型をキャストすることによって使えます)
single(float)と、long-doubleの精度のバージョンもインストールされていて、それらを使うためには、fftw_というプレフィクスをfftw_fやfftwl_に置き換え、リンク時に -lfftw3fや -lfftw3lを用いてください。ヘッダーファイルは同じもの(fftw3.h)で構いません。
FFTW_MEASUREとFFTW_ESTIMATE以外にも多数のフラグがあります。例えば、FFTW_PATIENTは、プラン生成に相当程度長い時間がかかる代わりに、実に速くプランが実行できます(4章の[[FFTW Reference>>FFTW Reference]]を参照)
プランを今後の為に保存することもできて、それについては3.3の[[Words of Wisdom-Saving Plans>>Words of Wisdom-Saving Plans]]を参照してください。
