はじめに
Pythonの線形回路解析パッケージであるLcapyを用いて回路の伝達関数を記号解析する方法について記事を書きました。
【回路×Python】Lcapyを使って回路の伝達関数を求める方法
伝達関数を求めたあとは、それを展開したり、数値解析したくなります。
これら方法についてこの記事では解説します。
環境
Python: 3.7.4、SymPy: 1.6.2、Lcapy: 0.67.0
解析する回路
オペアンプを抵抗RG, RFでフィードバックした反転増幅器について考えます。
オペアンプの特性は以下の式で表されるものとします。
$$A(f) = \frac{A_{OL}}{1+s/p}$$
この回路は前回の記事にも登場したもので、LTSPICE上では以下のように描きます。
伝達関数を求めるコードと実行結果は以下の通りです。
from lcapy import *
cct = Circuit("""
E1 0 OUT N002 0 1
C1 N002 0 {1/p}
R1 N002 0 1
R2 OUT INN RF
R3 INN N001 RG
V1 N001 0 s 1
G1 0 N002 0 INN AOL
""")
H = cct["OUT"].V(s).simplify()
H
実行結果：
$$- \frac{A_{\mathrm{OL}} R_{F} p}{A_{\mathrm{OL}} R_{G} p + R_{F} p + R_{G} p + s \left(R_{F} + R_{G}\right)}$$
１．伝達関数を展開する
伝達関数を展開する方法は公式ドキュメントのExpressions => Methodsに記載があります。
上の回路で使用した例を以下に示します。
このように、非常に少ないコードで簡単に式を展開・変形できます。
２．数値計算する
回路を記号解析して伝達関数を求めた後は、伝達関数に値を代入して数値解析してみましょう。
まず伝達関数を求めます。（コードは上に載せたものから一部修正しています。理由は下参照）
from lcapy import *
cct = Circuit("""
E1 0 OUT N002 0 1
C1 N002 0 {1/p}
R1 N002 0 1
R2 OUT INN RF
R3 INN N001 RG
V1 N001 0 step
G1 0 N002 0 INN AOL
""")
H = (cct["OUT"].V(s) / cct.V1.V(s)).simplify()
H
実行結果：
$$- \frac{A_{\mathrm{OL}} R_{F} p}{p \left(A_{\mathrm{OL}} R_{G} + R_{F} + R_{G}\right) + s \left(R_{F} + R_{G}\right)}$$
求めた伝達関数に値を代入するのは以下のコードです。
入力した値：AOL= 1000、p = 2π x 10kHz, RG=1kΩ、RF=2kΩ
H = H.subs("AOL", 1000).subs("p", 2*pi*1e4).subs("RG", 1000).subs("RF", 2000)
H
実行結果：
$$- \frac{40000000000.0 \pi}{3000 s + 20060000000.0 \pi}$$
伝達関数の特性（ここでは強度、単位dB）は以下のコードで計算・グラフ化できます。
import numpy as np
import matplotlib.pyplot as plt
freq = np.logspace(start=5, stop=9, num=(9-5)*21)
sm1 = H(f).dB.evaluate(freq)
fig, ax = plt.subplots()
ax.plot(freq,sm1)
ax.set_xscale("log")
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Signal Magnitude (dB)")
plt.show()
実行結果：
数値計算して配列を生成しているのは以下のコードです
sm1 = H(f).dB.evaluate(freq)
ここでは以下のような処理が行われています。
sに2πjfを代入 → これのデシベルを計算 → evaluate(freq)で数値計算
バグ？(2020/10/31)
「解析する回路」に記載したコードから数値計算をすると、エラーが発生します。
from lcapy import *
cct = Circuit("""
E1 0 OUT N002 0 1
C1 N002 0 {1/p}
R1 N002 0 1
R2 OUT INN RF
R3 INN N001 RG
V1 N001 0 s 1
G1 0 N002 0 INN AOL
""")
H = cct["OUT"].V(s).simplify()
H = H.subs("AOL", 1000).subs("p", 2*pi*1e4).subs("RG", 1000).subs("RF", 2000)
import numpy as np
import matplotlib.pyplot as plt
freq = np.logspace(start=5, stop=9, num=(9-5)*21)
sm1 = H(f).dB.evaluate(freq)
fig, ax = plt.subplots()
ax.plot(freq,sm1)
ax.set_xscale("log")
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Signal Magnitude (dB)")
plt.show()
""")
実行結果：
ValueError: Cannot convert non-causal s-expression to f domain
コード的にはおかしくないと思うのですが、エラーが発生します。
「２．数値計算する」のコードでも同じ伝達関数が得られるのですが、そちらは何故か計算できます。
私の環境では公式ドキュメントに記載されているコードでもエラーが発生しているので、バグじゃないかと思っています。
状況が変化したらこの箇所は修正します。