Category: General

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!

I’ve been studying electromagnetism recently, and consequently been doing more maths-by-hand. Every time I do this, I think about using computer algebra systems to check my working and help me with calculus. But the only computer algebra system I’ve used was back at uni (Maple, I think) – and whilst Mathematica looks whizzy it’s also quite pricey. So I thought I’d give a few ‘free’ ones a go.

I have a particular ‘test’ problem in mind, because I’ve just worked through it by hand whilst reading Feynman Vol 2. It’s calculating the Laplacian of a radially symmetric potential – which involves a nice mixture of partial and total derivatives, chain rules and product rules. Turns out, I do actually remember how to do all this stuff, but having done it by hand it makes a nice concrete example for doing it on computer.

Another motivation is that I hate the imprecise notation that physicists and some mathematician use. What is the derivative of phi(r) with respect to x? My brain says that phi is a unary function with parameter r, so the only thing that can cause it to change is changes in it’s input r. In physics you have to ‘understand’ that this means phi is really a function over R^3 (space) that’s defined like phi(x,y,z) = g(sqrt(x^2+y^2+z^2)) with the helper g(r)=… telling you how things change with distance. And quite often, phi will also quietly become a function of time too. In SICM they fix this problem by using a scheme-based computer algebra system. Here I’m trying to do the same thing but using a more mainstream maths package.

I just ran sagemath via docker, rather than worry about ‘installing it’:

docker run -it sagemath/sagemath:latest

And so we start to tell it about our world …

%display latex

# We'll use g(r) as our "how things change with distance" function, and we'll leave it abstrac

g = function("g")

# Now we can define phi to be a concrete function over all space that delegates to g

phi(x,y,z) = g(sqrt(x^2+y^2+z^2))
phi.differentiate(x,2)
 
# result
# x^2*D[0, 0](g)(sqrt(x^2 + y^2 + z^2))/(x^2 + y^2 + z^2) - x^2*D[0](g)(sqrt(x^2 + y^2 + z^2))/(x^2 + y^2 + z^2)^(3/2) + D[0](g)(sqrt(x^2 + y^2 + z^2))/sqrt(x^2 + y^2 + z^2)

Assuming that D[0] means ‘derivative with respect to 0th arg’ then this is right but boy is it ugly!

Specifically, we’re not giving a short name ‘r’ to the expression sqrt(x^2+y^2+z^2). We can try to cajoule sagemath along:

r = var('r')
phi.differentiate(x,2).subs({sqrt(x^2+y^2+z^2): r})

Sadly that only improves the g(r) expressions, but not the usages in the denominators.

Perhaps we should tell it about r^2 rather than r?

phi.differentiate(x,2).subs({(x^2+y^2+z^2): r^2})

This does affect the denominator, but weirdly leaves lots of sqrt(x^2) terms. Odd.

We can blunder on ahead and complete the Laplacian:

(phi.differentiate(x,2) + phi.differentiate(y,2) + phi.differentiate(z,2)).simplify_full()

.. which yields the right answer, but again in a not-very-simplified form. Still, this is a pretty impressive feat. It takes a lots of pencil + paper to do these three second derivatives, chain rules and everything else. Doing it in one line shows that sagemath is a pretty useful maths assistant!

So far, we’ve just modelled phi as a three-arg function – but sagemath actually knows what scalar and vector fields are, and can already do operations on them!

from sage.manifolds.operators import *
E.<x,y,z> = EuclideanSpace()
g = function("g")
phi = E.scalar_field(g(sqrt(x^2+y^2+z^2)), name='phi')
phi.laplacian().display()

… which gives the same answer as above, but now we’ve got “batteries included” and can calculate div/grad/curl and all that.

So this is now good enough to solve electrostatic problems, where the laplacian of the scalar potential is proportional to charge density. For dynamic cases however, we’ll need a function of x,y,z and also t. Not sure how that’ll work with the manifold support. Do I need a function from time to scalar fields? A field of time-to-scalar functions? A 4d spacetime structure?

This page is working with E/B fields using 4d manifolds. I’ve done courses on special relativity, so I can handle spacetime but the language of differential geometry (manifolds, 2-forms) is beyond me currently (although Sussman of SICP fame has a similarly computational approach to it). Perhaps it’s better to leave the sage.manifolds alone for now and return to explicit x/y/z/t approach (hey, it was good enough for Maxwell!).

When I was little, I played around with batteries and wires and various electrical things that I’d taken to bits – enough to both burn my fingers (from shorting the battery) and give myself shocks (using a transformer to transiently step up my 9v battery to much higher voltages). Consequently, I was familiar with the word “voltage” and the vague hand-wavy description of it “being like water pressure” or “being like height”.

But then when I learned physics, I found that the primary focus was on the electric (really electrostatic) field – being a 3d vector field. Voltage aka potential difference made an appearance here, thanks to the nice property that for electrostatic fields you can think in terms of a scalar field (“the potential”), with the electric field being its gradient.

But it seemed like a big gulf between the “circuit theory viewpoint” and the “electrostatic viewpoint”. For example, if I take a lump of neutral material and grab some of its electrons and move them over a bit, I end up with a positively charged lump and a negative lump, and the surrounding E-field is a particular curvy shape (aka a “dipole field”). If I then think of this as a battery, and try to connect some wires to it, then my “circuit theory” viewpoint says that the electric field inside the wire must point along the wire the whole way – even if I tangle the wire into crazy knots. But how do we get from the nice dipole field to the twisty turny field required to follow the wire? It seems like two different worlds. Fortunately, there is an explanation – just not one that my Uni physics book mentioned. The best description of this comes from Matter and Interaction textbook (eg. here). When the wires are connected, there is an initial transient in which charges are displaced by the E-field, in such a way that a charge imbalance is created on the surfaces of the wires which causes the local E-field to reorient itself along the wire.

This is good news, because it means that you can indeed use all your electrostatic intuition of voltage and electric fields when dealing with circuits with batteries and resistors. The electric vector fields are all of the simple variety where you can succinctly summarise them using scalar fields. One key consequence of these ‘simple’ fields is that the work done by (or against) the electric field between two points is the same regardless of the route you take. This means we can talk about “the” voltage between point A and point B even when there’s multiple paths between them, because even though the E-field might be different along different paths, the net work done by (or against) the E-field along the entire path is always the same. Going one step further, if we agree to use a particular point (often labelled “ground”) as our reference point, we can even talk about “the voltage at point A” – but we understand that to mean “the potential difference between A and our chosen reference point”.

Voltmeters essentially rely on this “path independence” to do their job. In our nice electrostatic world, it doesn’t matter that we connect our voltmeter to points A and B using long winding leads, because even though the E-field might vary in space, we KNOW that the net work done by the E-field on a charge going through our voltmeter is the same as would be done by a charge going through the resistor (or whatever component we’re using the voltmeter to measure). Note: most measuring devices are named alluringly to suggest they directly measure the property you care about, when in practise they measure something else and use it to infer the value of interest. For example, an analog voltmeter actually measures the current through a known resistance, and we trust/hope that this is the same as the “volts” across whatever we’ve connected the leads to).

However, this lovely simple land of electrostatic fields with their path-independence isn’t the whole story. You can also create an electric field from a time-varying magnetic field (Faraday’s Law) and they’re of a very different nature – they form loops, quite unlike electrostatic fields which never form loops. With this kind of field, the work done going from point A to point B now DOES depend on exactly which path you choose. These fields cannot be represented as the gradient of a scalar potential. Consequently, we lose the lovely simple world in which we could think of point A as just having a single number (the potential) and point B having some other number/potential, and the work done by moving between then is just potentialB-potentialA – all of this is gone!

We do still have some useful properties – we can relate the work done around any closed path to the changing magnetic field (its flux) over a surface ‘capping’ that loop. But the notion of voltage as some path-independent value between point A and point B is lost. If we connect a “voltmeter” to points A and B, we’ll measure something pertaining to the path taken by the voltmeter’s leads and the voltmeter itself. But out ability to say “and that will be the same as any other path (such as through the component under test)” is gone. If we physically move the voltmeter and its leads in space, we will get a different reading (depending on the changing flux within the ‘loop’ consisting of the voltmeter and its leads and the circuit under test) – even though it is still connected to the same pointA and pointB in the circuit.

This paper is a great summary of what’s going on, and I was relieved to finally find it after many years to failing to reconcile my “physics” and “circuit” viewpoints! This other paper is also very good.

Incidentally, this makes it clear what Kirchhoff’s voltage law is coming from. It’s just a restatement of the properties of an electrostatic field – if you go from pointA to pointA in an electrostatic field, there’s no work done on a charge. But Kirchhoff came up with his rule long before anyone started thinking in terms of electric fields. It was also well before Faraday started futzing with induction and creating electric fields for which that property doesn’t hold. Essentially, Kirchhoff’s voltage law is a special case of Faraday’s law when there is no changing magnetic flux. If there is a changing magnetic field, the work done around a loop (the emf) is non-zero and Kirchhoff’s Law no longer holds.