[Circuit x Python] How to expand and calculate transfer functions using Lcapy

Introduction

I wrote an article about how to symbolically analyze the transfer function of a circuit using Lcapy, a linear circuit analysis package in Python.

[Circuit x Python] How to find the transfer function of a circuit using Lcapy

After finding the transfer function, we want to expand it and perform numerical analysis. These methods are explained in this article.

environment

Python: 3.7.4、SymPy: 1.6.2、Lcapy: 0.67.0

Circuit to analyze

Consider an inverting amplifier that feeds back an operational amplifier with resistors RG and RF. The characteristics of the operational amplifier shall be expressed by the following equation. $A(f) = \frac{A_{OL}}{1+s/p}$

This circuit also appeared in Previous article, and is drawn as follows on LTSPICE.

The code to find the transfer function and the execution result are as follows.

python

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

Execution result: $- \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)}$

1. 1. Expand the transfer function

Instructions on how to expand the transfer function can be found in the official documentation Expressions => Methods (http://lcapy.elec.canterbury.ac.nz/expressions.html).

The example used in the above circuit is shown below.

In this way, you can easily expand and transform expressions with very little code.

2. Numerical calculation

After symbolically analyzing the circuit to find the transfer function, let's assign a value to the transfer function and perform numerical analysis.

First, find the transfer function. (The code has been partially modified from the one listed above. See below for the reason)

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

Execution result: $- \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)}$

The following code assigns a value to the obtained transfer function. Input values: 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

Execution result: $- \frac{40000000000.0 \pi}{3000 s + 20060000000.0 \pi}$

The characteristics of the transfer function (intensity, unit dB here) can be calculated and graphed with the following code.

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()

Execution result:

The following code is numerically calculated to generate an array sm1 = H(f).dB.evaluate(freq)

The following processing is performed here. Substitute 2πjf for s → Calculate the decibel of this → Numerical calculation with evaluate (freq)

bug? (2020/10/31)

An error will occur if you perform a numerical calculation from the code described in "Circuit to analyze".

python

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()
""")

Execution result: ValueError: Cannot convert non-causal s-expression to f domain

I don't think it's strange in code, but I get an error. The same transfer function can be obtained with the code of "2. Numerical calculation", but it can be calculated for some reason. In my environment, the code described in the official document also causes an error, so I think it is a bug. If the situation changes, I will correct this part.