Python design of Bessel filters of any order, with 3 common frequency normalizations

readme.md

This has been merged into SciPy now:

http://scipy.github.io/devdocs/generated/scipy.signal.besselap.html

http://scipy.github.io/devdocs/generated/scipy.signal.bessel.html

bessel_filters.py

# -*- coding: utf-8 -*-
"""
Generate Bessel filters directly, using mpmath

References:

http://en.wikipedia.org/wiki/Bessel_filter#Example

Thomson, W.E., "Delay Networks having Maximally Flat Frequency
    Characteristics", Proceedings of the Institution of Electrical Engineers,
    Part III, November 1949, Vol. 96, No. 44, pp. 487–490.

Bond, C.R., Bessel Filters Polynomials, Poles and Circuit Elements, 2003
www.crbond.com/papers/bsf2.pdf

Bond, C.R., Bessel Filter Constants
www.crbond.com/papers/bsf.pdf

Miller, Ray - A Bessel Filter Crossover, and Its Relation to Others
http://www.rane.com/pdf/ranenotes/Bessel_Filter_Crossover_&_its_Relation_to_Others.pdf

"""

from __future__ import division, print_function
from numpy import asarray
import numpy as np
from scipy.misc import factorial

try:
    from mpmath import mp
    mpmath_available = True
except ImportError:
    mpmath_available = False


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


def normfactor(a):
    """
    Frequency shift to apply to delay-normalized filter such that -3 dB point
    is at 1 rad/sec.

    First 10 values are listed in
    "Bessel Filters Polynomials, Poles and Circuit Elements 2003, C. Bond."
    """
    def G(w):
        """
        Gain of filter
        """
        return abs(a[-1]/mp.polyval(a, 1j*w))

    def cutoff(w):
        """
        When gain = -3 dB, return 0
        """
        return G(w) - 1/mp.sqrt(2)

    return mp.findroot(cutoff, 1.5)


"""
Some values:

1: 1,
2: 1.3616541287161305209588868444454448432347103729901,
3: 1.7556723686812106492036418865538548622632538627204,
4: 2.1139176749042158430196891286338507228884592920213,
5: 2.4274107021526281376663058876113906340172394775656,
6: 2.70339506120292187639578951630254258339151785719,
7: 2.9517221470387227719289052270107816308339931866838,
8: 3.179617237510651330501649153152461962340386280172,
9: 3.3916931389116601011122541953570649508638576678231,
10: 3.5909805945691634827892977108436993654835433664576,
11: 3.7796074164396200928908480462973530897057513767288,
12: 3.9591508211442853155392550478018629152262858202629,
13: 4.1308254993835359808364524050964109423143665202459,
14: 4.2955934095336375649772937469630970087603802703228,
15: 4.4542330216243774943378240757836612531623748693155,
16: 4.6073854654726479174415758515725454199893109030358,
17: 4.7555865489611477270354413366716975695016018946715,
18: 4.8992896772844880074416221966614330534420251746993,
19: 5.0388826814882076058089530669453332343792794904578,
20: 5.1747004417427074232967413801802587948461890472997,
21: 5.3070345313609172742698274606247947796240516171721,
22: 5.4361407032500359996281001872606472534988925645919,
23: 5.5622447837878781965332371305558405356769995361022,
24: 5.6855473712959635219603072771407213415136770804137,
25: 5.8062276237754185419638409233382918905704274993629,
"""


def _roots(a):
    """
    Find the roots of a polynomial, using mpmath.polyroots if available,
    or numpy.roots if not
    """

    N = (len(a) - 1)//2  # Order of the filter

    if mpmath_available:
        # Overkill: "The user may have to manually set the working precision
        # higher than the desired accuracy for the result, possibly much
        # higher."
        mp.dps = 150
        """
        bessel:
        with dps 50 and maxsteps = 100, it can do 39 order filters

        with dps 150 and maxsteps = 1000, it can do 70 order filters


        TODO: How many digits are necessary for float equivalence?  Does it
        vary with order?

        with dps 50 and maxsteps = 100, it can do __ order filters

        with dps 150 and maxsteps = 1000, it can do __>=40__ order filters
        """
        p, err = mp.polyroots(a, maxsteps=1000, error=True)
        if err > 1e-32:
            raise ValueError("Filter cannot be accurately computed "
                             "for order %s" % N)
        p = asarray(p).astype(complex)
    else:
        p = np.roots(a)
        if N > 25:
            # Bessel and Legendre filters seem to fail above N = 25
            raise ValueError("Filter cannot be accurately computed "
                             "for order %s" % N)
    return p


def besselap(N, norm='phase'):
    """
    Return (z,p,k) zero, pole, gain for analog prototype of an Nth order
    Bessel filter, normalized to a frequency of 1.

    `norm` is the frequency normalization:

    'phase': (default)

        The filter is normalized such that the filter asymptotes are the same
        as a Butterworth filter of the same order with an angular (e.g. rad/s)
        cutoff frequency of 1.

        The phase reaches its midpoint at this cutoff frequency for both
        low-pass and high-pass filters, so this is the "phase-matched" case.

        This is Matlab's normalization

    'delay':

        The filter is normalized such that the group delay in the passband
        is 1 (e.g. 1 second).  This is the "natural" type you get when you
        solve Bessel polynomials

    'mag':

        The filter is normalized such that the gain magnitude is -3 dB at
        frequency 1.  This is called "frequency normalization" by Bond

    """
    if N == 0:
        return asarray([]), asarray([]), 1

    # Find delay-normalized Bessel filter poles
    a = bessel_poly(N, reverse=True)
    p = _roots(a)

    if norm == 'delay':
        # Normalized for group delay of 1
        k = a[-1]
    elif norm == 'phase':
        # Phase-matched (1/2 max phase shift at 1 rad/sec)
        # Asymptotes are same as Butterworth filter
        p *= 1 / a[-1]**(1/N)
        k = 1
    elif norm == 'mag':
        # -3 dB magnitude point is at 1 rad/sec
        p *= 1 / normfactor(a)
        k = float(normfactor(a)**-N * a[-1])
    else:
        raise ValueError('normalization not understood')

    z = []
    p = p.astype(complex)

    return asarray(z), asarray(p), k


def tests():
    from scipy import signal
    from numpy.testing import (assert_allclose, assert_array_equal,
                               assert_array_almost_equal)
    from sos_stuff import cplxreal

    wikipedia_bessels = {
        0: [1],
        1: [1, 1],
        2: [3, 3, 1],
        3: [15, 15, 6, 1],
        4: [105, 105, 45, 10, 1],
        5: [945, 945, 420, 105, 15, 1]
        }

    for N in range(6):
        assert_array_equal(wikipedia_bessels[N], bessel_poly(N))
        assert_array_equal(wikipedia_bessels[N][::-1],
                           bessel_poly(N, reverse=True))

    a = bessel_poly(2, reverse=True)
    # Compare with symbolic result
    assert_allclose(float(normfactor(a)),
                    (np.sqrt(3)*np.sqrt(np.sqrt(5)-1)) / np.sqrt(2))

    # Values from Bond PDF
    Bond = {2: 1.36165412871613,
            3: 1.75567236868121,
            4: 2.11391767490422,
            5: 2.42741070215263,
            6: 2.70339506120292,
            7: 2.95172214703872,
            8: 3.17961723751065,
            9: 3.39169313891166,
            10: 3.59098059456916,
            }

    for N in range(2, 11):
        a = bessel_poly(N, reverse=True)
        assert_allclose(float(normfactor(a)), Bond[N])

    # 6th order
    bond_b = 10395
    bond_a = [1, 21, 210, 1260, 4725, 10395, 10395]
    b, a = zpk2tf(*besselap(6, 'delay'))
    assert_allclose(bond_b, b)
    assert_allclose(bond_a, a)

    bond_poles = {
        1: [-1.0000000000],
        2: [-1.5000000000 + 0.8660254038j],
        3: [-1.8389073227 + 1.7543809598j, -2.3221853546],
        4: [-2.1037893972 + 2.6574180419j, -2.8962106028 + 0.8672341289j],
        5: [-2.3246743032 + 3.5710229203j, -3.3519563992 + 1.7426614162j,
            -3.6467385953],
        6: [-2.5159322478 + 4.4926729537j, -3.7357083563 + 2.6262723114j,
            -4.2483593959 + 0.8675096732j],
        7: [-2.6856768789 + 5.4206941307j, -4.0701391636 + 3.5171740477j,
            -4.7582905282 + 1.7392860611j, -4.9717868585],
        8: [-2.8389839489 + 6.3539112986j, -4.3682892172 + 4.4144425005j,
            -5.2048407906 + 2.6161751526j, -5.5878860433 + 0.8676144454j],
        9: [-2.9792607982 + 7.2914636883j, -4.6384398872 + 5.3172716754j,
            -5.6044218195 + 3.4981569179j, -6.1293679043 + 1.7378483835j,
            -6.2970191817],
        10: [-3.1089162336 + 8.2326994591j, -4.8862195669 + 6.2249854825j,
             -5.9675283286 + 4.3849471889j, -6.6152909655 + 2.6115679208j,
             -6.9220449054 + 0.8676651955j]
        }

    for N in range(1, 11):
        p1 = np.sort(bond_poles[N])
        p2 = np.sort(np.concatenate(cplxreal(besselap(N, 'delay')[1])))
        assert_array_almost_equal(p1, p2, decimal=10)

    bond_poles = {
        1: [-1.0000000000],
        2: [-1.1016013306 + 0.6360098248j],
        3: [-1.0474091610 + 0.9992644363j, -1.3226757999],
        4: [-0.9952087644 + 1.2571057395j, -1.3700678306 + 0.4102497175j],
        5: [-0.9576765486 + 1.4711243207j, -1.3808773259 + 0.7179095876j,
            -1.5023162714],
        6: [-0.9306565229 + 1.6618632689j, -1.3818580976 + 0.9714718907j,
            -1.5714904036 + 0.3208963742j],
        7: [-0.9098677806 + 1.8364513530j, -1.3789032168 + 1.1915667778j,
            -1.6120387662 + 0.5892445069j, -1.6843681793],
        8: [-0.8928697188 + 1.9983258436j, -1.3738412176 + 1.3883565759j,
            -1.6369394181 + 0.8227956251j, -1.7574084004 + 0.2728675751j],
        9: [-0.8783992762 + 2.1498005243j, -1.3675883098 + 1.5677337122j,
            -1.6523964846 + 1.0313895670j, -1.8071705350 + 0.5123837306j,
            -1.8566005012],
        10: [-0.8657569017 + 2.2926048310j, -1.3606922784 + 1.7335057427j,
             -1.6618102414 + 1.2211002186j, -1.8421962445 + 0.7272575978j,
             -1.9276196914 + 0.2416234710j]
        }

    for N in range(1, 11):
        p1 = np.sort(bond_poles[N])
        p2 = np.sort(np.concatenate(cplxreal(besselap(N, 'mag')[1])))
        assert_array_almost_equal(p1, p2, decimal=10)

    for N in range(26):
        assert_allclose(sorted(signal.besselap(N)[1]), sorted(besselap(N)[1]))

    # Rane note examples
    a = [1, 1, 1/3]
    b2, a2 = zpk2tf(*besselap(2, 'delay'))
    assert_allclose(a[::-1], a2/b2)

    a = [1, 1, 2/5, 1/15]
    b2, a2 = zpk2tf(*besselap(3, 'delay'))
    assert_allclose(a[::-1], a2/b2)

    a = [1, 1, 9/21, 2/21, 1/105]
    b2, a2 = zpk2tf(*besselap(4, 'delay'))
    assert_allclose(a[::-1], a2/b2)

    a = [1, np.sqrt(3), 1]
    b2, a2 = zpk2tf(*besselap(2, 'phase'))
    assert_allclose(a[::-1], a2/b2)

    # TODO: Why so inaccurate?
    a = [1, 2.481, 2.463, 1.018]
    b2, a2 = zpk2tf(*besselap(3, 'phase'))
    assert_array_almost_equal(a[::-1], a2/b2, decimal=1)

    # TODO: Why so inaccurate?
    a = [1, 3.240, 4.5, 3.240, 1.050]
    b2, a2 = zpk2tf(*besselap(4, 'phase'))
    assert_array_almost_equal(a[::-1], a2/b2, decimal=1)

    # Phase to -3 dB factors:
    N, scale = 2, 1.272
    scale2 = besselap(N, 'mag')[1] / besselap(N, 'phase')[1]
    assert_array_almost_equal(scale, scale2, decimal=3)

    # TODO: Why so inaccurate?
    N, scale = 3, 1.413
    scale2 = besselap(N, 'mag')[1] / besselap(N, 'phase')[1]
    assert_array_almost_equal(scale, scale2, decimal=2)

    # TODO: Why so inaccurate?
    N, scale = 4, 1.533
    scale2 = besselap(N, 'mag')[1] / besselap(N, 'phase')[1]
    assert_array_almost_equal(scale, scale2, decimal=1)


if __name__ == "__main__":
    """
    run in Examples of besselap??
    """

    from scipy.signal import freqs, zpk2tf
    import matplotlib.pyplot as plt

    for norm in ('delay', 'phase', 'mag'):

        plt.figure(norm)
        plt.suptitle(norm + ' normalized')

        for N in (1, 2, 3, 4):
            z, p, k = besselap(N, norm)
            b, a = zpk2tf(z, p, k)
            w, h = freqs(b, a, np.logspace(-3, 3, 1000))
            plt.subplot(3, 1, 1)
            plt.semilogx(w, 20*np.log10(h))

            plt.subplot(3, 1, 2)
            plt.semilogx(w[1:], -np.diff(np.unwrap(np.angle(h)))/np.diff(w))

            plt.subplot(3, 1, 3)
            plt.semilogx(w, np.unwrap(np.angle(h)))

        plt.subplot(3, 1, 1)
        plt.title('Magnitude')
        plt.ylabel('dB')
        plt.axvline(1, color='red')
        plt.ylim(-40, 3)
        plt.grid(True, color='0.7', linestyle='-', which='major')
        plt.grid(True, color='0.9', linestyle='-', which='minor')

        plt.subplot(3, 1, 2)
        plt.title('Group delay')
        plt.ylabel('seconds')
        plt.axvline(1, color='red')
        plt.ylim(0, 3.5)
        plt.grid(True, color='0.7', linestyle='-', which='major')
        plt.grid(True, color='0.9', linestyle='-', which='minor')

        plt.subplot(3, 1, 3)
        plt.title('Phase')
        plt.ylabel('radians')
        plt.axvline(1, color='red')
        plt.ylim(-2*np.pi, 0)
        plt.grid(True, color='0.7', linestyle='-', which='major')
        plt.grid(True, color='0.9', linestyle='-', which='minor')

    plt.figure('mag')
    plt.subplot(3, 1, 1)
    plt.axhline(-3.01, color='red', alpha=0.5)

    plt.figure('delay')
    plt.subplot(3, 1, 2)
    plt.axhline(1, color='red', alpha=0.5)

    plt.figure('phase')
    plt.subplot(3, 1, 3)
    plt.axhline(-np.pi, color='red', alpha=0.5)
    plt.axhline(-3*np.pi/4, color='red', alpha=0.5)
    plt.axhline(-np.pi/2, color='red', alpha=0.5)

Normalized delay.png

Normalized magnitude.png

Normalized phase.png