Category: General

For variety, I switched from SageMath to Sympy since Sympy is a small more focused package. However, I’m struggling with it’s lack of seeming lack of equational reasoning. When you define an equation via Eq(lhs,rhs), it’s not a “first class” object – ie. you can’t divide it by 2 and expect both sides to be half as big, and you have to manually say “subs(titute) eq.lhs with eq.rhs” each time. Finally, I haven’t found good ways to say “substitute just this one occurrence of foo” – sympy wants to do them all. This is all totally anathema to how I think about maths – I’m used to pointing at some subexpression and saying “and use equation 2 to simplify that” etc.

To try out Sympy, I took up the challenge in Chapter 1 of Nahin’s Transients or Electrical Engineers book, in which he conjures up a second-order differential equation for a simple R/L/C ciruit and invites the reader to check he’s not lying. This was a decent workout, since 1) it involves first- and second-derivatives, 2) it involves plenty of equations and manipulating expressions. To be honest, I found this hard. I tried doing it on computer and got nowhere. I then did it by hand using paper and pen to figure out the steps. Then I went back to the computer and had to tackle several “how do I get the computer to do x” puzzles before I got anywhere.

In the end, I did indeed verify that Nahin is telling the truth (phew!) but it was hours and hour of headscratching work. It’s nice to have the computer verify I haven’t err’d in my derivation, but this was not “computer makes task easier” stuff – quite the opposite! I guess the nice bit comes later when I can get the computer to actually solve the differential equations for concrete cases. But this ‘abstract algebra manipulation’ stuff is painful.

In the end, the “proof” looks like..

from sympy import *
from sympy.abc import R,L,C,t

i = Function("i")(t)
i1 = Function("i1")(t)
i2 = Function("i2")(t)
v = Function("v")(t)
u = Function("u")(t)

# Chapter 1 of Transients For Electrical Engineers

# Equations which describe the given "R ser L par (R ser C)" circuit (ie. KVL/KCL)
eq1 = Eq(i, i1 + i2 )
eq2 = Eq(u, i*R + v)
eq3 = Eq(v, L * i1.diff())
eq4 = Eq(v, i2*R + (1/C) * Integral(i2,(t,0,t)))

# The book states a large differential equation with 'input' voltage u(t) + derivative on left,
# and current i(t) (+derivatives) on right
(u.diff(t,2) + (R/L) * u.diff() + 1/(L*C)*u)

# and challenges you to check the equality does actually follow from above 4 equations

# First step is easy; there's only one equation involving 'u' so let's apply that everywhere
_.subs( eq2.lhs, eq2.rhs).simplify()

# Now we have to be tactical; we have v and v' and v'' but we want to do different things to each
# I don't know how to tell sympy to "just substitutes the first instance, not all instances"
# So start by substituting for v'' (since that leaves v' and v alone for now) and use eq4 the capacitor eqn
_.subs(eq4.lhs.diff(t,2), eq4.rhs.diff(t,2)).simplify()

# And now we can tackle v' and then v - both use eq3 (the inductor one)
# This leaves us with just currents, albeit a mixture of i/i1/i2
_.subs(eq3.lhs.diff(t), eq3.rhs.diff(t)).simplify()
_.subs(eq3.lhs,eq3.rhs).expand().simplify()

# Now we want to group + collapse the i1's and i2's using eq1.
# But I don't know how to get sympy to "see" that we have lots of (i1+i2) patterns 
# So as a hack, we can always subtract something then add something equal to it and preserve equality, right?
# So we use that to turn i1+i2 into i ...
(_ + diff(eq1.rewrite(Add),t)/C).expand()

# Getting there .. now we need to sweep up the twice-differentiated cases..
(_ + R*eq1.rewrite(Add).diff(t,2)).expand()
# w^5 - which was what we wanted

Using Sagemath for electronics is nice because it forces you to be explicit about your assumptions and calculations. Ideally, you’d start by telling Sagemath about Maxwell’s laws and build up from there. For DC conditions, like calculating the bias for a common emitter amplifier we can assume no changing magnetic fields, so Faraday’s law of induction simplifies to Kirchhoff’s Voltage Law. We also need Ohm’s Law for the voltage-current relationship in resistors. I supplied the series/parallel resistance formula (but those too could come from KVL/KCL + Ohms law).

First let’s tell Sagemath about Kirchhoff’s Voltage Law:

from functools import reduce
from operator import add, neg

# Teach sagemath about physics (Faraday's law of induction in its KVL guise)

# Kirchhoff voltage law says sum of voltage drops around a loop is zero.
# Note: this is really just Faraday's law of inducion, for the special case of no changing magnetic fields
# Note: this is taking in symbolic expressions, and returning a symbolic expression
def KVL( voltages ): 
    return reduce(add,voltages) == 0

# Convenience case for where we have one voltage source, then we give remaining voltage drops as positive values 
# ie. KVL_vsource(5, [1,2,x] means a 5v source, then drop of 1v then drop of 2v then drop of x volts.
def KVL_vsource( vsource, vothers ):
     return KVL( [vsource] + list(map(neg,vothers)))

Next, let’s tell it about series/parallel resistors and Thevenin equivalent for voltage dividers (again, would be nice to derive this from KVL/KCL/Ohm)

# Kirchhoff's Current Law (aka conservation of charge) and Ohms Law give composition rules for impedance
# TODO: derive directly from KCL/Ohm.
recip(x) = 1 / x
ipar(a,b) = recip( recip(a) + recip(b))
iser(a,b) = a + b

# Thevenin equivalent voltage/resistance for divider (r1 and r2)
Vth_divider(v, r1, r2) = v * r2 / (r1+r2)
Rth_divider(r1, r2) = ipar(r1,r2)

Now we can pick some values for a common-emitter amplifier

# Transistor params
beta = 100
Vbe = 0.7

# Common emitter amplifier, biased using divider, with emitter feedback
Vcc = 5
Rtop,Rbot = 10e3, 3.3e3
Re = 470
Rc = 2e3

And calculate the quiescent voltages + currents:

# Turn the divider into thevenin equivalent
Vth = Vth_divider(Vcc, Rtop, Rbot )
Rth = Rth_divider(Rtop, Rbot)

# To get Ib, we use Kirchhoff's Voltage Law for loop around base/emitter, and solve for Ib
Ib = var("Ib")
eqn = KVL_vsource( Vth, [ Ib*Rth, Vbe, beta*Ib*Re])
print("KVL equation is ", eqn)
Ib = solve(eqn, Ib, solution_dict=True)[0][Ib].n()
print("So Ib is", Ib*1000, "mA")
# Finally, get voltage drop over Rc to get output voltage.
Ic = beta * Ib

Vc = Vcc - Ic * Rc
print("Vc is ", Vc)

var("Vce")
vce_eqn = KVL_vsource(Vcc, [Ic*Rc, Vce, Ic*Re])
Vce = solve(vce_eqn, Vce, solution_dict=True)[0][Vce].n()
print("Vce is ", Vce)

Here we’re just using Sagemath as a numeric calculator. It’s kinda neat to package up Kirchhoff’s Voltage Law so we can just state all the knowns and let Sagemath find the unknowns, without us having to do algebra. But fundamentally, this is stuff we could easily do in excel/python/etc.

But with Sagemath we can stick with symbolic expression. Let’s start by keeping the potential divider as-is, but we’ll make Rc and Re symbols.

# Common emitter amplifier, biased using divider, with emitter feedback
Vcc = 5
Rtop = 10e3
Rbot = 3.3e3
Re = var("Re")
Rc = var("Rc")

# Turn the divider into thevenin equivalent
# Note: some textbooks just ignore supply impedance, ie. assume it's a perfect voltage source which has fixed voltage
# regardless of how much current is drawn; this isn't the case for resistor-divider.  We can either make heuristic assumptions
# (eg. assume load has "high enough" impedance, so draws little current) or we can calculate exactly via Thevenin equivalent.
Vth = Vth_divider(Vcc, Rtop, Rbot )
Rth = Rth_divider(Rtop, Rbot)

# To get Ib, we use Kirchhoff's Voltage Law for loop around base/emitter, and solve for Ib
Ib = var("Ib")
eqn = KVL_vsource( Vth, [ Ib*Rth, Vbe, beta*Ib*Re])
Ib = solve(eqn, Ib, solution_dict=True)[0][Ib]
Ic = beta * Ib

Vc = Vcc - Ic * Rc
print("Vc is ", Vc)
solve([Vc == Vcc/2, Ic*Re==Vcc/10], [Rc,Re], solution_dict=True)

# Output:
Vc is  -719/10*Rc/(133*Re + 3300) + 5

So now the collector voltage is left as an expression, dependent on the (as yet unchosen values of Rc and Re).

Rather than picking them manually, we can express constraints:

solve([Vc == Vcc/2, Ic*Re==Vcc/10], [Rc,Re], solution_dict=True

# Output
[{Rc: 13750/9, Re: 2750/9}]

So we’ve effectively just stated the design rules of “collector should be midrail” and “Ve should be about 10% of supply”, and Sagemath has found the values of Rc and Re which “make it so”. Nice.

How about if we wanted IC=20mA and Vc at midrail?

solve([Vc == Vcc/2,  Ic==0.01], [Rc,Re], solution_dict=True)

# Output
[{Rc: 250, Re: 3890/133}]

This time we used Sagemath to model an LC bandpass filter, centered on 100MHz, and calculate the insertion loss in dB at various frequencies. It’s based on what the filter calculator defaults to for a 1st order shunt-first Chebyshev bandpass filter.

# Complex impedance for capacitor and inductor
C(c,f) = -I / (2 * pi * f * c)
L(l,f) = I * 2 * pi * f * l

# Composition rules for impedance
recip(x) = 1 / x
ipar2(a,b) = recip( recip(a) + recip(b))
ipar3(a,b,c) = recip( recip(a) + recip(b) + recip(c))
iser(a,b) = a + b

# A simple RC lowpass filter
V = 10
Rs = 50 # source impedance
Rl = 50 # load impedance

# These are just the default values from https://rf-tools.com/lc-filter/
# They give a bandpass filter centered at 100MHz, with 20MHz bandwidth.
C1 = 48.58e-12
L1 = 52.67e-9

# This is the overall LC bandpass filter, with source impedance and load.
impedance = iser( Rs, ipar3(C(C1,f),L(L1,f),Rl) )

# The filter website shows insertion loss, which looks at how the power in the load
# is affected by putting in the LC filter.

power(v,r) = (v^2 / r).abs()

# Without our filter ("device under test"), the power is simple - half the voltage
# is dropped over the source impedance, leaving half for the load.
power_without_dut = power(V/2,Rl)

# To calculate power in load with our filter, we need to know the current
current = V / impedance
V_over_load = V - (current * Rs)
power_in_load = power(V_over_load, Rl)

as_dB(frac) = log(frac,10) * 10
insertion_loss_dB = as_dB(power_in_load / power_without_dut)

insertion_loss_dB(51750194).n()  # ans: -3.29183565814456
insertion_loss_dB(41896869).n()  # ans: -5.05202056128256

These values for insertion loss match what the calculator website shows in the “S parameters” tab.

How to think about this? If you take away the LC bit, you’re just left with a “maximum power theorem” situation, ie. source impedance and load impedance are both constant and equal.

Now, let’s introduce the LC again but imagine that there’s some magic circumstance in which it’s impedance ends up zero – meaning we’ll still be in “maximum power theorem” situation. Since the L and C are of “opposite nature”, there’ll be some frequency where the L has reactance-blah reactance at the same time as the C has negative-blah impedance. With the L and C in parallel we sum their impedances (actually, reciprocal of sum of reciprocals) so they cancel out.

In all other circumstances, either the L will “win” with a larger reactance or the C will “win” and we’ll net out to the L-par-C having some non-zero impedance. This will pull us away from our perfect 50ohm “maximum power” point and the power in the load will be less.

At low frequencies, the inductor (L) will present a low-impedance path, meaning that the overall impedance will be much lower than 50. This results in a higher current, a greater voltage drop over the source impedance leaving our parallel RLC bundle seeing a lower voltage.

At high frequencies, it’s the capacitor which presents a low-impedance path, and by the same argument, we get a lower voltage over our load.

But when we’re dealing with RF we’re dealing with waves, and when there’s a change in impedance you get reflections. If our bandpass filter only had 50ohms impedance at it’s centre frequency, then it has a non-50ohm impedance everywhere else – which causes reflections. Actually, this is kinda obvious from conservation of energy: if there’s energy coming in, but it’s not ending up in the load then either it’s absorbed in the filter itself or it must be reflected. If our filter was absorbing energy, it would be heating up (unlikely, since we build filters with low-resistance components), or it would have to be radiating (unlikely unless the filter was physically very big).

If we’re using a bandpass filter at the input of a receiver, all the signals outside the passband are getting bounced back towards the antenna (the amount of energy here is tiny so it’s not a problem).

If we use a lowpass filter to “get rid of” harmonics at the output of a transmitter, the energy of those harmonics is heading straight back towards the power amp.

If we use an IF bandpass filter in a superhet receiver, then those reflected signals will be heading back towards a mixer where they can re-mix and cause problems.

If we’re worried by those reflections, what can we do? The answer is to use a diplexer – essentially a fork in the road leading to a lowpass and highpass filters respectively. Whichever part (low or high frequencies) we don’t want can be sent to a 50ohm terminator to be absorbed and turned into heat – and overall the whole thing looks like 50ohms across all frequencies This is something I’ve seen used in the “High Performance Direct Conversion Receivers” ARRL article. Also in this “Popcorn DC Receiver”, this one and this one.

When learning about amateur radio, there’s loads of discussion about reflections that happen when the impedance of your antenna isn’t 50ohm at the frequency you’re sending on. This leads to reflections back up your 50ohm coax, and the sum of the outgoing and reflected wave leads to standing waves. The standing waves can be characterised by the “standing wave ratio” (SWR) which is a measurable quantity, and can therefore be used to indicate that your antenna impedance isn’t right. However, if your feeder cable is attenuating the waves, the reflected wave might be weak by the time it gets back to the home end of the cable – which could wrongly lead you to infer that you have the impedance spot on.

Despite all this chatter about reflections from the antenna, I have seen almost no mention of reflections from filtering steps. I guess it’s just a question of degree. Reflections from antennas are bad because its potentially the full transmit power at your transmit frequency that’s bouncing back. Your lowpass filter is also after the power amp, but in normal operation the reflected power is much small since it’s only the power from the harmonics that bounces back. If your power amp was very very nonlinear you might have more power in your harmonics. Or if you had been sending at 7MHz with a 8MHz lowpass filter, and then switched up to sending at 14MHz but didn’t change to an appropriate lowpass filter you’d end up getting all your transmit power bouncing back.

I entered this particular rabbithole with a desire to understand how the various kinds of passive filters work. I’ve already physically constructed bandpass and lowpass filters, but I had to get the values for my capacitors and hand-wound inductors from a mysterious but helpful online filter calculator. Being a cynical kinda person, I first “built” each filter using ngspice (the circuit simulation tool) to check they looked sane before heating up my soldering iron.

But that online filter calculator is kinda mysterious. It lets you build various kinds of filters – Butterworth, Chebyshev, Elliptic – and each one comes in different flavours, such as shunt-first or series-first.

They’re all made from combinations of capacitors and inductors – the two standard frequency-dependent building blocks, which harness electric and magnetic fields respectively. I know how a capacitor behaves on its own, and I know how an inductor behaves on its own – but I don’t know what magic happens when you put a bunch of them together.

Buoyed by previous suggest in modelling pi-attenuators in Sagemath from first principles, I decided to go on a journey to also model various families of passive filters in Sagemath too. The hope is that I “tell” Sagemath about the axioms – what the impedance of each component is, and how series/parallel combinations work – along with the combinations of capacitors/inductors that make up a filter, and out should pop the shapes of whatever bandpass filter you like Of course, this is exactly what the SPICE simulator will do. But the attractive thing about Sagemath is that it’s tracking the entire equation relating everything, and hence you can ask many more questions of it.

I initially tried to dive straight in an model a 3rd order Chebyshev bandpass filter. I got kinda close, but couldn’t get to the point of having correct bandpass filter behaviour. So I took that as a hint that perhaps I should be starting with something simpler.

The simplest frequency-dependent circuit is an RC filter. I’ll use R=4700ohm and C=47nF. The resistor (with constant impedance) and capacitor (with impedance that decreases as frequency rises) act as a potential divider.

Now, a RC filter is simple enough that you can work it out by hand directly. Since I use emacs a lot, I often use emacs lisp as a kind of quick calculator, so here’s how I calculate the impedance (xc) and the output voltage (assuming 10v sine input):

(defun 1/ (x) (/ 1 x))
(defun sq (x) (* x x))

(let* ((f 10e3)
       (c 47e-9)
       (r 4700)
       (xc (1/ (* 2 pi f c))))
       (* 10 (/ xc (sqrt (+ (sq r) (sq xc))))))

;; Give 0.719, meaning at 10kHz the output is 0.719 volts.

I can also use ngspice to simulate the circuit, and do a sweep across frequencies:

RC filter

Vin 1 0 DC 0 AC 10
R1 1 2 4700
C1 2 0 47e-9

.control
AC LIN 10000 1 20000
run
set hcopypscolor
hardcopy rc.ps vm(2) xlimit 0 11k
shell evince rc.ps
.endc
.end

This plot visually confirms that at 10kHz, voltages is around 0.7v.

Now for the sagemath attempt. Sagemath knows about complex numbers, so we can represent impedances directly as complex numbers (sagemath uses uppercase I):

# Complex impedance for capacitor and inductor
C(c,f) = -I / (2 * pi * f * c)
L(l,f) = I * 2 * pi * f * l

# Composition rules for impedance
recip(x) = 1 / x
ipar(a,b) = recip( recip(a) + recip(b))
iser(a,b) = a + b

# If you have [v] volts across  a potential divider 
# consisting of impedance [a] on top, [b] below, 
# connected to [load], what voltage does the load see?
pd_vout(v, a, b, load) = v * ipar(b,load) / iser(a,ipar(b,load))

# Our example RC lowpass filter
R1 = 4700
C1 = 47e-9
overall = pd_vout(10, R1, C(C1,f), 10e9)

# What's the voltage at 10kHz?  Note: use n() to get 
# a number instead of a maths expression, and abs() to
# get the complex magnitude
overall(10e3).n().abs()

# Answer: 0.718621364643

It gets the right answer, but what’s interesting (to me) is that it’s clearly there’s no additional magic going on. Spice might be doing all sorts of elaborate things to get it’s answer. But in our Sagemath version, it’s obvious that the only ingredients are 1) impedance, represented as a complex number, 2) combinations of impedances, and 3) ohms law, as expressed in the potential divider equation.

In the next post, I’ll move onto something rather more useful – an LC bandpass filter.

Attenuators are one of those “needed in practise, but not in theory” devices. For example, my signal generator has a frequency/counter mode – but the manual warns that the maximum safe voltage is 5V. If my “device under test” is outputting a 9v signal, then I need to bring that down by half to measure its frequency safely.

This sounds simple – potential dividers are commonly used in circuits for exactly this purpose.

A simple potential divider – let’s say consisting of two 100 ohm resistors – will half the voltage … so long as you don’t connect any load. Once you connect a load (say a 1k resistor), it appears in parallel to the lower resistor turning your 100ohm resistor into a 91ohm resistor – and now you only get 47% of the initial voltage instead of 50%. This “sagging” effect can be reduced by using smaller resistors in the potential divider, but at the cost of a larger current flowing through them, wasting lots of energy.

But when dealing with radio frequencies, we have an additional concern: we need to ensure that we stick with a 50 ohm impedance everywhere – otherwise we’ll cause reflections. This means that 1) to our upstream, we need to look like a 50ohm load (assuming our downstream load is also 50ohm), and 2) to our downstream load we need to look like a 50ohm source (assuming our upstream is 50ohm).

The simple “potential divider” doesn’t work like that. If our load is 50ohm, then a divider that uses 100ohm resistors will appear as 133ohms overall – no good!

Is there some ‘magic’ value for the resistors that’ll make everything look like a nice 50ohms? Here’s some sagemath code to calculate the resistor value you’d need …

# Resistor equations
rser(a,b) = a + b
rpar(a,b) = 1/ ((1/a) + (1/b))
potentialDivider(r,load) = rser(r,rpar(r,load))

r = var('r')
result = solve(potentialDivider(r,50) == 50,r,solution_dict=True)
result[1][r].n()

    30.9016994374947

So a potential divider made from 30.9ohm resistors would make the source think we’re a well-behaved 50ohm load.

But what would the load see? It’ll see the source impedance in series with one 30.9 resistor, all in parallel with the other 30.9ohm resistor – namely 22.3ohms. No good!

So a simple potential divider, using 2 resistors, isn’t working out. We need to go up to a 3-component approach.

Textbooks will tell us that we can solve this problem with a pi network; namely a series resistor with two equal-sized shunt resistor before and after it. They’ll also give us a bunch of equations for finding the resistor values for various attenuations.

But wouldn’t it be fun if we could just describe the various contraints to a computer, and have it figure out all the equations for us?

General pi attenuator from first principles

First we have to explain to sagemath what a pi-network looks like: namely, how to calculate the impedance that the source and load will see. This allows us to express the “source/load impedance must match” criteria.

Then we need a couple of helpers to let figure out what the power through the load will be. This allows us to express the attenuation (ie. the ratio between source power and load power).

# Impedance of a general pi-network (a=shunt, b=series, c=shunt), as seen from the source
rpi_src(zload,a,b,c) = rpar(a,rser(b,rpar(zload,c)))

# Impedance of a general pi network, as seen from the load
rpi_load(zsrc,a,b,c) = rpar(c,rser(b,rpar(zsrc,a)))

# If you have [v] volts across  a potential divider consisting of resistor [a] on top, [b] below, connected to [load], what voltage does the load see?
pd_vout(v, a, b, load) = v * rpar(b,load) / rser(a,rpar(b,load))

# What's the power through a resistor [r] with [v] volts across it
power(v,r) = v**2 / r

pa = var('pa') # power attenuation

# source/load impedances must match, and power
# over load and source must match attenuation
pi_eqns = [ 
    rpi_src(zload,a,b,c)==zsrc, 
    rpi_load(zsrc,a,b,c)==zload, 
    (power(pd_vout(v,b,c,zload), zload)) / power(v,zsrc) == pa  ]

# Design a -3dB attenuator for 50ohm source and 50ohm load
soln = solve( pi_eqns + [
    zload == 50, 
    zsrc == 50, 
    pa == 1/2, # power attenuation is half, ie. 3dB
    # we want the load-impedance seen by the source 
    # to match the source impedance and vice-versa
    ],
    a,b,c,pa,zload,zsrc,solution_dict=True)

soln

{a: -100*sqrt(2) + 150,  # 291ohm
 b: -25/2*sqrt(2),       # 17.7ohm
 c: -100*sqrt(2) + 150,  # 291ohm
 pa: 1/2,
 zload: 50,
 zsrc: 50},

So now we know what values our resistors need to be! Interestingly, A and C turn out to be the same value (this happens whenever source and load impedance match) but we did not have to bake in this as a special case .. it just falls out of the equations.

Since we’re using an equation solver, we can use the same setup to find out what the attenuation would be for given resistor values. Here we provide values for a/b/c resistors but omit zload, zsrc and the pa (power attenuation) and the equation solver figures them out.

soln = solve(pi_eqns + [
    a == -100*sqrt(2) + 150,
    b == -25/2*sqrt(2),
    c == -100*sqrt(2) + 150],
    a,b,c,pa,zload,zsrc,solution_dict=True)
soln

  pa: 1/2, # ie. power is halved
  zload: 50,
  zsrc: 50

And if we pick some random values for a/b/c, we can see the consequences for the impedance and the power attenuation:

soln = solve( pi_eqns + [
    a == 10,
    b == 20,
    c == 30],
    a,b,c,pa,zload,zsrc,solution_dict=True)

{a: 10,
  b: 20,
  c: 30,
  pa: -4*sqrt(5) + 9,  # ie. 0.0557280900008408
  zload: 6*sqrt(5),    # ie. 13.4164078649987
  zsrc: 10/3*sqrt(5)}, # ie. 7.45355992499930

So if you ever need an 18x power attenuation for your 13.4ohm source and 7.45ohm load, you’re in luck!