# Chebyshev bandstop LC ladder design in Python

README.md

Chebyshev bandstop LC ladder design in Python

Simple script used to experiment with a design of a "FM trap" filter.

Install

python3 -m venv env
source env/bin/activate
pip3 install -r requirements.txt

bandstop_des.py

import math

import numpy as np
import matplotlib.pyplot as plt

from lcapy import Series, Shunt, L, C, R

import PySpice.Logging.Logging as Logging

logger = Logging.setup_logging()

from PySpice.Plot.BodeDiagram import bode_diagram_gain
from PySpice.Spice.Netlist import Circuit
from PySpice.Spice.BasicElement import Capacitor, Inductor
from PySpice.Unit import *

_DEBUG_ = False


def debug(msg):
    if _DEBUG_:
        print(msg)


# fmt: off
# R ratio = 1
cheb_1db_ripple_tab = {
    3: [2.023593, 0.994102, 2.023593],
    5: [2.134882, 1.091107, 3.000923, 1.091107, 2.134882],
    7: [2.166557, 1.111509, 3.093642, 1.173521, 3.093642, 1.111509, 2.166557],
    9: [2.179723, 1.119177, 3.121434, 1.189673, 3.174634, 1.189673, 3.121434, 1.119177, 2.179723],
    11: [2.186417, 1.122901, 3.133679, 1.195713, 3.197903, 1.204134, 3.197903, 1.195713, 3.133679, 1.122901, 2.186417],
    }

cheb_01db_ripple_tab = {
    3: [1.03156, 1.147397, 1.03156],
    5: [1.146813, 1.371213, 1.975003, 1.371213, 1.146813],
    7: [1.181178, 1.422806, 2.096671, 1.573401, 2.096671, 1.422806, 1.181178],
    9: [1.195672, 1.442601, 2.134551, 1.616723, 2.205369, 1.616723, 2.134551, 1.442601, 1.195672],
    11: [1.203089, 1.452293, 2.151438, 1.633242, 2.237794, 1.655924, 2.237794, 1.633242, 2.151438, 1.452293, 1.203089],
}
# fmt: on

# These are values available at the time
standardized_capacitance_values = [2.2, 12, 15, 33, 47, 100]


def norm_cap(c):
    c *= 1e12  # in pF
    cap_val = [
        c / 2 for c in standardized_capacitance_values
    ]  # admit series connections
    cap_val += [
        c + d
        for c in standardized_capacitance_values
        for d in standardized_capacitance_values
    ]  # admit parallel connections
    return 1e-12 * round(cap_val[np.argmin([abs(cs - c) for cs in cap_val])], 2)


def get_circuit(lc_ser, lc_par, Rs=50, Rp=50):
    circuit = Circuit("Cheb bandstop")
    circuit.SinusoidalVoltageSource("input", "inp", circuit.gnd, amplitude=1 @ u_V)

    circuit.R("s", "inp", 1, Rs)
    circuit.R("p", "outp", circuit.gnd, Rp)

    for i, (l, c) in enumerate(lc_par):
        p1 = i + 1
        p2 = "outp" if i + 1 == len(lc_par) else i + 2

        circuit.L(f"{i+1}p", p1, p2, (l * 1e9) @ u_nH)
        circuit.C(f"{i+1}p", p1, p2, (c * 1e12) @ u_pF)

    for i, (l, c) in enumerate(lc_ser):
        p1 = "outp" if i + 1 == len(lc_ser) else i + 1

        circuit.C(f"{i+1}s", p1, f"{i + 1}s", (c * 1e12) @ u_pF)
        circuit.L(f"{i+1}s", f"{i + 1}s", circuit.gnd, (l * 1e9) @ u_nH)

    return circuit


def pretty_circuit(circuit):
    lines = []
    for e in circuit.elements:
        if isinstance(e, Capacitor):
            v = round(e.capacitance.value, 2)
            u = e.capacitance.prefixed_unit
            lines.append(f"{e.name} = {v} {u}")
        elif isinstance(e, Inductor):
            v = round(e.inductance.value, 2)
            u = e.inductance.prefixed_unit
            lines.append(f"{e.name} = {v} {u}")

    return "\n".join(lines)


def get_sym_network(n):
    if n % 2 == 0:
        raise ValueError("Filter order must be odd")

    net = Shunt(C("C1s") + L("L1s"))
    for i in range(n - 1):
        if i % 2 == 0:
            net = net.chain(Series(C(f"C{i // 2 + 1}p") | L(f"L{i // 2 + 1}p")))
        else:
            net = net.chain(Shunt(C(f"C{i // 2 + 2}s") + L(f"L{i // 2 + 2}s")))

    return net


def cheb_bandstop_des(A_min, f_pl, f_ph, f_sl, f_sh, A_max=1, k_m=50):

    w_pl, w_ph, w_sl, w_sh = list(
        map(lambda f: 2 * math.pi * f, [f_pl, f_ph, f_sl, f_sh])
    )

    debug("--------------------------------")
    debug(f"w_pl = {w_pl} rad/s")
    debug(f"w_ph = {w_ph} rad/s")

    debug(f"w_sl = {w_sl} rad/s")
    debug(f"w_sh = {w_sh} rad/s")

    w_bw = w_ph - w_pl
    w_center = math.sqrt(w_ph * w_pl)

    debug(f"w_bw = {w_bw} rad/s ({w_bw / (2*math.pi)} Hz)")
    debug(f"w_center = {w_center} rad/s ({w_center / (2*math.pi)} Hz)")

    omega_p = w_bw / (w_ph - w_pl)
    omega_s = w_bw / (w_sh - w_sl)

    debug(f"omega_p = {omega_p}")
    debug(f"omega_s = {omega_s}")

    epsilon_p = math.sqrt(math.pow(10, A_max / 10) - 1)
    epsilon_s = math.sqrt(math.pow(10, A_min / 10) - 1)

    debug(f"epsilon_p = {epsilon_p}")
    debug(f"epsilon_s = {epsilon_s}")

    n = math.ceil(math.acosh(epsilon_s / epsilon_p) / math.acosh(omega_s / omega_p))
    if n % 2 == 0:
        n += 1

    print(f"n = {n}")

    if not n in cheb_1db_ripple_tab.keys():
        raise ValueError(f"Unsupported filter order: {n}")

    k_fbs = w_bw

    if A_max == 1:
        norm_values = cheb_1db_ripple_tab[n]
    elif A_max == 0.1:
        norm_values = cheb_01db_ripple_tab[n]
    else:
        raise ValueError("Unsupported A_max")

    ser_v = [norm_values[i] for i in range(0, len(norm_values), 2)]
    par_v = [norm_values[i] for i in range(1, len(norm_values), 2)]

    def ser_compute(v):
        C = (k_fbs * v) / (k_m * (w_center**2))
        L = k_m / (k_fbs * v)
        return L, C

    def par_compute(v):
        C = 1 / (k_m * k_fbs * v)
        L = (k_m * k_fbs * v) / (w_center**2)
        return L, C

    LC_ser = list(map(ser_compute, ser_v))
    debug(f"LC_serial = {LC_ser}")

    LC_par = list(map(par_compute, par_v))
    debug(f"LC_parallel = {LC_par}")

    return n, LC_ser, LC_par


f_pl, f_ph, f_sl, f_sh = [80e6, 116e6, 88e6, 108e6]  # MHz
n, lc_ser, lc_par = cheb_bandstop_des(45, f_pl, f_ph, f_sl, f_sh, 0.1)

circuit_ideal = get_circuit(lc_ser, lc_par, Rp=75)
pc = pretty_circuit(circuit_ideal)
print("Ideal filter values")
print(pc)

lc_ser = [(1e-9 * round(l * 1e9), norm_cap(c)) for l, c in lc_ser]
lc_par = [(1e-9 * round(l * 1e9), norm_cap(c)) for l, c in lc_par]

circuit = get_circuit(lc_ser, lc_par, Rp=75)
pc = pretty_circuit(circuit)
print("Rounded filter values")
print(pc)

net = get_sym_network(n)
net.cct.draw("filter_circuit.png", label_nodes=False, draw_nodes="connections")

net = Series(R("Rs")).chain(net.chain(Shunt(R("Rp"))))
s_circuit = net.cct

s_circuit = s_circuit.subs("Rp", 75)
s_circuit = s_circuit.subs("Rs", 50)

for i, (l, c) in enumerate(lc_ser):
    s_circuit = s_circuit.subs(f"L{i + 1}s", l)
    s_circuit = s_circuit.subs(f"C{i + 1}s", c)

for i, (l, c) in enumerate(lc_par):
    s_circuit = s_circuit.subs(f"L{i + 1}p", l)
    s_circuit = s_circuit.subs(f"C{i + 1}p", c)

s_circuit.draw(
    "filter_circuit_with_values.png",
    label_nodes=False,
    draw_nodes="connections",
    node_spacing=4.0,
)

figure, ax = plt.subplots(figsize=(20, 10))

simulator = circuit_ideal.simulator(temperature=25, nominal_temperature=25)
analysis = simulator.ac(
    start_frequency=f_pl * 0.5 @ u_Hz,
    stop_frequency=f_ph * 1.5 @ u_Hz,
    number_of_points=1000,
    variation="dec",
)

plt.title("Bode Diagram Gain")
bode_diagram_gain(
    axe=ax,
    frequency=analysis.frequency,
    gain=20 * np.log10(np.absolute(analysis.outp)),
    marker=".",
    color="blue",
    linestyle="-",
    label="Ideal",
)

simulator = circuit.simulator(temperature=25, nominal_temperature=25)
analysis = simulator.ac(
    start_frequency=f_pl * 0.5 @ u_Hz,
    stop_frequency=f_ph * 1.5 @ u_Hz,
    number_of_points=1000,
    variation="dec",
)

bode_diagram_gain(
    axe=ax,
    frequency=analysis.frequency,
    gain=20 * np.log10(np.absolute(analysis.outp)),
    marker=".",
    color="red",
    linestyle="-",
    label="Rounded values",
)

plt.text(0.05, 0.8, pc, fontsize="20", ha="left", va="top", transform=ax.transAxes)
plt.tight_layout()
plt.legend(fontsize="20", loc="lower right")

# H = s_circuit.transfer(0, 1, 2, 3)
# H.bode_plot((f_pl * 0.5, f_ph * 1.5), axes=ax, dbmin=None)

plt.show()

filter_circuit.png

filter_circuit_with_values.png

gain.png

lp_cheb_elements.py

import math


def cheb_lp_elements(n, A_p, R_l):
    if n < 0:
        raise ValueError(f"Order 'n' must be greater than zero")

    refl = math.sqrt(math.pow(10, A_p / 10) - 1)
    if n % 2 == 0:
        a = (1 + (refl**2)) * 4 * R_l / ((R_l + 1) ** 2)
        if a > 1:
            raise ValueError(f"Specification not implementable")
    else:
        a = 4 * R_l / ((R_l + 1) ** 2)

    x = math.pow((1 / refl) + math.sqrt((1 / (refl**2)) + 1), 1 / n)
    y = math.pow(
        math.sqrt((1 - a) / (refl**2)) + math.sqrt(((1 - a) / (refl**2)) + 1), 1 / n
    )

    x1 = x - 1 / x
    y1 = y - 1 / y
    num1 = math.pi / (2 * n)

    comp = (4 * math.sin(num1)) / (x1 - y1)

    elems_odd = []
    elems_even = []

    elems_odd.append(comp)

    for m in range(1, 1 + n // 2):
        b = (
            x1**2
            - 2 * x1 * y1 * math.cos((4 * m - 2) * num1)
            + (y1**2)
            + (2 * math.sin((4 * m - 2) * num1)) ** 2
        )
        comp = (16 * math.sin((4 * m - 3) * num1) * math.sin((4 * m - 1) * num1)) / (
            b * comp
        )

        elems_even.append(comp)

        if n <= 2 * m:
            break

        b = (
            x1**2
            - 2 * x1 * y1 * math.cos(4 * m * num1)
            + y1**2
            + (2 * math.sin((4 * m * num1))) ** 2
        )
        comp = (16 * math.sin((4 * m - 1) * num1) * math.sin((4 * m + 1) * num1)) / (
            b * comp
        )

        elems_odd.append(comp)

    elems = []
    while True:
        if len(elems_odd) == 0:
            break
        else:
            elems.append(elems_odd.pop(0))
        if len(elems_even) == 0:
            break
        else:
            elems.append(elems_even.pop(0))

    return elems


for n in range(1, 12):
    try:
        es = cheb_lp_elements(n, 1, 1)
        es = [round(e, 6) for e in es]
        print(f"{n}: {es},")
    except ValueError:
        pass

requirements.txt

asttokens==2.4.1
certifi==2024.2.2
cffi==1.16.0
charset-normalizer==3.3.2
contourpy==1.2.1
cycler==0.12.1
decorator==5.1.1
exceptiongroup==1.2.1
executing==2.0.1
fonttools==4.53.0
idna==3.7
ipython==8.25.0
jedi==0.19.1
kiwisolver==1.4.5
lcapy==1.22
matplotlib==3.2.0
matplotlib-inline==0.1.7
mpmath==1.3.0
networkx==3.3
numpy==1.26.4
packaging==24.0
parso==0.8.4
pexpect==4.9.0
pillow==10.3.0
ply==3.11
prompt_toolkit==3.0.45
property-cached==1.6.4
ptyprocess==0.7.0
pure-eval==0.2.2
pycparser==2.22
Pygments==2.18.0
pyparsing==3.1.2
PySpice==1.5
python-dateutil==2.9.0.post0
PyYAML==6.0.1
requests==2.32.3
scipy==1.13.1
six==1.16.0
stack-data==0.6.3
sympy==1.12.1
traitlets==5.14.3
typing_extensions==4.12.1
urllib3==2.2.1
wcwidth==0.2.13