f and Q values for implementing Bessel filter as second-order sections (and Bessel polynomials in Python)

All values.txt

Q for N = 
 1: --------------
 2: 0.57735026919 
 3: -------------- 0.691046625825
 4: 0.805538281842 0.521934581669
 5: -------------- 0.916477373948 0.563535620851
 6: 1.02331395383  0.611194546878 0.510317824749
 7: -------------- 1.12625754198  0.660821389297 0.5323556979  
 8: 1.22566942541  0.710852074442 0.559609164796 0.505991069397
 9: -------------- 1.32191158474  0.76061100441  0.589406099688 0.519708624045
10: 1.41530886916  0.809790964842 0.620470155556 0.537552151325 0.503912727276
11: -------------- 1.50614319627  0.858254347397 0.652129790265 0.557757625275 0.513291150482
12: 1.59465693507  0.905947107025 0.684008068137 0.579367238641 0.525936202016 0.502755558204
13: -------------- 1.68105842736  0.952858075613 0.715884117238 0.60182181594  0.540638359678 0.509578259933
14: 1.76552743493  0.998998442993 0.747625068271 0.624777082395 0.556680772868 0.519027293158 0.502045428643
15: -------------- 1.84821988785  1.04439091113  0.779150095987 0.648012471324 0.573614183126 0.530242036742 0.507234085654
16: 1.9292718407   1.08906376917  0.810410302962 0.671382379377 0.591144659703 0.542678365981 0.514570953471 0.501578400482
17: -------------- 2.0088027125   1.13304758938  0.841376937749 0.694788531655 0.609073284112 0.555976702005 0.523423635236 0.505658277957
18: 2.08691792612  1.17637337045  0.872034231424 0.718163551101 0.627261751983 0.569890924765 0.533371782078 0.511523796759 0.50125489338 
19: -------------- 2.16371105964  1.21907150269  0.902374908665 0.741460774758 0.645611852961 0.584247604949 0.544125898196 0.518697311353 0.504547600962
20: 2.23926560629  1.26117120993  0.932397288146 0.764647810579 0.664052481472 0.598921924986 0.555480327396 0.526848630061 0.509345928377 0.501021580965
21: -------------- 2.31365642136  1.30270026567  0.96210341835  0.787702113969 0.682531805651 0.613821586135 0.567286339654 0.535741766356 0.515281087097 0.503735024056
22: 2.38695091667  1.34368488961  0.991497755204 0.81060830488  0.701011199665 0.628878390935 0.57943181849  0.545207253735 0.52208637596  0.507736060535 0.500847111042
23: -------------- 2.45921005855  1.38414971551  1.02058630948  0.833356335852 0.71946106813  0.64404152916  0.591833716479 0.555111517796 0.529578662133 0.512723802741 0.503124630056
24: 2.53048919562  1.42411783481  1.04937620183  0.85593899901  0.737862159044 0.659265671705 0.604435823473 0.565352679646 0.537608804383 0.51849505465  0.506508536474 0.500715908905
25: -------------- 2.60083876344  1.46361085888  1.07787504693  0.878352946895 0.756194508791 0.674533119177 0.617161247256 0.575889371424 0.54604850857  0.524878372745 0.510789585775 0.502642143876

f multiplier to get same asymptotes as Butterworth (LPF and HPF are phase-matched), for N = 
 1: 1.0*          
 2: 1.0            
 3: 0.941600026533* 1.03054454544  
 4: 1.05881751607   0.944449808226 
 5: 0.926442077388* 1.08249898167   0.959761595159 
 6: 1.10221694805   0.977488555538  0.928156550439 
 7: 0.919487155649* 1.11880560415   0.994847495138  0.936949152329 
 8: 1.13294518316   1.01102810214   0.948341760923  0.920583104484 
 9: 0.915495779751* 1.14514968183   1.02585472504   0.960498668791  0.926247266902 
10: 1.1558037036    1.03936894925   0.972611094341  0.934100034374  0.916249124617 
11: 0.912906724455* 1.16519741334   1.05168282959   0.984316740934  0.942949951193  0.920193132602 
12: 1.17355271619   1.06292406317   0.99546178686   0.952166527388  0.92591429605   0.913454385093 
13: 0.91109146642*  1.18104182776   1.07321545484   1.00599396165   0.961405619782  0.932611355794  0.916355696498 
14: 1.18780032771   1.08266791426   1.0159112847    0.970477661528  0.939810654318  0.920703208467  0.911506866227 
15: 0.909748233727* 1.19393639282   1.09137890831   1.02523633131   0.97927996807   0.947224181054  0.925936706555  0.91372962678  
16: 1.19953740587   1.09943305993   1.03400291299   0.987760087301  0.954673832805  0.93169889496   0.917142770586  0.910073839264 
17: 0.908714103725* 1.20467475213   1.10690349924   1.04224907075   0.995895132988  0.962048803251  0.937756408953  0.921340810203  0.911830715155 
18: 1.20940734995   1.11385337718   1.05001339566   1.00367972938   0.969280631132  0.943954185964  0.926050333934  0.91458037566   0.908976081436 
19: 0.907893623592* 1.21378428545   1.12033730486   1.05733313013   1.01111885678   0.97632792811   0.950188250195  0.931083047917  0.918020722296  0.910399240009 
20: 1.21784680466   1.1264026302    1.0642432684    1.0182234301    0.983167015494  0.956388509221  0.936307259119  0.921939345182  0.912660567121  0.90810906098  
21: 0.907227918651* 1.22162983803   1.13209053887   1.0707761723    1.02500765809   0.989785773052  0.962507744041  0.941630539095  0.926182679425  0.915530785895  0.909284134383 
22: 1.22516318078   1.13743698469   1.07696149672   1.03148734658   0.996179868375  0.968513835739  0.946988833876  0.930637530698  0.918842345489  0.911176680853  0.90740679425  
23: 0.906504834665* 1.22847241656   1.14247346478   1.08282633643   1.03767860165   1.00235036472   0.974386503777  0.952333544924  0.93521970212   0.922499715897  0.913522360451  0.908537315966 
24: 1.23157964663   1.14722769588   1.0883951254    1.04359863738   1.00829752657   0.980122054682  0.957614328211  0.939902288094  0.926326979303  0.916382125274  0.91007413458   0.906792423356 
25: 0.90735557785*  1.23450407249   1.15172412685   1.09369031049   1.04926215576   1.01403218959   0.985701250074  0.962823028773  0.944627742698  0.930311087183  0.919369633361  0.912650977333  0.906702031339 

f multiplier to get -3 dB at fc, for N = 
 1: 1.0*           
 2: 1.27201964951  
 3: 1.32267579991*  1.44761713315  
 4: 1.60335751622   1.43017155999  
 5: 1.50231627145*  1.75537777664   1.5563471223   
 6: 1.9047076123    1.68916826762   1.60391912877  
 7: 1.68436817927*  2.04949090027   1.82241747886   1.71635604487  
 8: 2.18872623053   1.95319575902   1.8320926012    1.77846591177  
 9: 1.85660050123*  2.32233235836   2.08040543586   1.94786513423   1.87840422428  
10: 2.45062684305   2.20375262593   2.06220731793   1.98055310881   1.94270419166  
11: 2.01670147346*  2.57403662106   2.32327165002   2.17445328051   2.08306994025   2.03279787154  
12: 2.69298925084   2.43912611431   2.28431825401   2.18496722634   2.12472538477   2.09613322542  
13: 2.16608270526*  2.80787865058   2.55152585818   2.39170950692   2.28570254744   2.21724536226   2.178598197    
14: 2.91905714471   2.66069088948   2.49663434571   2.38497976939   2.30961462222   2.26265746534   2.24005716132  
15: 2.30637004989*  3.02683647605   2.76683540993   2.5991524698    2.48264509354   2.4013780964    2.34741064497   2.31646357396  
16: 3.13149167404   2.87016099416   2.69935018044   2.57862945683   2.49225505119   2.43227707449   2.39427710712   2.37582307687  
17: 2.43892718901*  3.23326555056   2.97085412104   2.7973260082    2.67291517463   2.58207391498   2.51687477182   2.47281641513   2.44729196328  
18: 3.33237300564   3.06908580184   2.89318259511   2.76551588399   2.67073340527   2.60094950474   2.55161764546   2.52001358804   2.50457164552  
19: 2.56484755551*  3.42900487079   3.16501220302   2.98702207363   2.85646430456   2.75817808906   2.68433211497   2.630358907     2.59345714553   2.57192605454  
20: 3.52333123464   3.25877569704   3.07894353744   2.94580435024   2.84438325189   2.76691082498   2.70881411245   2.66724655259   2.64040228249   2.62723439989  
21: 2.68500843719*  3.61550427934   3.35050607023   3.16904164639   3.03358632774   2.92934454178   2.84861318802   2.78682554869   2.74110646014   2.70958138974   2.69109396043  
22: 3.70566068548   3.44032173223   3.2574059854    3.11986367838   3.01307175388   2.92939234605   2.86428726094   2.81483068055   2.77915465405   2.75596888377   2.74456638588  
23: 2.79958271812*  3.79392366765   3.52833084348   3.34412104851   3.204690112     3.09558498892   3.00922346183   2.94111672911   2.88826359835   2.84898013042   2.82125512753   2.80585968356  
24: 3.88040469682   3.61463243697   3.4292654707    3.28812274966   3.17689762788   3.08812364257   3.01720732972   2.96140104561   2.91862858495   2.88729479473   2.8674198668    2.8570800015   
25: 2.91440986162*  3.96520496584   3.69931726336   3.51291368484   3.37021124774   3.25705322752   3.16605475731   3.09257032034   3.03412738742   2.98814254637   2.9529987927    2.9314185899    2.91231068194  

bessel_calculate.py

# -*- coding: utf-8 -*-
"""
Calculate the f multipler and Q values for implementing Nth-order Bessel
filters as biquad second-order sections.

Based on description at http://freeverb3.sourceforge.net/iir_filter.shtml
 
For highpass filter, use fc/fmultiplier instead of fc*fmultiplier

Originally I back-calculated from the denominators produced by the bessel()
filter design tool, which stops at order 25.

Then I made a function to output Bessel polynomials directly, which can be 
used to calculate Q, but not f.  (TODO: Create real Bessel filters directly 
and calculate f.)  The two methods disagree with each other somewhat above 
8th order.  I'm not sure which is more accurate.

Also, these are bilinear-transformed, so they're only good below fs/4
(https://gist.github.com/endolith/4964212)

Created on Mon Feb 18 21:34:15 2013
 
"""

from __future__ import division
from numpy import roots, reshape, sqrt, poly, polyval
from scipy.misc import factorial
from scipy.signal import bessel
from scipy.optimize import brentq
from sos import cplxpair # https://gist.github.com/endolith/4525003

def bessel_poly(n, reverse=False):
    """Return the Bessel polynomial of degree `n`
    
    If `reverse` is true, a reverse Bessel polynomial is output.

    Output is a list of coefficients:
    [1]                   = 1
    [1,  1]               = 1*s   +  1
    [1,  3,  3]           = 1*s^2 +  3*s   +  3
    [1,  6, 15, 15]       = 1*s^3 +  6*s^2 + 15*s   +  15
    [1, 10, 45, 105, 105] = 1*s^4 + 10*s^3 + 45*s^2 + 105*s + 105
    etc.

    Output is a Python list of arbitrary precision long ints, so n is only 
    limited by your hardware's memory.

    Sequence is http://oeis.org/A001498 , and output can be confirmed to 
    match http://oeis.org/A001498/b001498.txt :
    
    i = 0
    for n in xrange(51):
        for x in bessel_poly(n, reverse=True):
            print i, x
            i += 1
    
    """
    out = []
    for k in xrange(n + 1):
        num = factorial(2*n - k, exact=True)
        den = 2**(n - k) * factorial(k, exact=True) * factorial(n - k, exact=True)
        out.append(num // den)

    if reverse:
        return list(reversed(out))
    else:
        return out


"""
Original:
    
1: ---
2: 0.58
3: --- 0.69
4: 0.81 0.52
5: ---- 0.92 0.56
6: 1.02 0.61 0.51
"""
print "Q for N = "
for n in xrange(1, 9):
    print str(n).rjust(2) + ":",
#    b, a = bessel(n, 1, analog=True)
    a = bessel_poly(n, reverse=True)
    p = cplxpair(roots(a)) # Poles, sorted into conjugate pairs
    if n % 2:
        # Odd-order, has a 1st-order stage
        print '-' * 14, # 1st-order stages don't have a Q
        # Remove 1st-order stage (single real pole at the end)
        p = reshape(p[:-1], (-1, 2))
    else:
        # Even-order, is all 2nd-order stages
        p = reshape(p, (-1, 2))

    for section in reversed(p):
        a = poly(section) # Convert back into a polynomial
        """
        Polynomial is
             s^2 + wo/Q*s + wo^2 = 
        a[0]*s^2 + a[1]*s + a[2]
        so Q is:
        """
        print str(sqrt(a[2])/a[1]).ljust(14),
    print
print


"""
The f requires two steps.  First calculate the f multiplier for each biquad
that produces a normalized Bessel filter (normalized so the asymptotes match 
a Butterworth)

With these settings, an LPF and HPF have the same phase curve vs frequency 
("phase-matched")

Numbers with asterisks are 1st-order filters
"""
print "f multiplier to get same asymptotes as Butterworth (LPF and HPF phase-matched), for N = "
for n in xrange(1, 9):
    print str(n).rjust(2) + ":",
    b, a = bessel(n, 1, analog=True)
    p = cplxpair(roots(a)) # Poles, sorted into conjugate pairs
    if n % 2:
        # Odd-order, has a 1st-order stage
        print (str(abs(p[-1])) + '*').ljust(14), # 1st-order stage
        # Remove 1st-order stage (single real pole at the end)
        p = reshape(p[:-1], (-1, 2))
    else:
        # Even-order, is all 2nd-order stages
        p = reshape(p, (-1, 2))

    for section in reversed(p):
        a = poly(section) # Convert back into a polynomial
        """
        Polynomial is
             s^2 + wo/Q*s + wo^2 = 
        a[0]*s^2 + a[1]*s + a[2]
        so wo is sqrt(a[2]):
        """
        print str(sqrt(a[2])).ljust(15),
    print
print


"""
Second, measure the point at which the frequency response of the 
normalized filter = -3 dB, and calculate the frequency multiplier to 
shift it so that the -3 dB point is at the desired frequency.

This then matches the original values:

1: 1.00
2: 1.27
3: 1.32 1.45
4: 1.60 1.43
5: 1.50 1.76 1.56
6: 1.90 1.69 1.60
"""
print "f multiplier to get -3 dB at fc, for N = "
for n in xrange(1, 9):
    print str(n).rjust(2) + ":",
    b, a = bessel(n, 1, analog=True)
    p = cplxpair(roots(a)) # Poles, sorted into conjugate pairs

    # Measure frequency at which magnitude response = -3 dB
    def H(w):
        """Output 0 when magnitude of frequency response is -3 dB"""
        # From scipy.signal.freqs:
        s = 1j * w
        return abs(polyval(b, s) / polyval(a, s)) - 1/sqrt(2)
    
    w = brentq(H, 0, 5)
    
    # Invert to get frequency multiplier
    m = 1/w
    

    if n % 2:
        # Odd-order, has a 1st-order stage
        print (str(m*abs(p[-1])) + '*').ljust(15), # 1st-order stage
        # Remove 1st-order stage (single real pole at the end)
        p = reshape(p[:-1], (-1, 2))
    else:
        # Even-order, is all 2nd-order stages
        p = reshape(p, (-1, 2))

    for section in reversed(p):
        a = poly(section) # Convert back into a polynomial
        """
        Polynomial is
             s^2 + wo/Q*s + wo^2 = 
        a[0]*s^2 + a[1]*s + a[2]
        so wo is sqrt(a[2])

        then multiply by m so it's -3 dB
        """
        print str(m*sqrt(a[2])).ljust(15),
    print
print