フーリエ級数展開のデモンストレーションをPythonで書いた話

python

はじめに

東京大学の小山先生が、フーリエ級数展開のデモンストレーションをMATLABでお書きになった。

この素晴らしいアニメーションをPythonで再現するスクリプトを書いても良いのではないかと思い、今回の表題に至るわけである。
ちなみに再現したアニメーションは以下の通りである。グラフの軸ラベルがずっと固定であったり、描画範囲が微妙に異なるので完全再現ではないが、それなりに再現できていると思われる。

スクリプトの解説（描画まわり）

小山先生のMATLABスクリプトを参考に、以下のPythonスクリプトを書いた。
https://gist.github.com/tam17aki/4dc145b3f7d128fdb15ebf0e9137cd9b

基本的にはオリジナルのMATLABスクリプトを逐次Pythonに翻訳していく方針である。配列操作にはnumpy、またグラフの描画とアニメーション作成にはmatplotlibというライブラリを用いた。

アニメーションを作成する以前に、まず複数のグラフを一枚のキャンバスに重ね描きしないといけない。matplotlibはデフォルトで重ね描きが有効になっているので、MATLABのようにhold onとhold offを切り替える必要はない。

キャンバスを用意するにはまず以下を書く。

import matplotlib.pyplot as plt
fig = plt.figure(figsize=[8.0, 8.0])  # 800 x 800

これで図の外枠ができ上がるので、add_axes関数により実際のグラフ描画枠を配置していく。描画位置を指定するsubplotをより細かく制御する働きをしていると考えればよい。

ax1 = fig.add_axes([0.04, 0.85, 0.14, 0.14])
ax2 = fig.add_axes([0.24, 0.85, 0.74, 0.14])
ax3 = fig.add_axes([0.04, 0.65, 0.14, 0.14])
ax4 = fig.add_axes([0.24, 0.65, 0.74, 0.14])
...

add_axesの引数は「x軸の開始位置, y軸の開始位置, x軸の長さ, y軸の長さ」を、全体に対する比率で書いて与える。詳細な仕様は公式マニュアルを参考にしたほうが良い（URL）。

グラフのプロットは各描画枠（axes）ごとに行うことができる。例えば

ax2.plot(time_axis, sig, linestyle="-", color="b", linewidth=1.5)
ax4.plot(time_axis, coef[0] * numpy.sin(phi[:, 0]),
         linestyle="-", color="b", linewidth=1.5)

のような形である。

重ね描きは各axesごとにplotを並べるだけで実現できる。例えば以下のようにする。

ax3.plot(coef[0] * numpy.cos(ANGLE), coef[0] * numpy.sin(ANGLE),
         color="k", linewidth=1.5)
ax3.plot([0, circ[0, 0]], [0, circ[1, 0]],
         linestyle="-", color="b", marker="o",
         markerfacecolor="b", markersize=4)
ax3.plot([circ[0, 0], 0.6], [circ[1, 0], circ[1, 0]],
         linestyle=":", color="b", marker="o", markersize=4,
         markerfacecolor="b", linewidth=1)

以上の要素を用いることで、アニメーションの1コマに相当するグラフを描くことができるので、あとはそれらをくっつければ動画にすることができる。アニメーションの作成にはmatplotlibのanimationモジュールを用いた。

import matplotlib.animation as animation

1コマに相当する各プロットのグラフを集めて、それをリストに追加していくのが基本的な使い方である。ちなみにplotの戻り値は描画内容を格納する特殊なリストになっているので、「+=」で次々と追加していくことができる（URL）。

images = []  # アニメーションの総プロットを格納するリスト
...
for t0 in range(TIME_NUM):
...
    im = ax1.plot(coef[0] * numpy.cos(ANGLE),
                  coef[0] * numpy.sin(ANGLE),
                  color="k", linewidth=1.5)
    im += ax1.plot(coef[1] * numpy.cos(ANGLE) + circ[0, 0],
                   coef[1] * numpy.sin(ANGLE) + circ[1, 0],
                   color="k", linewidth=1.5)
...
    im += ax2.plot(time_axis, sig, linestyle="-", color="b", linewidth=1.5)
...
    images.append(im)  # 1コマ分の全プロットを集めたものをリストに追加

最終的に、以下のArtistAnimation関数を用いてアニメを作成し、save関数でmp4に書き出すことができる。

ANIME = animation.ArtistAnimation(fig, images, interval=40)
ANIME.save("plot_sawtooth.mp4", writer="ffmpeg", dpi=300)

ArtistAnimation関数の引数intervalには各フレーム間の遅延をmsecで与える。save関数の引数dpiには解像度を与える。今回はあまり高解像度にする必要はないだろう。mp4作成にあたって、事前にffmpegをインストールしておく必要がある。

おわりに

今回はmatplotlibの優秀さを実感した。アニメーション作成がやや重たいのが課題である。FuncAnimation関数を用いたとき、軽くなるかを一度検証してみる必要がある。

フーリエ級数展開はいいぞ。