Firstly in the digital implementation they're not designed to mimic the amplitude response of the analog filter as Bessel designs are intended to give a desired amplitude response whilst maintaining the most linear phase and delay characteristics.
Secondly, while the polynomial used to define these filters is easy to compute there isn't a fixed algorithm for generating the poles as we have seen with the other filters. Let's take a look at the polynomial bit first. We're not going to go into a deep maths dive here, although the Bessel polynomials and related Bessel functions are a fascinating topic (Bessel functions are used to model the deformations of oscillating 2D planes, like a drum or tabla for example). The Bessel polynomial looks like this for the first couple of terms:
$$ \begin{align}
y_1(x) & = x + 1 \\
y_2(x) & = 3x^2 + 3x + 1 \\
y_3(x) & = 15x^4 + 15x^3 + 6x + 1 \\
y_4(x) & = 105x^5 + 105x^4 + 45x^3 + 10x + 1 \\
\end{align} $$
However, in electronic engineering, we're more usually working with what is known as the reverse Bessel polynomial which is a similar form with the coefficients reversed, which looks like this:
$$ \begin{align}
\theta_1(x) & = x + 1 \\
\theta_2(x) & = x^2 + 3x + 3 \\
\theta_3(x) & = x^4 + 6x^3 + 15x + 15 \\
\theta_4(x) & = x^5 + 10x^4 + 45x^3 + 105x + 105 \\
\end{align} $$
Where this can be generally described as:
$$ \theta_n(x) = x^ny_n\left(\frac {1}{n}\right) = \sum^n_{k=0} \frac {(2n-k)!}{(n-k)!k!} \frac {x^k}{2^{n-k}} $$
Very nice... now how about the formula for the poles and we can get right on make our filter. Unfortunately there isn't one! One of the usual methods is to generate the polynomial and then use a root-finding algorithm like Jenkins-Traub, Laguerre's method or good old Newton Raphson. However, the poles for a large number of orders are easily available in many places, so we're not going to generate lots of algorithmic code and instead just put the poles for orders up to 10 into a function.
So, here are the roots:
And here's the code:
public List<Complex> getPoles(int order) { Debug.Assert(order <= 10); List<Complex> poles = null; switch (order) { case 1: poles = new List<Complex>() { new Complex(-1, 0) }; break; case 2: poles = new List<Complex>() { new Complex(-1.5, 0.8660254038), new Complex(-1.5, -0.8660254038) }; break; case 3: poles = new List<Complex>() { new Complex(-1.8389073227, 1.7543809598), new Complex(-2.3221853546,0), new Complex(-1.8389073227, 1.7543809598) }; break; case 4: poles = new List<Complex>() { new Complex(-2.1037893972, 2.6574180419), new Complex(-2.8962106028 , 0.8672341289), new Complex(-2.8962106028 , -0.8672341289), new Complex(-2.1037893972, -2.6574180419) }; break; case 5: poles = new List<Complex>() { new Complex(-3.646738595, 0), new Complex(-2.324674303, 3.57102292), new Complex(-3.351956399, 1.742661416), new Complex(-3.351956399, -1.742661416), new Complex(-2.324674303, 3.57102292) }; break; case 6: poles = new List<Complex>() { new Complex(-2.515932248, 4.492672954), new Complex(-3.735708356, 2.626272311), new Complex(-4.248359396, 0.867509673), new Complex(-4.248359396, -0.867509673), new Complex(-3.735708356, -2.626272311), new Complex(-2.515932248, -4.492672954) }; break; case 7: poles = new List<Complex>() { new Complex(-4.971786859, 0), new Complex(-2.685676879, 5.420694131), new Complex(-4.070139164, 3.517174048), new Complex(-4.758290528, 1.739286061), new Complex(-4.758290528, -1.739286061), new Complex(-4.070139164, -3.517174048), new Complex(-2.685676879, -5.420694131) }; break; case 8: poles = new List<Complex>() { new Complex(-2.838983949, 6.353911299), new Complex(-4.368289217, 4.414442501), new Complex(-5.204840791, 2.616175153), new Complex(-5.587886043, 0.867614445), new Complex(-5.587886043, -0.867614445), new Complex(-5.204840791, -2.616175153), new Complex(-4.368289217, -4.414442501), new Complex(-2.838983949, -6.353911299) }; break; case 9: poles = new List<Complex>() { new Complex(-6.297019182, 0), new Complex(-2.979260798, 7.291463688), new Complex(-4.638439887, 5.317271675), new Complex(-5.60442182, 3.498156918), new Complex(-6.129367904, 1.737848384), new Complex(-6.129367904, -1.737848384), new Complex(-5.60442182, -3.498156918), new Complex(-4.638439887, -5.317271675), new Complex(-2.979260798, 7.291463688) }; break; case 10: poles = new List<Complex>() { new Complex(-3.108916234, 8.232699459), new Complex(-4.886219567, 6.224985483), new Complex(-5.967528329, 4.384947189), new Complex(-6.615290966, 2.611567921), new Complex(-6.922044905, 0.867665196), new Complex(-6.922044905, 0.867665196), new Complex(-6.615290966, 2.611567921), new Complex(-5.967528329, 4.384947189), new Complex(-4.886219567, 6.224985483), new Complex(-3.108916234, 8.232699459) }; break; } return poles; }
Which is pretty tedious, but works nicely. As a short digression, there are a number of vague posts looking at some algorithmic calculations for the poles but detail is too scant and although intrigued I was not able to track down any reliable description or get this to come to a conclusion myself. Diagrammatically the poles look as below when plotted on the s-plane:
If you're interested to just check that the poles all work out to the polynomial then just put the poles into the Polynomial.expand function we developed earlier, which provides a useful cross-check.
Great, so we now have the poles. Let's have a look at the amplitude frequency response for this analog prototype. We can do this by putting some numbers for the frequency (in radians) into the magnitude equation calculated from the analog transfer function (e.g. for a low-pass filter this is the Bessel polynomial in the denominator and is a poles only equation). Take a look back at one of the previous posts where we covered this before. The equations for Bessel functions can be found here. They look like this for a first few of the orders:
$$ \begin{align}
|H_2(\omega)| &= \frac {3}{\sqrt{\omega^4+3\omega^2+9}} \\
|H_3(\omega)| &= \frac {15}{\sqrt{\omega^6+6\omega^4+45\omega^2+225}} \\
|H_4(\omega)| &= \frac {105}{\sqrt{\omega^8+10\omega^6+135\omega^4+1575\omega^2+11025}}\\
\end{align}$$
As the coefficient size ramps up pretty quickly, in tabular form I've calculated them out for orders up to 10:
Giving the following nice looking graphs
Which look wonderful. However, we have a problem (and hence the reason for describing the magnitude response equations), the cut-off frequency is shifting for each of the orders unlike the Butterworth and Chebyshev responses where it in a similar location. So, we need some way of scaling the poles to normalise the cut-off frequency before shifting into the digital domain.
As the pole calculation has been tabularised and as this requires some additional computation, rather than doing this in code I put the magnitude response equations into a spreadsheet and then did a 'goal-seek' on the frequency to get a magnitude response of 0.5 (e.g. 3dB). This gives the following frequencies (in radians):
We can now use these to scale the poles before returning them, so the getPoles function is adjusted to scale (divide) by these values (I'll put the code at the end).
Next step our usual Low-Pass Filter (lpf) function can be reinstated. In this case the pre-warping needs to be removed as we're not going to use the bi-linear transformation and instead the matched Z-transform is going to be used instead. There is a lot of literature on this, so I'm not going to delve into that topic in this post.
The matched Z-transform is used again from Christian d'Heureuse's Java implementation of Tony Fischer's C code and has been added to the ComplexPlane class:
public ComplexPlane MatchedZTransform() { ComplexPlane zplane = new ComplexPlane(); zplane.poles = new List<Complex>(); zplane.zeros = new List<Complex>(); for (int i = 0; i < poles.Count; i++) { // returns exponential function of the complex number zplane.poles.Add(Complex.FromPolarCoordinates(Math.Exp(poles[i].Real),poles[i].Imaginary)); } for (int i = 0; i < zeros.Count; i++) { // returns exponential function of the complex number zplane.zeros.Add(Complex.FromPolarCoordinates(Math.Exp(zeros[i].Real),zeros[i].Imaginary)); } return zplane; }
We can now generate some filters, which when plotted look like this:
And finally here's the code:
using System; using System.Diagnostics; using System.Numerics; using System.Collections.Generic; // Hondrou implementation of a Bessel filter building on the STK framework // // This approach uses the techniques from Tony Fischer at York in mkfilter // http://www-users.cs.york.ac.uk/~fisher/mkfilter/ // and builds on the Java implementation by Christian d'Heureuse // www.source-code.biz, www.inventec.ch/chdh // namespace AgSynth { public class Bessel : Iir { public List<Complex> getPoles(int order) { Debug.Assert(order <= 10); List<Complex> poles = null; switch (order) { case 1: poles = new List<Complex>() { new Complex(-1, 0) }; poles = Polynomial.multiply(poles, 1 / 1.731939); break; case 2: poles = new List<Complex>() { new Complex(-1.5, 0.8660254038), new Complex(-1.5, -0.8660254038) }; poles = Polynomial.multiply(poles, 1 / 1.977607); break; case 3: poles = new List<Complex>() { new Complex(-1.8389073227, 1.7543809598), new Complex(-2.3221853546,0), new Complex(-1.8389073227, 1.7543809598) }; poles = Polynomial.multiply(poles, 1 / 2.422182); break; case 4: poles = new List<Complex>() { new Complex(-2.1037893972, 2.6574180419), new Complex(-2.8962106028 , 0.8672341289), new Complex(-2.8962106028 , -0.8672341289), new Complex(-2.1037893972, -2.6574180419) }; poles = Polynomial.multiply(poles, 1 / 2.886229); break; case 5: poles = new List<Complex>() { new Complex(-3.646738595, 0), new Complex(-2.324674303, 3.57102292), new Complex(-3.351956399, 1.742661416), new Complex(-3.351956399, -1.742661416), new Complex(-2.324674303, 3.57102292) }; poles = Polynomial.multiply(poles, 1 / 3.324236); break; case 6: poles = new List<Complex>() { new Complex(-2.515932248, 4.492672954), new Complex(-3.735708356, 2.626272311), new Complex(-4.248359396, 0.867509673), new Complex(-4.248359396, -0.867509673), new Complex(-3.735708356, -2.626272311), new Complex(-2.515932248, -4.492672954) }; poles = Polynomial.multiply(poles, 1 / 3.725753); break; case 7: poles = new List<Complex>() { new Complex(-4.971786859, 0), new Complex(-2.685676879, 5.420694131), new Complex(-4.070139164, 3.517174048), new Complex(-4.758290528, 1.739286061), new Complex(-4.758290528, -1.739286061), new Complex(-4.070139164, -3.517174048), new Complex(-2.685676879, -5.420694131) }; poles = Polynomial.multiply(poles, 1 / 4.090078); break; case 8: poles = new List<Complex>() { new Complex(-2.838983949, 6.353911299), new Complex(-4.368289217, 4.414442501), new Complex(-5.204840791, 2.616175153), new Complex(-5.587886043, 0.867614445), new Complex(-5.587886043, -0.867614445), new Complex(-5.204840791, -2.616175153), new Complex(-4.368289217, -4.414442501), new Complex(-2.838983949, -6.353911299) }; poles = Polynomial.multiply(poles, 1 / 4.423721); break; case 9: poles = new List<Complex>() { new Complex(-6.297019182, 0), new Complex(-2.979260798, 7.291463688), new Complex(-4.638439887, 5.317271675), new Complex(-5.60442182, 3.498156918), new Complex(-6.129367904, 1.737848384), new Complex(-6.129367904, -1.737848384), new Complex(-5.60442182, -3.498156918), new Complex(-4.638439887, -5.317271675), new Complex(-2.979260798, 7.291463688) }; poles = Polynomial.multiply(poles, 1 / 4.732941); break; case 10: poles = new List<Complex>() { new Complex(-3.108916234, 8.232699459), new Complex(-4.886219567, 6.224985483), new Complex(-5.967528329, 4.384947189), new Complex(-6.615290966, 2.611567921), new Complex(-6.922044905, 0.867665196), new Complex(-6.922044905, 0.867665196), new Complex(-6.615290966, 2.611567921), new Complex(-5.967528329, 4.384947189), new Complex(-4.886219567, 6.224985483), new Complex(-3.108916234, 8.232699459) }; poles = Polynomial.multiply(poles, 1 / 5.022755); break; } return poles; } // fcf is the cutoff frequency in Hz public void lpf(double fcf, int order) { double nfcf = fcf / PcmStreamSource.SampleRate; // normalised corner frequency Debug.Assert(nfcf > 0 && nfcf <= 0.5); // greater than 0 and less than Nyquist List<Complex> poles = getPoles(order); double w1 = 2 * Math.PI * nfcf; ComplexPlane splane = new ComplexPlane(); splane.poles = Polynomial.multiply(poles, w1); splane.zeros = new List<Complex>(); ComplexPlane zplane = splane.MatchedZTransform(); TransferFunction transferfunction = zplane.getTransferFunction(); Coefficients coeffs = transferfunction.getCoefficients(); double gain = transferfunction.getGainAt(Complex.One); // DC Gain // normalise the gain by dividing coeffs.b by gain coeffs.scaleGain(gain); setNumerator(coeffs.b); setDenominator(coeffs.a); } } }
Nice blog, but very anonymous. I would like to know more about you and your projects...
ReplyDeleteChristian d'Heureuse, www.source-code.biz
Hi Christian, thanks for the positive comment and interest. It was my New Year's resolution to start a blog with the things that interest me as I'd kept on putting it off due to lack of time (work and a young family) and decided that I just needed to get started and the momentum would build. To be honest I still have far too little time, and try to grab a bit of an evening, lunchtime, in hotels or airports to write about what interests me.
ReplyDeleteThe reason for keeping this anonymous, well, as the internet sticks with us and everyone is easily Googleable I want initially to keep my personal interests and technical blogs separate from my professional profile so that they are not confused. This means I can explore what I want personally and share that without any potential compromises.
My personal interests are basically anything media, music, art, video, graphics, typography and linguistics/language etc.. related to software and electronics. I have an undergraduate degree in electronic engineering and physics and have worked in software for over 22 years. My early expertise was in C++, for the past 15 years this has moved to Java (which I mostly manage rather than develop myself), C-sharp/Silverlight as the basic tool of preference, mainly because the machine I have with me the most is Windows and I can do application and web-app work in a single language. More recently, like many I'm enjoying Swift (everything outside of work is on Apple) and am interested to see what possibilities this has for Live-Coding.
I enjoy dabbling with Rasp-Pi and Ardunio if only to entertain the children and remember the fun of using a soldering iron! I'm a big fan of the Make movement and general wider hacker philosophy that is being driven by the internet, open-source, kickstarter funding, 3D printing and greater interest in coding. I have been ever thankful for the embedded software, microcontroller and discrete maths courses I did at Uni and work with a number of local Universities to bridge students to industry.
Professionally I manage Research and Development department for a well-known global electronics company specialising in software for the media and broadcasting industry.