gnuplotとFFTW
勉強のため、gnuplot（gnuplot homepage）とc言語のFFTW（Fastest Fourier Transform in the West.）ライブラリの環境を作ってみます。
gnuplot
gnuplotはコマンドラインから使えるグラフィックユーティリティです。科学者や学生向けに数学の関数をヴィジュアル化する目的で作られました。
homebrewでインストールします。こちらを参考にしました。（gnuplotをHomebrewからインストールするときの手順）
$ brew install gnuplot --with-aquaterm
aquatermはGUIです。gnuplotの出力を表示するために、インストール時点でオプションを追加します。
インストール完了後gnuplotコマンドで対話環境が起動します。まずはaquatermでsin(x)を表示してみましょう。
gnuplot> set terminal aqua gnuplot> plot sin(x)
出力方式、出力先を設定することで、画像ファイルやhtmlのcanvasとして吐き出すことができます。
gnuplot> set terminal png gnuplot> set output 'sin_x.png' gnuplot> plot sin(x)
関数をデータとして吐き出したい場合は、
gnuplot> set table "gnuplot.sample" gnuplot> set yrange [0:1] gnuplot> set xrange [-5:5] gnuplot> plot exp(- x * x / 2) / sqrt(2 * pi) gnuplot> unset table
とすると、gnuplot.sampleファイルに関数の値が書き込まれます。
# Curve 0 of 1, 100 points # Curve title: "exp(- x * x / 2) / sqrt(2 * pi)" # x y type -5 1.48672e-06 i -4.89899 2.45106e-06 i -4.79798 3.99989e-06 i -4.69697 6.46117e-06 i -4.59596 1.0331e-05 i -4.49495 1.6351e-05 i -4.39394 2.56161e-05 i -4.29293 3.97238e-05 i -4.19192 6.09759e-05 i -4.09091 9.26476e-05 i -3.9899 0.000139341 i -3.88889 0.00020744 i -3.78788 0.000305686 i -3.68687 0.00044589 i -3.58586 0.000643795 i -3.48485 0.000920105 i -3.38384 0.00130165 i -3.28283 0.00182273 i -3.18182 0.0025265 i -3.08081 0.00346644 i -2.9798 0.00470779 i -2.87879 0.00632878 i -2.77778 0.00842153 i -2.67677 0.0110926 i -2.57576 0.0144624 i -2.47475 0.0186646 i -2.37374 0.0238433 i -2.27273 0.0301496 i -2.17172 0.0377369 i -2.07071 0.0467541 i -1.9697 0.057338 i -1.86869 0.069604 i -1.76768 0.0836362 i -1.66667 0.0994771 i -1.56566 0.117117 i -1.46465 0.136486 i -1.36364 0.157443 i -1.26263 0.179775 i -1.16162 0.20319 i -1.06061 0.227324 i -0.959596 0.251742 i -0.858586 0.275953 i -0.757576 0.299423 i -0.656566 0.32159 i -0.555556 0.341892 i -0.454545 0.359787 i -0.353535 0.374774 i -0.252525 0.386423 i -0.151515 0.394389 i -0.0505051 0.398434 i 0.0505051 0.398434 i 0.151515 0.394389 i 0.252525 0.386423 i 0.353535 0.374774 i 0.454545 0.359787 i 0.555556 0.341892 i 0.656566 0.32159 i 0.757576 0.299423 i 0.858586 0.275953 i 0.959596 0.251742 i 1.06061 0.227324 i 1.16162 0.20319 i 1.26263 0.179775 i 1.36364 0.157443 i 1.46465 0.136486 i 1.56566 0.117117 i 1.66667 0.0994771 i 1.76768 0.0836362 i 1.86869 0.069604 i 1.9697 0.057338 i 2.07071 0.0467541 i 2.17172 0.0377369 i 2.27273 0.0301496 i 2.37374 0.0238433 i 2.47475 0.0186646 i 2.57576 0.0144624 i 2.67677 0.0110926 i 2.77778 0.00842153 i 2.87879 0.00632878 i 2.9798 0.00470779 i 3.08081 0.00346644 i 3.18182 0.0025265 i 3.28283 0.00182273 i 3.38384 0.00130165 i 3.48485 0.000920105 i 3.58586 0.000643795 i 3.68687 0.00044589 i 3.78788 0.000305686 i 3.88889 0.00020744 i 3.9899 0.000139341 i 4.09091 9.26476e-05 i 4.19192 6.09759e-05 i 4.29293 3.97238e-05 i 4.39394 2.56161e-05 i 4.49495 1.6351e-05 i 4.59596 1.0331e-05 i 4.69697 6.46117e-06 i 4.79798 3.99989e-06 i 4.89899 2.45106e-06 i 5 1.48672e-06 i
ここでは一列目がx軸、二列目がy軸になっています。三列目はちょっと謎...
このデータを使ってグラフを表示するのも簡単にできます。
gnuplot> plot "gnuplot.sample"
FFTW
FFTWは離散フーリエ変換を計算するためのc言語のライブラリです。変換する関数は次元、入力サイズ、実数or虚数別が選べます。
$ brew install fftw
FFTWには離散フーリエ変換を扱うための関数や複素数データ型が用意されています。
FFTW does not use a fixed algorithm for computing the transform, but instead it adapts the DFT algorithm to details of the underlying hardware in order to maximize performance.
FFTWでは入力データに対して、毎回決まった変換アルゴリズムを適応するのではなく、最適化した変換アルゴリズムを用いて離散フーリエ変換を行う仕組みになっています。そのため、変換実行前に「プラン」というものを作成します。複素一次元離散フーリエ変換のプランを作成する関数を見てみましょう。
fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);
fftw_plan_dft_1dはfftw_planを戻り値にもつ関数です。引数として次のものをとります；
int n： データサイズfftw_complex *in：入力データ配列fftw_complex *out：出力データ配列int sign：フーリエ変換（-1）、逆フーリエ変換（+1）かの設定。定数が用意されています。unsigned flags：通常はFFTW_ESTIMATEをつかい、サイズが同じ変換を何回もし、初期化コスト(?)を気にしないならばFFTW_MEASUREを使うみたいです。
感想
力が尽きてしまい、フーリエ変換前後のグラフを出力することは断念しました。sinをフーリエ変換して、一本だけ立つことを確認するとか、関数の和がどんどん短形関数に近づいていく様をみるとか、やりたいです。
