がれすたさんのDIY日記

電子回路、Python、組み込みシステム開発、自作エフェクターを語るblog

Python:scipy.signal.bodeでボード線図を引く



今回はScipyのsignal.bodeを使って周波数応答を描画してみる。
使うメソッドは

  • scipy.signal.lti
  • scipy.signal.bode

である。

scipy.signal.lti

線形時不変システムを作成できるメソッド。引数の渡し方によって伝達関数だったり状態空間モデルとかいろいろ定義できる。
詳しくは公式ドキュメントが参考になる。
scipy.signal.lti — SciPy v1.2.0 Reference Guide

G(s) = \frac{10000}{s+10000}

というのを作りたかったら

from scipy import signal

# transferfunc = 10000 / s+10000
s = signal.lti([10000], [1.0,10000])

とすれば良い

scipy.signal.bode

次にボード線図の方だ
これもドキュメント読んでもらえればわかるんだけど(scipy.signal.bode — SciPy v1.2.0 Reference Guide

f:id:gsmcustomeffects:20190113112621p:plain
bodeメソッドの引数、返り値

基本的にLTIで作ったインスタンスを入れてあげればw,mag,phaseが帰ってくる感じ.
見た感じARMAモデルそのまま打ち込んでも表示してくれるみたい

w[radian/sec]はグラフのx軸に当たる部分で正規化周波数で帰って来るのでfにしたい場合は2piで割れば良い

実際の使い方はこうなる

w, mag, phase = signal.bode(s)

bode plot

最後に全体のコードを張っておく

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# transferfunc = 10000 / s+10000
s = signal.lti([10000], [1.0,10000])
w, mag, phase = signal.bode(s)

f = w/(2.0*np.pi)
fig, (ax0, ax1) = plt.subplots(nrows=2, sharex=True)
ax0.semilogx(f, mag, 'b-')
ax0.set_ylabel("Magnitude[dB]")
ax1.semilogx(f, phase, 'r-')
ax1.set_ylabel("Phase[degree]")
ax1.set_xlabel("Freq[Hz]")
plt.show()

f:id:gsmcustomeffects:20190113113807p:plain
bode plot