Formula for Bessel low-pass filter coefficients
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:


Generate the coefficients of a reverse Bessel polynomial
Find the roots of the polynomial (= poles of the Bessel filter) using the Aberth-Ehrlich method
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)$:

 

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...):