読者です 読者をやめる 読者になる 読者になる

ushumpei’s blog

生活で気になったことを随時調べて書いていきます。

gnuplotとFFTW

勉強のため、gnuplotgnuplot 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)

f:id:ushumpei:20160604174532p:plain

出力方式、出力先を設定することで、画像ファイルや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"

f:id:ushumpei:20160616000500p:plain

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_1dfftw_planを戻り値にもつ関数です。引数として次のものをとります;

  • int n: データサイズ
  • fftw_complex *in:入力データ配列
  • fftw_complex *out:出力データ配列
  • int signフーリエ変換(-1)、逆フーリエ変換(+1)かの設定。定数が用意されています。
  • unsigned flags:通常はFFTW_ESTIMATEをつかい、サイズが同じ変換を何回もし、初期化コスト(?)を気にしないならばFFTW_MEASUREを使うみたいです。

感想

力が尽きてしまい、フーリエ変換前後のグラフを出力することは断念しました。sinをフーリエ変換して、一本だけ立つことを確認するとか、関数の和がどんどん短形関数に近づいていく様をみるとか、やりたいです。