gnuplot の基本的なテクニック - フーリエ展開の可視化をやってみる

# はじめに
`gnuplot` はサイエンスでも使われる有名なプロットソフトウェアです．最近では `Python`+`matplotlib` などが主流になりつつありますが，`gnuplot` 単体でも多少がんばればそれなりにのプロットを作成することもできるものです．そこで，今回は `gnuplot` 初心者に向けて，`gnuplot` がちょっと便利になる基本的なテクニックを紹介します．

## 目標
今回はこんな感じのプロット作成するためのテクニックを紹介します．元ネタは Twitter で見た["これ"](motoneta)です．矩形波をフーリエ級数展開して可視化しています．使用したバージョンは `gnuplot-5.0.0` です．ただし，本稿で使用しているほとんどの機能は `gnuplot-4.6.x` でも使用可能です．

![c.gif](https://qiita-image-store.s3.amazonaws.com/0/46396/1d8f29fc-6fb1-176b-5966-0e09a89e4aaf.gif)

## 再帰関数 or `sum` 記法
矩形波をフーリエ級数展開すると以下の式で表現できます．

```math
f(x) = \sum_{n=1}^\infty \frac{4}{\pi(2n-1)}\sin\left((2n-1)x\right)
```

無限大まで足すわけにはいかないので，これを適当なところまで足した関数 $f_N(x)$ を作成することにします．

```math
f_N(x) = \sum_{n=1}^N \frac{4}{\pi(2n-1)}\sin\left((2n-1)x\right)
```

このように和の形式で定義されている関数を `gnuplot` で用いるためには，まず再帰を利用する手法が考えられます．関数 `f(x,n)` は以下のように定義できます．

```gnuplot
f(x,n) = (n==0)?0:4.0/pi/(2*n-1)*sin((2*n-1)*x)+f(x,n-1)
```

![a.png](https://qiita-image-store.s3.amazonaws.com/0/46396/61d3e926-c7ac-e6f3-ee9f-b756417a0995.png)

この定義では `n` に非負整数以外を入れると再帰が止まらずに stack overflow を起こしてしまうことに注意してください．同様の関数は `sum` 記法を用いても定義できます．`sum` 記法は最近導入された手法ですが，現在比較的広く利用されている `gnuplot-4.6` にはすでに導入されています． `sum` 記法を使った場合には以下の表記をします．

```gnuplot
f(x,n) = sum [m=1:n] 4.0/pi/(2*m-1)*sin((2*m-1)*x)
```

数学的な定義とかなり近い表記になりました．どちらを使用してもいいと思います．なお，この `sum` 記法は関数定義以外でも用いることができます．以下は $\pi$ の近似値を計算したものです．

```gnuplot
print sum [i=1:1000] 4.0*(-1)**(i+1)/(2*i-1)
# 3.14059265383979
```

これで矩形波のフーリエ級数展開を表現した関数の定義ができました．次は回転する円の部分に注目しましょう．

## `set parametric`
円周上を移動する点を表現するにはどうしても二価関数を使うことになります．`parametric` を設定して $x,y$ 座標を別々に与えましょう．

```gnuplot
set parametric # dummy 変数が t になる
```

`parametric` を設定するとダミー変数が `t` になります．プロットでは $x,y$ 座標を表す関数(数値)を独立に与えます．$x$ 座標は何でもいいのですが元ネタに合わせて位相を $\pi/2$ だけずらした関数にすることにしました．

```gnuplot
# わかりやすくするために関数を定義
y(t,n) = f(t,n)
x(t,n) = f(t-pi/2.0,n)

set size ratio -1 # x,y 座標を同じスケールにする
set sample 1e3    # サンプル点を 1000 個にする

plot x(t,10),y(t,10) w line
```
以下のような正方形に近い図形がプロットされるはずです．正方形の上を動点が等速で移動することを考えると，動点の $x,y$ 座標が矩形波に対応していることが分かります．

![a.png](https://qiita-image-store.s3.amazonaws.com/0/46396/f3c46ab9-59fc-b417-0c06-294f6040b207.png)


## `do` 記法による繰り返し
それでは動点を動かしてみましょう．そのためには `do` 記法を使うと便利です．`do` 記法も `gnuplot-4.6` 以降であれば問題なく利用できます．

```gnuplot
# 1 点だけプロットするので範囲を明示する必要がある
set xrange [-1.5:1.5]
set yrange [-1.5:1.5]
do for [T=0:49] {
  plot x(2*pi*T/50.,5), y(2*pi*T/50.,5) w p lw 2
}
```

動点が領域内を動き回る様子が確認できたでしょうか？


## gif アニメーション出力
せっかくなので gif アニメーションに書き出してみましょう．`set terminal` で `gif` を指定，オプションに `animate` を指定することで，連続して書きだしたプロットが gi作成したプロット
f アニメーションになります．

```gnuplot
set terminal gif anim
set output "anim.gif"
do for [T=0:49] { plot x(2*pi*T/50.,5), y(2*pi*T/50.,5) w p lw 2 }
unset output
```

![x.gif](https://qiita-image-store.s3.amazonaws.com/0/46396/4ffddcd6-8545-84d2-b54c-e2f5eb4cc9dd.gif)


## `plot` 内部の `for` 記法でまとめてプロット
動点を複数の `n` についてプロットしてみましょう．

```gnuplot
plot \
  x(2*pi*T/100.,1), y(2*pi*T/100.,1) w p, \
  x(2*pi*T/100.,2), y(2*pi*T/100.,2) w p, \
  x(2*pi*T/100.,3), y(2*pi*T/100.,3) w p, \
  x(2*pi*T/100.,4), y(2*pi*T/100.,4) w p, \
  x(2*pi*T/100.,5), y(2*pi*T/100.,5) w p
```

もちろん上記の方法でもいいですが，似たようなプロット命令を何度も繰り返すのは不格好ですし，ちょっとした修正が必要になったときに大変面倒になります．ここでは `plot` 内部で `for` 記法を用いることで表記を簡略化しましょう．

```gnuplot
plot for [n=1:5] x(2*pi*T/100.,n), y(2*pi*T/100.,n) w p
```

`plot` 内 `for` 記法によって比較的スマートに記述することができます．先ほどの `do` 記法と組み合わせることによって多数の動点が一度に動きまわる様子を gif アニメーションで可視化できます．

```gnuplot
set terminal gif amimate
set output "amim.gif"
do for [T=0:99] {
  plot for [n=1:5] x(2*pi*T/50.,n), y(2*pi*T/50.,n) w p pt 1 lw 2 notitle
}
unset output
```

![x.gif](https://qiita-image-store.s3.amazonaws.com/0/46396/d2fafa32-a9d7-2b56-2ff9-64b0efbc27ea.gif)


## まとめ
以上のテクニックを組み合わせてみましょう．冒頭の gif アニメーションは以下のようなコマンドによって生成することができます．`set object` でちょっとした装飾を加えていますが，上記のテクニックを組み合わせるだけで描くことができます．

```gnuplot
set term gif amim size 640,300
set output "anim.gif"

set parametric # ダミー変数は t です
set xrange [-2:12]
set yrange [-2:2]
set trange [0:10]
set size ratio -1
set sample 1e3

# ここでは第 8 項まで足します
N = 8
do for [T=0:99] {
  do for [m=1:N] {
    # 動点に合わせて背景の円を描きます
    set object m circle at x(2*pi*T/100.,m-1),y(2*pi*T/100.,m-1) \
      size 4.0/pi/(2*m-1) fc rgb "gray" fs empty border
    }
    # 横方向のラインを描きます
    set arrow 1 nohead from x(2*pi*T/100.,N),y(2*pi*T/100.,N) \
      to 2.0,y(2*pi*T/100.,N) lc rgb "coral" lw 2
    plot \
      for [n=1:N] \
        x(2*pi*T/100.,n),y(2*pi*T/100.,n) w p pt 7 ps .5 notitle, \
      t+2, y(2*pi*T/100.-t,N) w l lw 3 lc 1 notitle
  }
}
unset output
```

上記のコマンドで作成したプロットを以下に貼り付けます．

![amim.gif](https://qiita-image-store.s3.amazonaws.com/0/46396/94da0151-ffd7-7d7c-4724-37d16730a64d.gif)


## 参考資料
- [uncanny wally on Twitter: "http://t.co/mdM8u3Vis3"](motoneta)
- [gnuplot homepage](gnuplot)
- [Demo scripts for gnuplot version 5](demo)

[https://twitter.com/ehthayer/status/612290883187286016]: motoneta
[http://www.gnuplot.info/]: gnuplot
[http://gnuplot.sourceforge.net/demo_5.0/]: demo