Skip to content

Instantly share code, notes, and snippets.

@endolith

endolith/All values.txt

Last active February 8, 2019 06:23
Show Gist options
  • Save endolith/4982787 to your computer and use it in GitHub Desktop.
Save endolith/4982787 to your computer and use it in GitHub Desktop.
Download ZIP
f and Q values for implementing Bessel filter as second-order sections (and Bessel polynomials in Python)
Raw
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
Raw
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment