Trouble with filter using Laplace transfer function in LTSpice

Question

I would like to introduce an ideal 5th-order Bessel filter to my simulation. My initial thought is to use a voltage-dependent voltage source and Laplace transform. I copied the Bessel filter coefficients from Wikipedia. Bessel LPF filter coefficients.

s+1

s^2+3s+3

s^3+6s^2+15s+15

s^4+10s^3+45s^2+105s+105

s^5+15s^4+105s^3+420s^2+945s+945

I used coefficients to write down the transfer function: How it is written in LTspice:

Laplace = 945 / (pow((s/w0),5)+15 * pow((s/w0),4)+105 * pow((s/w0),3)+420 * pow((s/w0),2)+945 * (s/w0)+945)

[Spaces between integers and “pow” are there so the “*” sign appears in the LTspice formula spaces are not present.]

945 in the numerator corresponds to a DC gain of 0 dB.

ω = 2πf(-3 dB)

f(-3dB) = 3.5 kHz

ω = 2 * π * 3500

However, after simulating the circuit I have noticed that it does not behave as it should. My passband is going far beyond 3.5 kHz. My -3dB point is at 8.5 kHz. However, the behavior (even group delay in passband) is of the Bessel filter. Therefore, I scaled down the ω by the ratio of 3.5 kHz/8.5 kHz: ω = 2 * π * 1440. This resulted in a -3 dB point at 3.5 kHz and even group delay in the passband. This approach however is a bodged one to fit the desired result without fixing the fundamental problem somewhere in the transfer function.

I welcome all the suggestions as to where the mistake might be.

Answer 1

Save some grief: set \$\omega=1\$

Just save yourself some algebra grief and set \$\omega=1\$ when focused on a particular filter shape (like Bessel.) The filter shape can be moved around freely (at least mathematically, if not in practical terms) on the frequency axis.

So just give yourself a break. Don't fret \$\omega\$. Set it to 1.

That allows you to focus on the remaining important details without getting bogged down with extra symbols and it will be easier to think more clearly. (Super-genius types exempted.)

(NOTE: In the following, I'll use freely available SymPy and SageMath and Python.)

Butterworth at \$\omega=1\$ is \$-3.0103\:\text{dB}\$ regardless of order

It's possible that you became confused by the Butterworth polynomial family. It has a unique property that makes it special in terms of its \$-3.0103\:\text{dB}\$ cutoff.

Let's look at 2nd, 3rd, 4th, 5th, and 6th order Butterworth polynomials in order to illustrative this unique property. I'm going to ask SymPy to tell me what value of \$\omega\$ provides the \$-3.0103\:\text{dB}\$ cutoff:

# Butterworth
#
solve(Eq(abs(expand(prod(Butterworth(2))).subs({s:I*w})),sqrt(2)),w)[0].n()
1.00000000000000
solve(Eq(abs(expand(prod(Butterworth(3))).subs({s:I*w})),sqrt(2)),w)[0].n()
1.00000000000000
solve(Eq(abs(expand(prod(Butterworth(4))).subs({s:I*w})),sqrt(2)),w)[0].n()
1.00000000000000
solve(Eq(abs(expand(prod(Butterworth(5))).subs({s:I*w})),sqrt(2)),w)[0].n()
1.00000000000000
solve(Eq(abs(expand(prod(Butterworth(6))).subs({s:I*w})),sqrt(2)),w)[0].n()
1.00000000000000

As you can see, it doesn't matter what order I choose. A Butterworth polynomial always has a \$-3.0103\:\text{dB}\$ corner right at \$\omega=1\$.

It's my guess that you've been given to believe that all filter polynomials have this same feature. But the reality is that no other filter polynomial does.

However, a note: That doesn't mean that other filter polynomials cannot be frequency-scaled so that they do meet this behavior. But then the polynomial becomes order-dependent and will no longer look the same as its canonical polynomial.

Canonical Bessel polynomials

Before I apply SymPy on the canonical Bessel polynomials you provided, let's first take a look at curves (partly from Desmos and partly using Paint) that plot out where the canonical Bessel functions intersect with a \$-3.0103\:\text{dB}\$ horizontal line:

It's plain to see above that different Bessel orders intersect at different values of \$\omega\$.

Now, we are ready to let SymPy quickly do the calculations for these canonical Bessel polynomials you provided:

# Bessel
#
solve(Eq(abs(((s**2+3*s+3)/3).subs({s:I*w})),sqrt(2)),w)[0].n()
1.36165412871613
solve(Eq(abs(((s**3+6*s**2+15*s+15)/15).subs({s:I*w})),sqrt(2)),w)[0].n()
1.75567236868121
solve(Eq(abs(((s**4+10*s**3+45*s**2+105*s+105)/105).subs({s:I*w})),sqrt(2)),w)[0].n()
2.11391767490422
solve(Eq(abs(((s**5+15*s**4+105*s**3+420*s**2+945*s+945)/945).subs({s:I*w})),sqrt(2)),w)[0].n()
2.42741070215263

(I skipped order 1 because it is exactly the same shape as the Butterworth.)

Note that none of these came up with \$\omega=1\$. And the last one says that \$\omega=2.42741070215263\$ when the transfer function is \$-3.0103\:\text{dB}\$.

That should be at \$f=0.386334412161760\$. So let's take a quick peek at what LTspice says:

Yeah. I'd say that's a match.

Your expectations about these canonical Bessel polynomials are the problem. That's all. (Likely because the canonical Butterworth polynomials led you into conflated thinking about all filter shapes.)

Frequency-scaling Bessel polynomials

Practical filter designers focused on solving a problem at hand and looking for pre-computed filter tables to save time will want those tables set up so that they provide \$-3.0103\:\text{dB}\$ at \$\omega=1\$ and not where the canonical form would instead put it.

So, in response, textbook authors typically provide table values for each order with filter-order-based adjustments already applied so that the expected \$-3.0103\:\text{dB}\$ cutoff does occur at \$\omega=1\$.

But such adjustments also leave you with a different polynomial, as well. It will no longer match up with the canonical form.

To make this point, let's look at one of those practitioner-tables from the 2nd edition of "Intuitive Analog Circuit Design" by Dr. Marc T. Thompson:

(I had to use Paint to fix errors in the table layout before posting it.)

Let's pick the 3rd order entry and take their parameters to see what we get as the polynomial:

expand((s-(-1.0509-I*1.0025))*(s-(-1.0509+I*1.0025))*(s-(-1.327)))
s**3 + 3.4288*s**2 + 4.89848566*s + 2.79916989862

Hmm. Now, that's definitely not the polynomial you listed for 3rd order!

(It's also slightly wrong. More on that, later.)

To see what the table's author did in preparing that table, let's take those \$\omega\$-factors I computed with SymPy above and use them to make filter-order-based adjustments that I think the author applied to create that table.

Here's those frequency-adjusted Bessel functions where they have been scaled (stretched on the frequency axis) so that they all wind up crossing at \$\omega=1\$:

These are still Bessel functions appropriate for their order. But they are no longer in their canonical form. Each order has been stretched on their frequency axis, and differently for each based upon its order, so that they wind up crossing where a practical designer would likely want.

But as I earlier showed with the 3rd order case from the table, the resulting polynomial isn't canonical anymore. It's been moved around on the frequency axis. Which explains the differences in constants.

If all that is true, then how do you know it's really a Bessel anymore?

Well, there are two constants present in the 2nd order version. It's how they relate together that tell you what kind of filter shape you have. And for higher order, that's still true. It's just that you have more constants and how they all relate to each other still matters.

It's time to investigate that.

Filter shape

Filter shape is easier to follow with 2nd order because a 2nd order filter has only two basic ideas: \$\omega_{_0}\$ and either \$Q\$ or \$\zeta\$ (each determines the other one.) \$\omega_{_0}\$ just moves the shape left or right on the frequency axis. It doesn't change the shape. So the only parameter determining shape with 2nd order is \$Q\$ or \$\zeta\$. That makes it easy.

Let's test this idea with the 2nd order Bessel. I'll use the canonical values and then also the re-positioned values that I used in the above plot.

Let's give it a whirl:

cbessel = 3/(s**2+3*s+3)                    # canonical 2nd order Bessel

k0 = 2/(-1 + sqrt(5))
k1 = sqrt(6)/sqrt(-1 + sqrt(5))
sbessel = k0/(s**2+k1*s+k0)                 # scaled 2nd order     
sbessel.n()
1.61803398874989/(s**2 + 2.20320266118432*s + 1.61803398874989)

nsimplify(tf2(cbessel)[zeta])
sqrt(3)/2                                   # canonical 2nd order shape

nsimplify(tf2(sbessel)[zeta])
sqrt(3)/2                                   # scaled 2nd order shape

You can see that the shape didn't change. It's the exact same shape. So it is still a Bessel. But the constants are certainly different!! That's because the scaled 2nd order Bessel was frequency-stretched in just the right way to get the cross-over where, as a practical designer, you'd likely want it placed.

I hope that's enough to help you at least accept the idea here; that polynomials may need to be scaled/stretched on the frequency axis to help designers achieve particular goals. And more, that textbook authors will often reflect such desires when laying out design tables.

But at least one thing should now be evident: mathematical families of polynomials may have canonical forms more commonly used by mathematicians and frequency-stretched versions more commonly used by electronics designers and they may look differently -- especially for any family other than Butterworth (which is special in this way.)

How to develop a frequency-stretched Bessel polynomial

So, this beggars another question:

How did I come up with \$\frac2{\sqrt{5}\,-\,1}\$ and \$\sqrt{\frac6{\sqrt{5}\,-\,1}}\$ for a frequency-stretched 2nd order Bessel?

Well, I am not reading this from any textbook (nor have I ever.) In fact, I don't even have an idea what textbook to look in. It's just how I'd handle it. Maybe it works for you to see it in action:

$$\begin{align*} H\left(s\right)&=\frac3{s^2+3\,s+3} \\\\&=\frac{\omega_{_0}^{\,2}}{s^2+2\zeta\,\omega_{_0}\,s+\omega_{_0}^{\,2}}, \text{ where }\omega_{_0}=\sqrt{3} \text{ and } \zeta=\frac{\sqrt{3}}2 \end{align*}$$

Note that the shape matches earlier computation using SymPy.

For the frequency-scaling, we must first find out how much to scale the frequency. I already did that far above from here. But let's do it again, but this time exactly:

nsimplify(solve(Eq(abs(((s**2+3*s+3)/3).subs({s:I*w})),sqrt(2)),w)[0])
sqrt(-3/2 + 3*sqrt(5)/2)

Let's scale all frequencies by that factor, calling it \$K_f=\sqrt{\frac32\left[\vphantom{\frac32}\sqrt{5}-1\right]}\$ so that \$\hat{s}=K_f\cdot s\$.

$$\begin{align*} H\left(\hat{s}\right)&=\frac{\omega_{_0}^{\,2}}{\hat{s}^2+2\zeta\,\omega_{_0}\,\hat{s}+\omega_{_0}^{\,2}} \\\\ &=\frac{\omega_{_0}^{\,2}}{\left(K_f\,s\right)^2+2\zeta\,\omega_{_0}\,\left(K_f\,s\right)+\omega_{_0}^{\,2}} \\\\ &=\frac{\omega_{_0}^{\,2}}{K_f^{\,2}\,s^2+2\zeta\,\omega_{_0}\,K_f\,s+\omega_{_0}^{\,2}} \\\\ &=\frac{ \frac1{K_f^{\,2}} }{ \frac1{K_f^{\,2}}}\cdot\frac{\omega_{_0}^{\,2}}{K_f^{\,2}\,s^2+2\zeta\,\omega_{_0}\,K_f\,s+\omega_{_0}^{\,2}} \\\\ &=\frac{\left(\frac{\omega_{_0}}{K_f}\right)^2}{s^2+2\zeta\left(\frac{\omega_{_0}}{K_f}\right)s+\left(\frac{\omega_{_0}}{K_f}\right)^2} \\\\ &=\frac{\frac2{\sqrt{5}\,-\,1}}{s^2+\sqrt{\frac6{\sqrt{5}\,-\,1}}\,s+\frac2{\sqrt{5}\,-\,1}} \end{align*}$$

That last denominator factors as \$s_\pm=\frac14\sqrt{2\vphantom{\sqrt{5}+1}}\sqrt{\sqrt{5}+1}\left[\vphantom{\sqrt{\sqrt{5}+1}}-\sqrt{3}\pm j\right]\$. There's only a factor of \$\sqrt{3}\$ difference between the two constants: \$s_\pm=-1.10160133059216\pm j\, 0.636009824757034\$.

This is a good moment to re-visit the table I pulled from a textbook where the author supplied \$-1.1030\pm j\, 0.6368\$. These constants also differ from each other by about a factor of \$\sqrt{3}\$. And that's important.

But while they are similar in possessing the same \$\sqrt{3}\$ ratio, the constant magnitudes are not the same. And that does have a (slight) impact on frequency scaling.

Let's perform a simple test to see which of the two is right.

(tf2(((c1**2+c2**2)/((s-c1+I*c2)*(s-c1-I*c2))).subs({c1:-1.103,c2:0.6368}))[zeta]).n()
0.866031301401309                  # zeta from textbook table

(tf2(((c1**2+c2**2)/((s-c1+I*c2)*(s-c1-I*c2))).subs({c1:-1.10160133059216,c2:0.636009824757034}))[zeta]).n()
0.866025403784438                  # zeta from my process

(sqrt(3)/2).n()
0.866025403784439                  # true zeta for 2nd order Bessel

Both are close and especially also share the same relative ratio.

The key here is how one might frequency-scale a filter so that its \$-3.0103\:\text{dB}\$ cutoff occurs when \$\omega=1\$.

As a final note in this section, let's see if there is a simpler way to frequency-scale. I'll use the 3rd order Bessel for that. Look back up to find the scaling factor:

# Compute new factors (s^3 factor is always 1):
(1/1.75567236868121)**1 * 6           # s^2 factor is scaled once
3.41749412192832
(1/1.75567236868121)**2 * 15          # s factor is scaled twice
4.86636086392275
(1/1.75567236868121)**3 * 15          $ units factor is scaled thrice
2.77179327460633

# Plug these in and work out the roots:
roots(s**3 + 3.41749412192832*s**2 + 4.86636086392275*s + 2.77179327460633,s)
{-1.32267579991045: 1,
 -1.04740916100894 - 0.999264436280635*I: 1,
 -1.04740916100894 + 0.999264436280635*I: 1}

That's it. We get about \$-1.3227\$ and \$-1.04741\pm j\,0.9993\$.

By comparison, the textbooks's table values are \$-1.327\$ and \$-1.0509\pm j\,1.0025\$. These also don't match in magnitude.

It's time to discuss errors.

Note about errors

The above errors in the textbook table are minor. And while getting errors in the 3rd decimal place seems a bit odd to me there's probably no real harm here.

A point I'd like to make is that textbooks can make errors that have greater impact. Sometimes, the errors are just wrong. So it helps to be able to find your own way, if and when needed.

(I have tracked these constants at least going back to a 1997/1998 edition of "Filter Design" by Steve Winder [Newnes Press 1998, also EDN 1997, I think.] But it may go well before then, too.)

Although these errors are minor, let's use LTspice to make a quick comparison. (It's useful to know how to do this in LTspice, so that errors may be quickly identified.)

I'll test the textbook version against the one I developed above:

As you can see, it's not a lot of difference.

Still? Never rest easy. Trust, but verify.

Answer 2

I suppose that the mentioned frequency (f=3.5 kHz) is the desired 3dB-frequency (your requirement), right?

In this case, it is not true that the corresponding w value (2Pi*f) would appear in the denominator of the 5th-order transfer function. Instead, it would be necessary to find the magnitude of the function and solve it for the 3-dB point - a very involved calculation.

Therefore: The classical strategy is to use a chain of 2nd-order stages and (if necessary) a first-order section. For such a configuration you can find all necessary information (pole data) in normalizes filter tables which can be applied for a certain passband edge (3dB frequency).

Answer 3

The Bessel coefficients from Wikipedia yield a filter with unit group delay, not unit cutoff frequency. It is thus not surprising that you don't get the cutoff frequency you expect.

What is your application? Bessel filters are commonly employed when smoothing signals in the time domain. If frequency domain filter requirements are driving the design, a different form of filter response may be more appropriate.