Modelling passive ladder topology circuits with Maxima

A ladder topology is essentially a cascade of voltage divider circuits; or, equivalently, a string of series impedances interspersed with shunt admittances, starting with either. We will only consider the so-called Π-topology (shunt—series—shunt—…), as the dual T-topology can be modelled by simply substituting zero for the first series impediance. (Alternatively, the voltage—current duality also means that these two topologies are mathematically the same.)

We model the circuit by finding out its ABCD parameters matrix A, which, by the definition of said matrix, is the matrix product of the matrices of individual series impedances and shunt admittances of the ladder, which may also be derived from the definition and have the following forms:

AshY (Y) := matrix ([ 1, 0 ], [ Y, 1 ]);
AseZ (Z) := matrix ([ 1, Z ], [ 0, 1 ]);

Note that the Maxima code can be extracted from this document using, for example, the following Bash command: $ pcregrep -Mao1 -- $'<pre class="maxima-code"\n>([[:print:][:cntrl:]]*?)</pre>' < 019.ladder.en.xhtml .

A ladder can then be produced from the array YZY of alternating shunt admittances and series impedances by the following function:

AladPi (YZY) := block ([ A, i ],
    i : length (YZY),
    A : matrix ([ 1, 0 ], [ 0, 1 ]),
    loop,
    A : (if oddp (i)    then AshY (YZY[i])
                        else AseZ (YZY[i])) . A,
    i : -1 + i,
    if (i > 0) then go (loop) else A);

Alternatively, we could’ve used the transfer parameters matrix T which has the same cascading property, but it neither refers to ordinarily measurable values directly nor provides a simple conversion to such. For the A matrix, the conversion to the measurable scattering parameters matrix S is as follows (where R is a 2-element array of the circuit’s characteristic impedances):

SfromA (M, R) :=  block ([ A, B, C, D, p, q ],
    A : M[1, 1], B : M[1, 2],
    C : M[2, 1], D : M[2, 2],
    p : (B - C * R[1] * R[2]),
    q : (A * R[2] - D * R[1]),
    (1 / (B + C * R[1] * R[2] + A * R[2] + D * R[1])
     *  matrix ([ p + q, 2 * R[1] * (A * D - B * C) ],
                [ 2 * R[2], p - q ])));

Let us consider a prototype Chebyshev low-pass filter of order n and ripple factor e. The coefficients of such a filter are given by:

cheblG (k, n, e) :=
    if (not integerp (k) or not integerp (n)
        or 0 > k or k > 1 + n) then false
    elseif  k > n then  if oddp (n) then 1
                        else coth (cheblb (e) / 4) ^ 2
    elseif  k > 1 then ((4 * cheblA (k - 1, n)
                         * cheblA (k, n))
                        / (cheblB (k - 1, n, e)
                           * cheblG (k - 1, n, e)))
    elseif  k = 1 then (2 / cheblg (n, e)) * cheblA (1, n)
    else 1;
cheblg (n, e)
    := sinh (cheblb (e) / (2 * n));
cheblb (e)
    := log (coth (e));
cheblA (k, n)
    := sin ((2 * k - 1) / (2 * n) * %pi);
cheblB (k, n, e)
    := cheblg (n, e) ^ 2 + sin (k / n * %pi) ^ 2;

The S matrix for a filter of order 3 and the ripple factor of 0.01 dB is then:

declare (s, complex, %omega, real, f, real);

cheblr (dB)
    := dB * log (10) / 40;
cheb_pro_YZY
    :  makelist (cheblG (i, 3, cheblr (.01)), i, 1, 3);
S_cheb_pro_3_1
    :  optimize (SfromA (AladPi ('s * cheb_pro_YZY),
                         [ 1, 1 ]))$
S_cheb_pro_3 (s)
    := subst ('s = s, S_cheb_pro_3_1);

With the actual reflection and transmission values being:

mag_dB  (v) := (20 / log (10)) * log (cabs (v));
arg_deg (v) := (180 / %pi) * carg (v);

for u : -1 thru 1 step .2 do block ([ %omega, Sm ],
    %omega : (2 ^ u),
    Sm : float (ev (S_cheb_pro_3 (%i * %omega), eval)),
    printf (true, "~5f   ~7f   ~5f   ~7f   ~5f~%", %omega,
            mag_dB (Sm[1, 1]), arg_deg (Sm[1, 1]),
            mag_dB (Sm[2, 1]), arg_deg (Sm[2, 1])))$

The characteristic frequency response shape obtained suggests the correctness of the equations used so far.

Let us obtain an expression for an arbitrary 8-element ladder topology.

S_gen_8
    :  optimize (SfromA (AladPi ([ Y[1], Z[1], Y[2], Z[2],
                                   Y[3], Z[3], Y[4], Z[4] ]),
                         [ R[1], R[2] ]))$

fortran (subst (R = R_c, S_gen_8));

As we plan to use Gnuplot’s fit command in the next section, we then edit the resulting pseudo-Fortran expression into valid Gnuplot code, as well as add the relevant supporting definitions (the specific forms for admittances and impedances used are obtained elsewhere):

j = { 0, 1 } ;
isi (a, b) = 1. / (1. / a + 1. / b) ;
array R_c[2] ;
Y_RLC (s, R, L, C) =       1. / (R + s * L) + (s * C) ;
Z_RLC (s, R, L, C) = 1. / (1. / (R + s * L) + (s * C)) ;

array v[18] ;
array Y[4] ;
array Z[4] ;

S_YZ_generic (s) = ( \
    v[1] = Z[3] * Y[4], \
    v[2] = Y[3] * (v[1] + 1), \
    v[3] = Z[2] * (v[2] + Y[4]), \
    v[4] = Y[2] * (v[3] + v[1] + 1), \
    v[5] = Z[1] * (v[4] + v[2] + Y[4]) + v[3] + v[1] + 1, \
    v[6] = R_c[2] * v[5], \
    v[7] = Y[1] * v[5] + v[4] + v[2] + Y[4], \
    v[8] = - R_c[1] * R_c[2] * v[7], \
    v[9] = Y[4] * Z[4], \
    v[10] = Z[3] * (v[9] + 1), \
    v[11] = Y[3] * (v[10] + Z[4]), \
    v[12] = Z[2] * (v[11] + v[9] + 1), \
    v[13] = Y[2] * (v[12] + v[10] + Z[4]), \
    v[14] = Z[1] * (v[13] + v[11] + v[9] + 1), \
    v[15] = v[14] + v[12] + v[10] + Z[4], \
    v[16] = Y[1] * v[15] + v[13] + v[11] + v[9] + 1, \
    v[17] = R_c[1] * v[16] , \
    v[18] = 1 / (R_c[1] * R_c[2] * v[7] + v[6] + v[17] \
                 + v[14] + v[12] + v[10] + Z[4]), \
    0) ;

Note that the Gnuplot code can be extracted from this document using, for example, the following Bash command: $ pcregrep -Mao1 -- $'<pre class="gnuplot-code"\n>([[:print:][:cntrl:]]*?)</pre>' < 019.ladder.en.xhtml | sed -e /^\<\!--/,/--\>/d .

Using the Gnuplot definitions in the previous section, we now implement the following.

S_YZ_1_1 (s) = ( \
    Y[1] = s * C1, \
    Z[1] =  isi (s * L2, 1. / (s * C2)), \
    Y[2] = s * C3, \
    Z[2] = (isi (s * L4, 1. / (s * C4)) \
            + 1. / (s * C5)), \
    Y[3] = s * C7 + 1. / (s * L6  + 1. / (s * C6)), \
    Z[3] = 1. / (s * C8), \
    Y[4] = s * C9 + 1. / (s * L10 + 1. / (s * C10)), \
    Z[4] = 1. / (s * C11), \
    0) ;

S_YZ_2_1 (s) = ( \
    S_YZ_1_1 (s), \
    Z[2] = (isi (Z_RLC (s, RL4, CL4, L4), 1. / (s * C4)) \
            + 1. / (s * C5)), \
    Y[3] = (isi (1. / Z_RLC (s, RL6, CL6, L6),  s * C6) \
            + 1. / (s * C7)), \
    0) ;

S_YZ_3_1 (s) = ( \
    S_YZ_2_1 (s), \
    Z[1] =  isi (Z_RLC (s, RL2, CL2, L2), 1. / (s * C2)), \
    Y[4] = s * C9 + isi (1. / Z_RLC (s, RL10, CL10, L10), s * C10), \
    0) ;

S_YZ_1 (s) = (S_YZ_1_1 (s), S_YZ_generic (s)) ;
S_YZ_2 (s) = (S_YZ_2_1 (s), S_YZ_generic (s)) ;
S_YZ_3 (s) = (S_YZ_3_1 (s), S_YZ_generic (s)) ;

S_YZ (s) = S_YZ_1 (s) ;

S11_YZ (s) = ( \
    S_YZ (s), \
    v[18] * (v[8] + v[6] - R_c[1] * v[16] \
             + v[14] + v[12] + v[10] + Z[4])) ;
S12_YZ (s) = ( \
    S_YZ (s), \
    2 * R_c[1] * (v[5] * v[16] - v[7] * v[15]) * v[18]) ;
S21_YZ (s) = ( \
    S_YZ (s), \
    2 * R_c[2] * v[18]) ;
S22_YZ (s) = ( \
    S_YZ (s), \
    v[18] * (v[8] - R_c[2] * v[5] + v[17] + v[14] + v[12] + v[10] + Z[4])) ;

We can now plot our example filter’s transfer function:

dB (S) = (20 / log (10)) * log (abs (S)) ;

set logscale x ;

S_YZ (s) = S_YZ_1 (s) ;
u = f_example_i (0) ;
plot [ f = f_0 / 3 : f_0 * 3 ] dB (S21_YZ (2 * pi * j * f)) ;

So far we’ve observed that the sheer number of variables that become closely intertwined in the model’s final expression make reliable retrieval of said variables from measurements largely an impossibility. A possible alternative approach is to model the transfer function as a ratio of two polynomials and attempt to retrieve poles and zeros of the function instead. An obvious problem there is that while zeros correspond directly to the network’s resonant RLC circuits, the interpretation of deviation of poles from their desired location is by no means as straightforward.

Following the principle of knowing one’s tools (that is, taking deep understanding of the operation of the software used as an imperative, in place of allowing for mindless application of said software until some result is obtained) and noting that Maxima may be an overkill, investigating the applicability of the smaller Jacal system for the task at hand appears promising.

For completeness, here are the Maxima expressions for modelling elliptic 5-th order bandpass filters discussed above.

ilsi (L) := 1 / lsum (1 / x, x, L);

S_ell_8_1_1 : optimize (subst ([
    Y[1] = s * C1,
    Z[1] =  ilsi ([ s * L2, 1 / (s * C2) ]),
    Y[2] = s * C3,
    Z[2] = (ilsi ([ s * L4, 1 / (s * C4) ])
            + 1 / (s * C5)),
    Y[3] = s * C7 + 1 / (s * L6  + 1 / (s * C6)),
    Z[3] = 1 / (s * C8),
    Y[4] = s * C9 + 1 / (s * L10 + 1 / (s * C10)),
    Z[4] = 1 / (s * C11) ], S_gen_8))$

[Chebyshev]

Chebyshev filter // Wikipedia. — URI: http://en.wikipedia.org/wiki/Chebyshev_filter

[Elliptic]

Elliptic filter // Wikipedia. — URI: http://en.wikipedia.org/wiki/Elliptic_filter

[Gnuplot]

Gnuplot. — URI: http://gnuplot.info/

[Jacal]

The JACAL symbolic math system. — URI: http://people.csail.mit.edu/jaffer/JACAL

[Maxima]

Maxima, a Computer Algebra System. — URI: http://maxima.sourceforge.net/

[N-port]

N-port matrix conversions // Qucs technical papers. — URI: http://web.archive.org/web/20170129063841/http://qucs.sourceforge.net/tech/node98.html

[Rational]

Elliptic rational function // Wikipedia. — URI: http://en.wikipedia.org/wiki/Elliptic_rational_function

[Scattering]

Scattering parameters // Wikipedia. — URI: http://en.wikipedia.org/wiki/Elliptic_filter