1
$\begingroup$

When I am filtering a signal in python, it has a built in function to generate bessel filter coefficients given a cutoff ratio and a filter order (number of poles). I am trying to translate this to C code, but I cannot seem to find the formula that is used to calculate the coefficients.

Can someone point me to the formulas to calculate the IIR filter coefficients given a cutoff ratio and order?

I am trying to read the scipy source on github but I am having a very hard time of it...

$\endgroup$
3
$\begingroup$

Looking at the source code for scipy.signal.bessel, you're out of luck for finding a formula: they don't use one.

They just have a big if \ elif sequence for various values of filter order:

if N == 0:
    p = []
elif N == 1:
    p = [-1]
elif N == 2:
    p = [-.8660254037844386467637229 + .4999999999999999999999996j,
         -.8660254037844386467637229 - .4999999999999999999999996j]
elif N == 3:
    p = [-.9416000265332067855971980,
         -.7456403858480766441810907 - .7113666249728352680992154j,
         -.7456403858480766441810907 + .7113666249728352680992154j]
elif N == 4:

The Wikipedia page has some formulae for Bessel filters in the continuous domain: $$ \theta(s) = \sum_{k=0}^n a_k s^k $$ where $$ a_k = \frac{(2n-k)!}{2^{n-l}k!(n-k)!} \mbox{ for } n = 0,1,\ldots,n. $$ To get discrete-time near-equivalents, you'll need to do a continuous-to-discrete transformation. Note that this will not necessarily preserve the nice phase properties of the continuous-time Bessel filter.

$\endgroup$
  • $\begingroup$ I don't mind using a lookup table for my own code. I need to dig through the rest of their code to figure out to transform for a given cutoff. Those arrays are for the general case and have not yet been shifted for a given cutoff, I think. $\endgroup$ – KBriggs Apr 25 '16 at 19:28
  • 1
    $\begingroup$ Yes. I haven't looked at that code in detail, but usually a unit cutoff is assumed for the general case. $\endgroup$ – Peter K. Apr 25 '16 at 19:38
  • 1
    $\begingroup$ @KBriggs The frequency shifting is done in the function called _zpklp2lp in that source code. $\endgroup$ – Peter K. Apr 26 '16 at 16:02
  • $\begingroup$ Thanks, I managed to pull out the relevant parts of the scipy code, so I now have my own functional python implementation. Translating it to C is going to be a bit of a pain because of the complex numbers, but it all looks doable. $\endgroup$ – KBriggs Apr 26 '16 at 16:05
  • $\begingroup$ Actually you're looking at an outdated version. I rewrote it in Sept 2015 github.com/scipy/scipy/commit/… You have to look at the "master" branch to see the latest code. $\endgroup$ – endolith May 5 '17 at 23:30
3
$\begingroup$

Hi, I wrote this. It's probably more complicated/inefficient than it needs to be. :D

Practical Bessel filter design involves root-finding of a polynomial to generate second-order sections; I don't believe there is any simple formula. The way SciPy stores filter coefficients is pole and zero locations, so the code finds these locations numerically.

The steps I originally used are:

  1. Generate the coefficients of a reverse Bessel polynomial
  2. Find the roots of the polynomial (= poles of the Bessel filter) using the Aberth-Ehrlich method
  3. Scale the poles inward or outward to normalize either the -3 dB frequency, phase, or group delay.

Step 2 generates approximate roots from Campos-Calderón 2011 as the starting points for the root-finding, but I don't know if that's really necessary. Most root-finding algorithms just use asymmetrical starting points on a spiral, etc. and still find the answer fine. I assumed it would be faster this way.

Aberth-Ehrlich is good for quickly finding multiple single roots simultaneously, by modelling them as point charges that repel each other, so multiple test points don't fall into the same holes.

Eventually I ended up eliminating Step 1 altogether. For the Newton's method part of Aberth-Ehrlich, it actually evaluates $K_{n+\frac 1 2}\left(\frac 1 x \right)$ instead, which is a totally different (compiled) function that has the same roots, and has a relatively simple derivative. This is possibly crazy and less efficient than it could be, but it works. (kve is actually zbesk from AMOS Fortran library.)

$\theta_{5}(x)$ vs $K_{5.5}\left(\frac 1 x \right)$:

Complex plot of 5th order reverse Bessel polynomial, with amplitude as lightness and phase as color Complex plot of 5.5th order modified Bessel function of 1/x, with amplitude as lightness and phase as color

Then in step 3, I don't even use the Bessel polynomial to normalize the function, since I eliminated the need for all but the last coefficient and it's faster to just generate that one by itself.

There's also this algorithm I found later, which could probably be made more efficient: Orchard 1965 - The Roots of the Maximally Flat-Delay Polynomials (Bessel filters) though it has some numerical error at high N and doesn't give an algorithm for finding the new approximate roots from the old (just "eyeball it on graph paper").

But now you can generate 500th order Bessel filters with 3 different normalizations, instead of just reading the coefficients from a table (not that you would ever need to...):

Poles of a 500th-order Bessel filter

$\endgroup$
  • 1
    $\begingroup$ the wikipedia section. just a disclaimer that it doesn't have any more claim to truth than my comments do :) $\endgroup$ – endolith May 6 '17 at 2:14
  • 1
    $\begingroup$ oh, i don't think either of us would write anything untruthful in either Wikipedia or SE. :-) $$ $$ and i've got the scars (from Wikipedia and from the Christianity.SE) to show for it. funny, i was fighting atheist POV-pushers at Wikipedia and Christians at the SE page. both groups seem to highly overestimate their own objectivity. $\endgroup$ – robert bristow-johnson May 6 '17 at 2:15
  • 1
    $\begingroup$ Thanks for the thorough explanation! If you're curious, my C port is here: github.com/shadowk29/CUSUM/blob/master/bessel.c $\endgroup$ – KBriggs May 12 '17 at 14:12
  • 1
    $\begingroup$ @KBriggs You saw the above discussion about Thiran filters, right? You shouldn't use bilinear-transformed Bessel unless the cutoff is significantly below fs/4 dsp.stackexchange.com/a/8420/29 I do have a working Thiran filter design function, but I need to clean it up a lot and convert it to using poles/zeros natively before putting it in SciPy. $\endgroup$ – endolith May 12 '17 at 14:21
  • 1
    $\begingroup$ I suspect it won't matter too much, but I'll look into it to be sure. Thanks for the heads up, that's new information for me. I'll look into Thiran when it hits scipy :) $\endgroup$ – KBriggs May 12 '17 at 14:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.