Circuit simulation using Julia (and ACME.jl)

We revisit filter circuits using ACME.jl
filter
vcf
mki-x-esEDU
module
julia
Published

June 15, 2023

1 Intro

Julia is a high-level, general-purpose dynamic programming language. Its features are well suited for numerical analysis and computational science.

I’ve been curious about Julia for a while, and came across a neat package for it called ACME.jl:

“a package for the simulation of electrical circuits, focusing on audio effect circuits”

It has a clever, compact syntax for specifying circuits, and supports a wide range of components. I thought it’d be interesting to run some of the circuits from the mki-x-ESEDU VCF series through it (see Part 1, Part 2, Part 3)!

2 Parsing netlists

I’ve setup a script that converts a standard spiceworld Circuit object to a ACME.jl listing. It first lists all the components, then for each node, lists all the connects that occur at that node. As a simple example the following netlist:

# python code
def rc_1st_order_passive():
    netlist = f"""
        Vin 1 0 saw
        R1 1 2 1000
        C1 0 2 1.0e-6
        J1 2 0
    """

    circuit = parse_netlist(netlist, "RC1stOrderPassive")
    output_path = acme_from_circuit(circuit, "rc_1st_order_passive")

produces the following Julia code:

# julia code
using ACME
using Plots
fs = 44100

function saw(t; freq=10, ppv=10)
    x = @. ppv/2*(2*mod(t*freq, 1)-1)
    return reshape(x, 1, length(x))    
end;

function rc_1st_order_passive(::Type{Circuit}; )
    @circuit begin
    # elements
    jin = voltagesource()
    r1 = resistor(1000)
    c1 = capacitor(1e-06)
    jout_1 = voltageprobe()

    # node connections
    #= node 00 =# jin[-] == c1[1] == jout_1[-] == gnd
    #= node 01 =# jin[+] == r1[1]
    #= node 02 =# r1[2] == c1[2] == jout_1[+]
    end
end

model = DiscreteModel(rc_1st_order_passive(Circuit), 1//fs)

ts = [n/fs for n in 0:10000]
x = saw(ts, freq=25)
y = run!(model, x)
print("Running rc_1st_order_passive complete!\n")

plot(ts, x[1, :], label="V_in", size=(750,350))
plot!(ts, y[1, :], label="V_out", size=(750,350))

Note that it includes some plotting code. We’ll setup some common functions too:

using ACME
using Plots
fs = 44100

function sine(t; freq=10, ppv=10)
    x = @. ppv/2 * sin(2*pi*t*freq)
    return reshape(x, 1, length(x))    
end;

function saw(t; freq=10, ppv=10)
    x = @. ppv/2*(2*mod(t*freq, 1)-1)
    return reshape(x, 1, length(x))    
end;

function square(t; freq=10, ppv=10)
    x = @. ppv/2*(2.0*(mod(t*freq, 1) > 0.5) - 1)
    return reshape(x, 1, length(x))
end;

3 Examples

3.1 Passive RC, 2nd Order

function rc_2nd_order_passive(::Type{Circuit}; )
    @circuit begin
    # elements
    jin = voltagesource()
    r1 = resistor(1000)
    c1 = capacitor(1e-06)
    r2 = resistor(1000)
    c2 = capacitor(1e-06)
    jout_1 = voltageprobe()

    # node connections
    #= node 00 =# jin[-] == c1[1] == c2[1] == jout_1[-] == gnd
    #= node 01 =# jin[+] == r1[1]
    #= node 02 =# r1[2] == c1[2] == r2[1]
    #= node 03 =# r2[2] == c2[2] == jout_1[+]

    end
end

model = DiscreteModel(rc_2nd_order_passive(Circuit), 1//fs)

ts = [n/fs for n in 0:10000]
x = saw(ts, freq=25)
y = run!(model, x)

plot(ts, x[1, :], label="V_in", size=(750,350))
plot!(ts, y[1, :], label="V_out", size=(750,350))

3.2 Active RC, 2nd order

function rc_2nd_order_active(::Type{Circuit}; Rf=5000)
    @circuit begin
    # elements
    j0 = voltagesource()
    r1 = resistor(Rf)
    c1 = capacitor(1e-06)
    opa = opamp()
    r2 = resistor(Rf)
    c2 = capacitor(1e-06)
    jout_1 = voltageprobe()

    # node connections
    #= node 00 =# j0[-] == c1[2] == c2[2] == jout_1[-] == opa["out-"] == gnd
    #= node 01 =# j0[+] == r1[1]
    #= node 02 =# r1[2] == c1[1] == opa["in+"]
    #= node 03 =# opa["in-"] == opa["out+"] == r2[1]
    #= node 04 =# r2[2] == c2[1] == jout_1[+]

    end
end

model = DiscreteModel(rc_2nd_order_active(Circuit), 1//fs)

ts = [n/fs for n in 0:10000]
x = saw(ts, freq=25)
y = run!(model, x)

plot(ts, x[1, :], label="V_in", size=(750,350))
plot!(ts, y[1, :], label="V_out", size=(750,350))

3.3 Fixed Resonance, RC 2nd order

function rc_2nd_order_fixed_resonance(::Type{Circuit}; Rf=10000.0)
    @circuit begin
    # elements
    j0 = voltagesource()
    r1 = resistor(Rf)
    c1 = capacitor(4.7e-08)
    opa1 = opamp()
    r2 = resistor(Rf)
    c2 = capacitor(4.7e-08)
    opa2 = opamp()
    jout_1 = voltageprobe()

    # node connections
    #= node 00 =# j0[-] == c2[2] == jout_1[-] == opa1["out-"] == opa2["out-"] == gnd
    #= node 01 =# j0[+] == r1[1]
    #= node 02 =# r1[2] == c1[2] == opa1["in+"]
    #= node 03 =# opa1["in-"] == opa1["out+"] == r2[1]
    #= node 04 =# r2[2] == c2[1] == opa2["in+"]
    #= node 05 =# c1[1] == opa2["in-"] == opa2["out+"] == jout_1[+]

    end
end

model = DiscreteModel(rc_2nd_order_fixed_resonance(Circuit), 1//fs)

ts = [n/fs for n in 0:10000]
x = saw(ts, freq=25)
y = run!(model, x)

plot(ts, x[1, :], label="V_in", size=(750,350))
plot!(ts, y[1, :], label="V_out", size=(750,350))

3.4 Variable Resonance, RC 2nd order

function rc_2nd_order_variable_resonance(::Type{Circuit}; y=0.5,Rf=500)
    @circuit begin
    # elements
    j0 = voltagesource()
    r1 = resistor(100000)
    r2 = resistor(33000)
    opa1 = opamp()
    r3 = resistor(Rf)
    opa2 = opamp()
    r4 = resistor(Rf)
    opa3 = opamp()
    p = potentiometer(100000, (y == nothing ? () : (y,))...)
    opa4 = opamp()
    r5 = resistor(68000)
    r6 = resistor(100000)
    c1 = capacitor(1e-06)
    c2 = capacitor(1e-06)
    jout_1 = voltageprobe()

    # node connections
    #= node 00 =# j0[-] == r2[2] == p[3] == r5[2] == c2[2] == jout_1[-] == opa1["out-"] == opa2["out-"] == opa3["out-"] == opa4["out-"] == gnd
    #= node 01 =# j0[+] == r1[1]
    #= node 02 =# r1[2] == r2[1] == opa1["in+"]
    #= node 03 =# opa1["in-"] == opa1["out+"] == r3[1]
    #= node 04 =# r3[2] == opa2["in+"] == c1[2]
    #= node 05 =# opa2["in-"] == opa2["out+"] == r4[1]
    #= node 06 =# r4[2] == opa3["in+"] == c2[1]
    #= node 07 =# opa3["in-"] == opa3["out+"] == p[1] == jout_1[+]
    #= node 08 =# p[2] == opa4["in+"]
    #= node 09 =# opa4["in-"] == r5[1] == r6[1]
    #= node 10 =# opa4["out+"] == r6[2] == c1[1]

    end
end


plots = []
for y in [0.1, 0.2, 0.5, 0.8]
    model = DiscreteModel(rc_2nd_order_variable_resonance(Circuit, y=y), 1//fs)
    
    ts = [n/fs for n in 0:10000]
    x = saw(ts, freq=25)
    vout = run!(model, x)
    
    p = plot(ts, dropdims(x, dims=1), label="V_in", title=string("y=", y))
    plot!(ts, dropdims(vout, dims=1), label="V_out")
    push!(plots, p) 
end
plot(plots..., size=(750,350))

Note the unbounded oscillation that grows with ideal (i.e. non voltage clamped) opamps!

3.5 Diode Ladder Full

This produces quite a long code listing, expand below!

Code 
function diode_ladder_full(::Type{Circuit}; ygain=0.5,yoffset=0.5,ycv1=0.5,ycv2=0.5,yr=0.5)
    @circuit begin
    # elements
    j0 = voltagesource()
    pgain = potentiometer(100000, (ygain == nothing ? () : (ygain,))...)
    cdc = capacitor(1e-06)
    r1 = resistor(100000)
    r2 = resistor(1000)
    opa1 = opamp()
    r3 = resistor(33000)
    r4 = resistor(33000)
    r5 = resistor(2000)
    opa2 = opamp()
    r6 = resistor(33000)
    r7 = resistor(33000)
    r8 = resistor(33000)
    opa3 = opamp()
    r9 = resistor(33000)
    r10 = resistor(33000)
    jneg1 = voltagesource(-12.0)
    r11 = resistor(100000)
    r12 = resistor(27000)
    dcv = diode()
    opa4 = opamp()
    r13 = resistor(14000)
    jpos = voltagesource(12.0)
    jneg = voltagesource(-12.0)
    poffset = potentiometer(100000, (yoffset == nothing ? () : (yoffset,))...)
    r14 = resistor(27000)
    r15 = resistor(100000)
    jcv1 = voltagesource()
    pcv1 = potentiometer(100000, (ycv1 == nothing ? () : (ycv1,))...)
    r16 = resistor(68000)
    jcv2 = voltagesource()
    pcv2 = potentiometer(100000, (ycv2 == nothing ? () : (ycv2,))...)
    r17 = resistor(68000)
    d1 = diode()
    c1 = capacitor(1e-09)
    d2 = diode()
    c2 = capacitor(1e-09)
    d3 = diode()
    c3 = capacitor(1e-09)
    d4 = diode()
    c4 = capacitor(1e-09)
    d5 = diode()
    c5 = capacitor(1e-09)
    d6 = diode()
    opa5 = opamp()
    pr = potentiometer(100000, (yr == nothing ? () : (yr,))...)
    opa6 = opamp()
    r18 = resistor(100000)
    dp1 = diode()
    dp2 = diode()
    rtrim = resistor(20000)
    r19 = resistor(1000)
    opa7 = opamp()
    r20 = resistor(33000)
    cac = capacitor(1e-06)
    r21 = resistor(100000)
    opa8 = opamp()
    jout_1 = voltageprobe()

    # node connections
    #= node 00 =# j0[-] == pgain[1] == r2[2] == opa2["in+"] == r8[2] == jneg1[-] == r12[2] == opa4["in+"] == jpos[-] == jneg[-] == jcv1[-] == pcv1[1] == jcv2[-] == pcv2[1] == c2[2] == c3[2] == c4[2] == pr[1] == rtrim[2] == opa7["in+"] == r21[2] == jout_1[-] == opa1["out-"] == opa2["out-"] == opa3["out-"] == opa4["out-"] == opa5["out-"] == opa6["out-"] == opa7["out-"] == opa8["out-"] == gnd
    #= node 01 =# j0[+] == pgain[3]
    #= node 02 =# pgain[2] == cdc[1]
    #= node 03 =# cdc[2] == r1[1]
    #= node 04 =# r1[2] == r2[1] == opa1["in+"]
    #= node 05 =# opa1["in-"] == opa1["out+"] == r3[1] == r9[1]
    #= node 06 =# r3[2] == r4[1] == opa2["in-"] == r6[2]
    #= node 07 =# r4[2] == r5[1] == d1[+]
    #= node 08 =# r5[2] == opa2["out+"]
    #= node 09 =# r6[1] == r7[1] == r11[2] == r12[1] == dcv[-]
    #= node 10 =# r7[2] == r8[1] == opa3["in+"]
    #= node 11 =# opa3["in-"] == r9[2] == r10[1]
    #= node 12 =# opa3["out+"] == r10[2] == d6[-]
    #= node 13 =# jneg1[+] == r11[1]
    #= node 14 =# dcv[+] == opa4["out+"] == r13[1]
    #= node 15 =# opa4["in-"] == r13[2] == r15[1] == r16[2] == r17[2]
    #= node 16 =# jpos[+] == poffset[3]
    #= node 17 =# poffset[1] == r14[1]
    #= node 18 =# jneg[+] == r14[2]
    #= node 19 =# poffset[2] == r15[2]
    #= node 20 =# jcv1[+] == pcv1[3]
    #= node 21 =# pcv1[2] == r16[1]
    #= node 22 =# jcv2[+] == pcv2[3]
    #= node 23 =# pcv2[2] == r17[1]
    #= node 24 =# d1[-] == c1[1] == d2[+]
    #= node 25 =# d2[-] == c2[1] == d3[+]
    #= node 26 =# d3[-] == c3[1] == d4[+] == opa5["in+"]
    #= node 27 =# d4[-] == c4[1] == d5[+]
    #= node 28 =# d5[-] == c5[1] == d6[+]
    #= node 29 =# opa5["in-"] == opa5["out+"] == pr[3] == r19[1]
    #= node 30 =# pr[2] == opa6["in+"]
    #= node 31 =# opa6["in-"] == r18[1] == dp1[+] == dp2[-] == rtrim[1]
    #= node 32 =# c1[2] == c5[2] == opa6["out+"] == r18[2] == dp1[-] == dp2[+]
    #= node 33 =# r19[2] == opa7["in-"] == r20[1]
    #= node 34 =# opa7["out+"] == r20[2] == cac[1]
    #= node 35 =# cac[2] == r21[1] == opa8["in+"]
    #= node 36 =# opa8["in-"] == opa8["out+"] == jout_1[+]

    end
end

for yr in [0.4, 0.5, 0.6, 0.7]
    model = DiscreteModel(diode_ladder_full(Circuit, ygain=1.0, yoffset=0.5,ycv1=0.25,ycv2=0.5,yr=yr), 1//fs)
    
    ts = [n/fs for n in 0:20000]
    vin = saw(ts, freq=25)
    vcv1in = sine(ts, freq=1)
    vcv2in = [0 for t in ts]
    x = [vin[1, :], vcv1in[1, :], vcv2in]
    x = permutedims(hcat(x...))   
    y = run!(model, x)
    
    plot(ts, x[1, :], label="V_in", size=(750,350), title=string("yr = ", yr))
    plot!(ts, x[2, :], label="V_cvin1")
    display(plot!(ts, y[1, :], label="V_out"))
end

4 Conclusions

It was super easy to get up and running with the library, and it handled quite complex circuits with no issue at all! The full diode ladder circuit takes about an hour to symbolically solve (numerical solution is much faster of course), so the near instantaneous result here is impressive. The library will be a useful way to cross check my results against.