約 2,247,622 件
https://w.atwiki.jp/amaeda/pages/21.html
/* * test of self compiled fftw3 * 2D version * @author maeda * @date 2008/09/09 */ #include stdio.h #define _USE_MATH_DEFINES #include math.h #include "fftw3.h" #pragma comment( lib, "fftw3.lib" ) #define SIZEX 16 #define SIZEY 8 #define SIZE (SIZEX*SIZEY) int main( void ){ fftw_complex *in = NULL; fftw_complex *out = NULL; fftw_plan p = NULL; int i,j,idx; size_t mem_size = sizeof(fftw_complex) * SIZE; in = (fftw_complex*)fftw_malloc( mem_size ); out = (fftw_complex*)fftw_malloc( mem_size ); if( !in || !out ){ fprintf( stderr, "failed to allocate %d[byte] memory(-.-)\n", (int)mem_size ); return false; } // !! row-major alignment is recommended, but here, column-major. p = fftw_plan_dft_2d( SIZEY, SIZEX, in, out, FFTW_FORWARD, FFTW_ESTIMATE ); // input data creation printf("----- INPUT -----\n"); for( j=0; j SIZEY; j++ ){ for( i=0; i SIZEX; i++ ){ idx = SIZEX*j+i; // column-major alignment in[idx][0] = 1 + 2*sin(2*M_PI*i/SIZEX) + sin(4*M_PI*j/SIZEY); in[idx][1] = 0; printf("%d %d %lf %lf\n", i, j, in[idx][0], in[idx][1] ); } } fftw_execute(p); // output is DC exchanged and scaled. double scale = 1. / SIZE; printf("\n----- RESULT -----\n"); for( j=0; j SIZEY; j++ ){ for( i=0; i SIZEX; i++ ){ idx = SIZEX*j+i; printf("%d %d %lf %lf\n", i, j, out[idx][0]*scale, out[idx][1]*scale ); } } if( p ) fftw_destroy_plan(p); if( in ) fftw_free(in); if( out ) fftw_free(out); return true; }
https://w.atwiki.jp/amaeda/
- 番目のお客様 (^-^)/~~ FFTWのページ このページは、フリーのFFTライブラリFFTWが有効に活用されるように作られたページです。 マニュアルの日本語訳や、使用上の注意、サンプルプログラムがあります。 未完成のところも多いですが、常時更新していきます。 FFTWはバージョン2と3では大幅に仕様が異なるので、ご利用のバージョンにご注意ください。 ご意見や、感想を掲示板のほうへお寄せください。 動かない、分からないなどのお問い合わせは takaidohigasi あっと gmail.com へ！(回答あるかどうかは、保証しかねますが。。。) FFTW3.1.2 マニュアル FFTWとは？？ FFTWの使い方 (簡単な使い方)複素１次元離散フーリエ変換 複素多次元離散フーリエ変換 実データ１次元離散フーリエ変換 実データ多次元離散フーリエ変換 その他重要なことデータの並び 多次元配列のフォーマット 知の利用-プランの保存 知の利用上の注意 FFTWリファレンス (詳細なマニュアル) データタイプとファイル プランの利用 ベーシックインターフェース サンプルプログラム programs 掲示板 BBS リンク 本家FFTW 本家マニュアル(英語) 修正/更新 実データ１次元離散フーリエ変換の対称性を利用した、計算機資源節約の節で、hermitian redundancyをハミルトニアン冗長性と誤訳していました。正しくはエルミート～。 4.2節プランの利用を追加しました。 4.3節ベーシックインターフェース4.3.2プランのフラグのところまで追加しました。 位相を求めるサンプルプログラムの間違えを修正しました。ご指摘ありがとうございます！
https://w.atwiki.jp/amaeda/pages/12.html
/* * test of self compiled fftw3 * * @author maeda * @date 2008/09/09 */ #include stdio.h #define _USE_MATH_DEFINES #include math.h #include "fftw3.h" #pragma comment( lib, "fftw3.lib" ) #define SIZE 64 int main( void ){ fftw_complex *in = NULL; fftw_complex *out = NULL; fftw_plan p = NULL; int i; size_t mem_size = sizeof(fftw_complex) * SIZE; in = (fftw_complex*)fftw_malloc( mem_size ); out = (fftw_complex*)fftw_malloc( mem_size ); if( !in || !out ){ fprintf( stderr, "failed to allocate %d[byte] memory(-.-)\n", (int)mem_size ); return false; } p = fftw_plan_dft_1d( SIZE, in, out, FFTW_FORWARD, FFTW_ESTIMATE ); // input data creation printf("----- INPUT -----\n"); for( i=0; i SIZE; i++ ){ in[i][0] = 1 + 2*sin(2*M_PI*i/SIZE) + sin(4*M_PI*i/SIZE); in[i][1] = 0; printf("%d %lf\n", i, in[i][0] ); } fftw_execute(p); // output is DC exchanged and scaled. double scale = 1. / SIZE; printf("\n----- RESULT -----\n"); for( i=0; i SIZE; i++ ){ printf("%d %lf %lfi\n", i, out[i][0]*scale, out[i][1]*scale ); } if( p ) fftw_destroy_plan(p); if( in ) fftw_free(in); if( out ) fftw_free(out); return true; }
https://w.atwiki.jp/amaeda/pages/20.html
チュートリアル Tutorial この章では、FFTWの基本的な使い方、つまりどうやって単一の配列のフーリエ変換を求めるのか、について述べます。 この章では、全ての事柄に対しては言及しません。 詳しくいえば、ここに載っていない追加的なルーチンやフラグがFFTWには実装されています。それについては、追加的になにができるのかが記してある場所をできる限り、明示したいとおもいます。より詳細な情報については4章のFFTW Referenceを参照してください。 プログラムを使用するためには、FFTWのインストールとコンパイルが必要です。インストールの詳細については、8章のInstallation and Customizationを参照してください。 チュートリアルは、順番に読んでいくことをおすすめします。例え、最初の節で説明される一次元DFT以外に興味を持っていたとしても、少なくとも、他の節を読む前に、最初の節(Complex One-Dimensional DFTs)は目を通してください。 FFTWのバージョン2以前のユーザーは7章のUpgrading FFTW version2を参照すると良いでしょう。 複素１次元離散フーリエ変換(Complex One-Dimensional DFTs) 複素多次元離散フーリエ変換(Complex Multi-Dimensional DFTs) 実データ１次元離散フーリエ変換(One-Dimensional DFTs of Real Data) 実データ多次元離散フーリエ変換(Multi-Dimensional DFTs of Real Data)
https://w.atwiki.jp/amaeda/pages/49.html
FFTWマニュアル - FFTWリファレンス FFTWリファレンス FFTW Reference この章では、連続（シングルプロセッサの）FFTWの関数すべての、完全なマニュアルを扱います。並列変換については、5章のParallel FFTWをご覧くさい データタイプとファイル 複素数 精度 メモリの確保 プランの利用 ベーシックインターフェース
https://w.atwiki.jp/amaeda/pages/14.html
Introduction このマニュアルはFFTW(the Fastest Fourier Transform in the West)バージョン3.1.2のものです.FFTWは離散フーリエ変換を計算する包括的な高速Cルーチン集であり、以下のような特徴があります. 複素数データ,実数データ,偶対称,奇対称な実数データのDFT(*1),実数データの離散Hartley変換(DHT)を計算可能 入力データは任意のサイズがとれる.O(nlogn)アルゴリズムを,素数を含む,すべてのサイズで採用. 計算の時間の次元がデータ数nに対して,nlogn 任意の多次元データをサポート SSE/SSE2/3DNow!/Altivecの命令セットをサポート FFTW3.1.2は共有メモリシステムに対しては、並列(マルチスレッド)変換が可能。FFTW3.1.2は分散メモリ並列変換ができないが、MPIバージョンをじきに実装予定 当面の間は,FFTW2.1.3に実装されたMPIが利用可能 ここでは、どんなアプリケーションに利用するかに関わらず、DFTの性質および使い方を熟知しているものとして話を進めます。もしそうでなければ、例えば、"The Fast Fourier Transform and Its Applications"(*2),を参照してください. 我々のWebページにもFFT関係の情報のリンク集があります。 FFTWを効果的に使うためには、FFTWの内部構造の基本的なコンセプトについて知っている必要があります。すなわち、FFTWはDFT変換を計算する時に固定のアルゴリズムではなく、パフォーマンスを最大化するために、細かいハードウェア条件に合わせたアルゴリズムを用います。つまり、計算は大きくわけて2段階にわけることができます。 第一にFFTWのプランナーはあなたのマシンで計算する時にもっとも速くなる方法を”学習”します。プランナーはこの情報を含んだ「プラン」と呼ばれるデータ構造体を生成します。 続いて、プランによって示された入力データ配列を変換するために、プランが実行されます。 このプランは何回でも再利用可能です。高性能なアプリケーションでは、同じサイズの変換が何度も計算されます。そのため、初期化の計算コストが比較的高くても、許容可能でしょう。 一方で、特定のサイズの一回限りの変換を行いたい時には、プラン生成にかかる時間が重要になります。この場合は、FFTWはヒューリスティックにあるいは前回計算されたプランに基づいて、高速にプランを計算します。 FFTWは任意の長さ、ランク、次元そして、一般的なメモリレイアウトの変換をサポートしています。しかし、簡単なケースではこの様に一般性を持たせることは不必要で、混乱を招くことでしょう。そういうわけで、我々は三つのレベルの一般性に対し、それぞれインターフェースを作りました。 ベーシックインターフェース 連続データの単一の変換を行う場合 アドバンストインターフェース 多次元配列、規則的な配列の変換を行う場合 グルインターフェース 最も一般的なデータ配置、次元、規則を持つ場合の変換を行う場合 我々は、大多数のユーザーがベーシックインターフェースで最も良い結果を得られることを期待しています。グルインターフェースは、諸問題を避けるために注意深くドキュメントを読む必要があります。 プランナーは、自動でプランの最適化を行ってくれる一方で、上級ユーザ向けの機能として、マニュアルでカスタマイズすることも可能です。例えば、コードの大きさが懸案事項である場合には、FFTWをリンクする際に、そのアプリケーションに必要な部分のみをリンクするツールも用意しています。反対に、標準のFFTWでは機能的に不十分で、拡張の必要があるかもしれません。例えば、標準のFFTWでは、配列サイズが、小さい素数(2,3,5,7)で素因数分解される場合に最も効率的に機能します。一方で、その他の場合には、一般的な利用を想定した遅いルーチンを使います。もし、その他の配列サイズに対し、効率的な計算をする必要があれば、FFTWのコードジェネレータが役に立ちます。コードジェネレータは高速なCプログラム(コードレット)を任意の配列サイズに対し生成します。例えば、513(=19*33)のサイズの変換がしたい場合、FFTWを因子19が効率的に計算できるようカスタマイズできます。 FFTWに関する詳細の情報に関しては、"FFTW An adaptive software arrchitecture for the FFT”(*3)や"The Fastest Fourier Transform in the West"(*4)を参照してください。 コードジェネレータについては、"A fast Fourier transform compiler"(*5)を参照してください。 これらの論文は、最新バージョンのFFTW,FAQ,性能テスト,その他のリンクと合わせてFFTWのホームページ(*6)から入手可能です。 FFTWの現在のバージョンは、過去30年間にわたるFFTの文献の様々なアイデアを反映しています。 FFTWには、何らかの形で、Cooley-Tukeyアルゴリズム(素因数のアルゴリズム)、Raderのアルゴリズム(素数サイズのアルゴリズム)、そしてsplit-radixアルゴリズム(Dan Bernsteinによる改良版)が盛り込まれてます。FFTWのコードジェネレータも我々は完全に理解できていないアルゴリズムを実装しています。適宜、引用文献のこと。 (*1)これらの対称な変換はそれぞれ、一般的に離散コサイン/サイン変換として知られている (*2) "The Fast Fourier Transform and Its Applications",E. O. Brigham, Prentice-Hall, Englewood Cliffs, NJ, 1988 (*3) FFTW An adaptive software architecture for the FFT,M. Frigo and S. G. Johnson, 23rd International Conference on Acoustics, Speech, and Signal Processing (Proc. ICASSP 1998 3, p. 1381) (*4) The Fastest Fourier Transform in the West, M. Frigo and S. G. Johnson, technical report MIT-LCS-TR-728 (Sep. ’97) (*5) A fast Fourier transform compiler, M. Frigo, in the Proceedings of the 1999 ACM SIGPLAN Conference on Programming Language Design and Implementation (PLDI), Atlanta, Georgia, May 1999 (*6) http //www.fftw.org/
https://w.atwiki.jp/amaeda/pages/17.html
全般 # スケーリング 変換の結果は配列のスケール倍される。 プラン生成 プラン生成は、原則的にデータの”メモリ空間”に対して行われ、 (これはFFTWのバージョン３になってから大きく変更されたポイントである。) プラン生成時には、計算時間の最適化などを行うので、多少時間がかかる。 したがって、繰り返し計算を行う時には、一度生成したプランを繰り返し使うことが推奨されている。 すると、例えば、同じサイズの入力データが複数あるとき、一度プランを生成してから、プラン生成時に指定したメモリ空間にデータをコピーして演算を行うことが考えられる。 この際、プラン用のメモリ空間が、”バッファ”として余分に消費されることになるので、多くのメモリを必要とする場合には検討が必要だ。 また、ここまででプランは特定の”メモリ空間”に対して生成されると述べたが、データサイズが同じ、などの制約下であれば、他のメモリ空間に対してプランを再利用する関数も存在する。 また、プラン生成時に入出力配列は上書きされることがある。フラグの設定をよく確認しよう。面倒だったら、入力配列はプラン生成後の初期化するよう心がけよう。 データのメモリ空間上の並び(行メジャー) FFTWのデータはrow-major(行メジャー)で考えられて作られている。 行メジャーとは、最後の配列が、もっともはやく変化する順序である。 int a[SIZEX][SIZEY]; というものがあったとすると、メモリ空間では、 a[0][0],a[0][1],・・・,a[1][0],・・・,a[SIZEX-1][SIZEY-1] のように並ぶ。 利用上は、column-major(列メジャー)でも構わないが注意が必要である。 すると、利用者は、int a[SIZEY][SIZEX]のつもりなのに、FFTWは、int a[SIZEX][SIZEY]だと思っているので、 のようにパラメータのヒントなどで誤解を招く表現が生じる。 要注意である。 このページでは、画像データなどを扱うことを考え、FFTWの推奨ではないのを承知で、列メジャーのデータを扱うコードに統一することとする。 行メジャー(row-major)ーと列メジャー(column-major) 実数- 複素数/複素数- 実数変換
https://w.atwiki.jp/amaeda/pages/31.html
FFTWマニュアル - データの並び - SIMDの並びとfftw_malloc SIMDのならびとfftw_malloc SIMD alignment and fftw_malloc "Single Instruction Multiple Data"(単一命令複数データ流)を意味するSIMDは、 いくつかの、単一命令で複数回(通常２回または４回)同時に操作を行えるプロセッサでサポートされています。SIMD浮動小数点命令は、いくつかのポピューラーなCPUで使えます。たとえば、SSE/SSE2(単精度/倍精度)は、Pentium Ⅲ/Ⅳ以上、3DNow!(単精度)はAMD K7以上、そしてAltiVec(単精度)はPowerPC G4以上で利用できます。FFTWは、これらのいかなるシステムの対してもSIMD命令をサポートするようにコンパイルできます。 FFTWライブラリをSIMDサポートをするようにコンパイルして、リンクすると、ほとんどの複素数変換およびr2c/c2r変換で、劇的なスピードの向上が得られます。しかしながら、このスピードアップを得るためには、FFTWにわたされる複素数(または実数)データの配列は、メモリ空間で特別にならんでなければいけません(典型的には16-byte揃えです)。そして、しばしばこのならびは、mallocなどの通常のメモリ確保ルーチンよりも制約がきびしくなっています。 したがいまして、SIMDに適切な並びを保証するためには、変換に用いるデータをfftw_mallocで確保し、fftw_freeで解放するようにすることをお勧めします。これらの関数はmalloc/freeと、返り値のポインタが、必要とされる並びを保証していること以外は、まったく同じインターフェースおよび動作をします。(これは、memalignだとか、個々のOSのそれに相当する関数を呼ぶことで実現されています) しかし、絶対にfftw_mallocでメモリを確保しなければならないわけではありません。あなたの好きなように、mallocやnew(C++)や、固定サイズの配列の宣言によって確保してもよいでしょう。もし配列が適切に並んでいないようであれば、FFTWはSIMDの拡張命令を使わなくなります。
https://w.atwiki.jp/amaeda/pages/16.html
複素１次元離散フーリエ変換 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_destroy_plan(plan)を呼び出すことによって、メモリを解放します。 void fftw_destroy_plan(fftw_plan plan); fftw_mallocで確保された配列は、通常のfreeでなく、fftw_freeで開放してください。(ましてや、deleteなどしてはいけません) 離散フーリエ変換DFTの結果は、配列outに順番に格納され、周波数がゼロの直流成分(DC)がout[0]に入ります。もし、入出力配列が異なる配列なら(in!=out)、入力配列は変更されません。 そうでなければ、入力の配列は変換によって上書きされます。 ユーザーは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を参照) プランを今後の為に保存することもできて、それについては3.3のWords of Wisdom-Saving Plansを参照してください。
https://w.atwiki.jp/ketcindy/pages/103.html
逆変換を説明する図を描く． inverse.zip ／／ 点Aと点Bをとり，中間あたりに点Cをとる． Bspline("1",[A,C,B]); ／／ Bspline曲線を描く．点Cを動かして調節する． Arrowhead(A,A-C); ／／ 矢じりをつける． Ptsize(5); Drawpoint([A,B]); ／／ 大き目の点を打つ． Expr([A,"se","\mathrm{P}(x,\ y)"]); Expr([B,"n","\mathrm{P}`(x`,\ y`)"]); Expr([C,"se","f^{-1}"]); ／／ 点や変換の名前を書く． Setcolor([0,0,1]); Bspline("2",[A,D,B]); Arrowhead(B,B-D); Expr([D,"nw","f"]); Setcolor([0,0,0]); ／／ 青色で描く．