An example of Bessel filter calculation and
Inverse Bessel polynomial roots by the DKA (Durand-Kerner-Aberth) method


Takayuk HOSODA
Apr. 11, 2024

Preface

It is now spring 2024, and Generative Pre-trained Transformers and cherry blossoms are in full bloom, I found a design table of active filters in the handbook of a technical magazine supplement. I know that flowers and GPTs are for viewing and not for use, but I was concerned because the design values for the Bessel filter ^[1][2] in the numerical table seemed to be "slightly" different.

I don't want to be serendipitous about it, but I was curious and did a bit of research. I couldn't trace the source of the table, but it seems to match the published values by about five digits when compared with the values calculated with the strict -3 dB cut-off frequency instead of √0.5, which seems like a total oversimplification ^(*1) The reason for about five digits is that the dawn of modern network theory in Japan was in the 1950s, and the original calculations may have been due to the single precision of the available computing environment.

In actual circuit design, pure Bessel filters are rarely employed, and even if they were, it is rare to design from a table of numbers now that a variety of design applications are available. However, there are also cases where a Bessel filter is required in terms of mathematical understanding and intrinsic nature, even if it is not a Bessel filter itself, from the point of view of error rates and EMC. For these reasons, this is a good opportunity to summarise the calculation examples, calculation programs and design tables for Bessel filters once again.

Example calculation of a normalised Bessel filter

The transfer function of a third-order Bessel filter TBF(s, 3) can be approximated by a serial fraction expansion of the time delay transfer function ^[2][3] as follows

The cut-off frequency f_c at which the absolute value of the amplitude of TBF(s, 3) is 1 / √2 is calculated as follows

f_c = solve(abs(TBF(s, 3), s) - sqrt(0.5) → f_c ≃ 1.75567236868121064920

Where solve(f(s), s) means solve ^(*2) the function f(s) for s.

The normalised transfer function TBFN(s, 3) of the third-order Bessel filter normalised by f_c approximately becomes

The poles p of TBFN(s, 3) are obtained by solving the denominator as follows.

solve(15 / TBFN(s, 3), s) → p ≃ {-1.322675800, -1.047409161 ± i 0.9992644365 }

The natural frequency ω and quality factor Q for one pole p of p are defined as follows

Therefore, the normalised design parameters for a third-order Bessel filter in a two-stage configuration are calculated as follows

Examples of Bessel filter characteristics

Frequency versus amplitude characteristics of the normalised Bessel filters (1st to 12th order)

fig.1 Frequency response of the normalized Bessel filters (1st - 12th order)

Results of SPICE3 simulations of the normalised Bessel filters (1st to 12th order)

fig.2 Group delay characteristics of the normalized Bessel filters (1st to 12th order)
AC simulation results of the normalised Bessel filter of active filters with ideal op-amps using SPICE3.


fig.3 Transcient response of the normalized Bessel filters (1st to 12th order)
TRAN simulation results of the normalised Bessel filter of active filters with ideal op-amps using SPICE3.

Bessel Filter Specifications

The transfer function of the Bessel filters (1st to 12th order)

The cut-off frequency of the Bessel filters (1st to 12th order)

Poles of the normalised Bessel filter (2nd to 12th order)

fig.4 Poles of the normalized Bessel filters (1st - 16th order)
Poles are connected by lines for each group.

Roots of inverse Bessel polynomials

Source code (C99)

bf_dka -- A solver for reverse Bessel polynomials using the Durand-Kerner-Aberth method.

bf_dka.c ◆ click to unfold/fold

/* bf_dka -- A solver for the reverse Bessel polynomials using the Durand-Kerner-Aberth method.
 * Rev. 1.21 (extended precision, minimal edition) (May 31, 2024) (c) 2024, Takayuki HOSODA
 * The calculation results are valid for 14-digits upto 12-th order, 12-digits upto 16-th order
 * and 11-digis at 17-th order. C99 (e.g. gcc 4.2 or later) is required for compilation.
 * http://www.finetune.co.jp/~lyuka/technote/herb/bessel.html
 * DISCLAIMER : This program is free software. It is provided 'as is', without any warranty,
 * express or implied, under the extent permitted by applicable law. You are free
 * to redistribute and/or modify it.
 * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
 * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 */
#include <float.h>
#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef __FreeBSD__
#include <floatingpoint.h>
#endif
#ifndef M_PI
#define M_PI            3.14159265358979323846
#endif
#ifndef M_PI_2
#define M_PI_2          1.57079632679489661923
#endif
#ifndef LDBL_EPSILON
#define LDBL_EPSILON    1.0842021724855044340E-19L
#endif
#define DECELERATE    ((LDBL_EPSILON) * (long double)(1LL << 32))
#define THRESHOLD     ((LDBL_EPSILON) * (long double)(1LL << 27))
#define MAXN           17
#define MAX_ITTR      100
#define F_NORMALIZE     1
#define E_HELP          1
#define E_OPTION        2
#define E_DOMAIN        3

typedef struct fltr {
    int     order;
}   filter;

/* Though omitted for simplicity, the function prototypes should better be written here. */

int    flags = 0;
int    order;
char * myname = "bf_dka_mini";
char * descr = "A root-finding program for Bessel polynomials using the DKA method.";
char * version = "Rev. 1.21 (extended precision, minimal edition) (May 31, 2024)";
char * copyright = "(c) 2024, Takayuki HOSODA";
char * disclaim  = "is free software with ABSOLUTELY NO WARRANTY.";
long double coef[MAXN + 1];  /* global work variables to be used by ostl() */

/* complex polynomial */
long double complex poly(long double complex z, long double *c, int n) {
    long double complex p = 0.0;
    for (int i = 0; i <= n; i++) p = z * p + c[i];
    return p;
}

/* normalize coefficients of the polynomial */
int normalizePoly(long double *a, long double h, int n) {
    long double hi = h;
    for (int i = n; i >= 0; i--, hi *= h) a[i] *= hi;
    return n;
}

/* the function to be solved by ostl() */
long double magfc(long double x) {
    extern long double coef[];
    extern int    order;
    return cabsl(coef[order] / poly(x * I, coef, order)) - sqrtl(0.5L);
}

/* ostl - Solve the equation f(x) = 0 to within desired error value for x.
 * Based on ost Rev.2.3 (Mar. 3, 2010) (c) 2010 Takayuki HOSODA
 * http://www.finetune.co.jp/~lyuka/technote/fract/ost.html
 */
int ostl(long double x1, long double x2, long double (*func)(long double), long double eps, int ittr, long double *x) {
    long double     a, b, c, d, e, f, t;
    int             i;
    if (func == NULL) return 10;
    if (x == NULL) return 9;
    if (0 >= ittr) { return 8; }
    a = x1; b = x2;
    b = (a == b) ? a == 0 ? 0.01 : a * 1.01 : b;
    c = (a + b) * 0.5;
    if (0.0 == (d = (*func)(a))) { *x = a; return 0; }
    if (0.0 == (e = (*func)(b))) { *x = b; return 0; }
    if (0.0 == (f = (*func)(c))) { *x = c; return 0; }
    if (0.0 == (f - e) * (f - d)) {
        c = (a + 3 * b) * 0.25;
        if (0.0 == (f = (*func)(c))) { *x = c; return 0; }
        if (0.0 == (f - e) * (f - d)) { *x = c; return 7; }
    }
    for (i = 0; i < ittr; i++) {
        if (0.0 == f) { *x = c; return 0; }
        if (fabsl(f) <= eps)  { *x = c; return 1; }
        if (0.0 == (c - a)) { *x = c; return 2; }
        if (0 == (f - e)) { *x = c; return 3; }
        if ((fabsl(f) > fabsl(e)) && (fabsl(e) <= sqrtl(eps))) { *x = b; return 4; }
        /* Ostrowski's method :
        t &=& \frac{x_\mathrm{n+2} - x_\mathrm{n+1}}{x_\mathrm{n+2} - x_\mathrm{n}}  
              \frac{f(x_\mathrm{n+2}) - f(x_\mathrm{n})}{f(x_\mathrm{n+2}) - f(x_\mathrm{n+1})} 
              \frac{f(x_\mathrm{n+1})}{f(x_\mathrm{n})}\\
        x_\mathrm{n+3} &=& \frac{x_\mathrm{n+1} - tx_\mathrm{n}}{1 - t} */
        t = e * (f - d) * (c - b);
        t /= d * (f - e) * (c - a);
        if (t == 1.0) { *x = b; return 6; }
        t = (t * a - b) / (t - 1.0);
        a = b; b = c; c = t;
        d = e; e = f; f = (*func)(c);
    }
    *x = c; return 5;
}

/* EhrilichAberth method :
\widehat{x_i} &=& x_i - \frac{{p_n(x_i)}}{{{p'_n(x_i)} - {p_n(x_i)} \displaystyle\sum_{\substack{j = 0\\j  \neq i}}^ {n-1}  \frac{1}{x_i - x_j}}} */
long double complex diffEhrilichAberth(long double complex *x, long double *a, int k, int n) {
    long double complex s = 0.0;
    long double complex p = 0.0;
    long double complex d = 0.0;
    long double complex z = x[k];
    for (int i = 0; i < n; i++) {
        p = z * p + a[i];
        d = z * d + a[i] * (n - i); /* d = p'_n(z_i) / p_n(z_i) */
    }
    p = a[n] + z * p;
    for (int i = 0; i < n; i++)
        if (i != k) s += 1.0 / (z - x[i]); /* s = \sum_{j = 0, j \neq i}^{n - 1}(1 / (z_i - z_j) */
    z = (p == 0.0) ? 0.0 : p / (d - p * s); /* workarround for complex division of zero by infinity */
    return z;
}

/* maximum */
long double maxfabsl(long double *a, int n) {
    long double max = 0;
    for (int i = 0; i <= n; i++) if (max < fabsl(a[i])) max = fabsl(a[i]);
    return max;
}

/* magnitude maximum */
long double maxcabsl(long double complex *c, int n) {
    long double max = 0;
    for (int i = 0; i <= n; i++) if (max < cabsl(c[i])) max = cabsl(c[i]);
    return max;
}

/* set Aberth's initial values for c[] to x[] and returns its radious r0 */
/* valid only for the reverse Bessel polynomials (upto at least n = 17) */
long double setAberth(long double complex *x, long double *a, int n) {
    long double complex g;
    long double r0;
    g = - a[1] / n; /* g = -\frac{a_1}{n} */
    r0 = fabsl(g) + pow(maxfabsl(a, n), 1.0 / n); /* r_0 = |g| + \sqrt[\displaystyle  n]{\max_{0\le k \le n}|a_k|} */
    for (int k = 0; k < n; k++)
        x[k] = g + r0 * cexp(I * (2 * M_PI * k + M_PI_2 ) / n);
        /* x_k = g + r_0\exp\left\{\left(\frac{2\pi k}{n}+\frac{\pi}{2n}}\right)\mathrm{i}\right\}\;|\;0\le k<n,\;\mathrm{i}^2=-1 */
    return r0;
}

/* integer reduce */
long long reduce(long long *a, long long *b) {
    long long t, p, q;
    p = *a;
    q = *b;
    do {
        if (p <= q) {
            t = p; p = q; q = t;
        }
        p %= q;
    } while (p);
    *a /= q;
    *b /= q;
    return q;
}

/* factorial */
long long fact(long long ld) {
    long long m;
    if (ld <= 0) return 1LL;
    m = ld--;
    while(ld > 1LL)
        m *= ld--;
    return m;
}

/* partial factorial */
long long pfact(long long j, long long k) {
    long long m;
    if (j <= 0) return 1LL;
    m = j--;
    while(j > k)
        m *= j--;
    return m;
}

/* Bessel polynomial
\displaystyle{a_{i}={\frac {(2n-i)!}{2^{n-i}i!(n-i)!}}\quad i=0..n} */
long long rBessCoef(long double *a, int n) {
    long long coef = 1LL;
    long long p, q, r, b;
    for (int i = n; i >= 0; i--) {
        if ((2 * n - i) >= 22) {
            q = pfact((2 * n - i), 21); p = fact(i);
            r = pfact(21, (n - i));     b = 1LL << (n - i);
            reduce(&q, &p); reduce(&r, &b);
            reduce(&q, &b); reduce(&r, &p);
            coef = (q * r) / (p * b);
            a[n - i] = (long double)coef;
        } else {
            a[n - i] = (long double)(coef = (pfact((2LL * n - i), (n - i)) / (fact(i) * (1LL << (n - i)))));
        }
    }
    return coef;
}

/* the comparator function of qsort() */
int cCompar(const void * c1, const void * c2) {
    if (cimagl(*(long double complex *)c1) > cimagl(*(long double complex *)c2)) {
        return -1;
    } else if (cimagl(*(long double complex *)c1) < cimagl(*(long double complex *)c2)) {
        return 1;
    } else {
        return 0;
    }
}

/* polynomial solver by Durand-Aberth-Kerner (Ehrilich-Aberth) method */
long double solvePoly(long double complex *x, long double *a, int n) {
    long double ep;
    long double er = DBL_MAX;
    long double decelerate;
    long double complex xn[n];
    long double complex dx[n];
    for (int k = 1; k < MAX_ITTR; k++) {
        decelerate = er <= DECELERATE ? 0.25 : 1.0;                        /* Decelerate near convergence */
        ep = er;
        for (int i = 0; i < n; i++)
            xn[i] = x[i] - (dx[i] = decelerate * diffEhrilichAberth(x, a, i, n));   /* Ehrilich-Aberth method */
        er = maxcabsl(dx, n - 1);
        for (int i = 0; i < n; i++) x[i] = xn[i]; /* copy array */
        if (er < (DBL_EPSILON * (long double)(2 * n)) || ((ep <= er) && (er < THRESHOLD))) break; /* convergence check */
    }
    return er;
}

/* print the roots of the reverse Bessel polynomials solved, i.e. poles of the Bessel filters */
void printResults(long double complex *x, int n) {
    qsort(x, n, sizeof(long double complex), cCompar);       /* sort x[] in descending order of imaginary part */
    for (int i = 0; i < n; i++) {                            /* print results with Q and W */
        if (cimagl(x[i]) > DBL_EPSILON) {
            printf("# %sPole%d[%d] = %15.11Lf +-%15.11Lf * I (Q = %15.11Lf, w = %15.11Lf)\n",
                flags & F_NORMALIZE ? "n" : "", n, i, creall(x[i]), cimagl(x[i]),\
                cabsl(x[i]) / (-2 * creall(x[i])), cabsl(x[i]));
        } else if (cabsl(cimagl(x[i])) <= DBL_EPSILON) {
            printf("# %sPole%d[%d] = %15.11Lf                       (Q = %15.11Lf, w = %15.11Lf)\n",
                flags & F_NORMALIZE ? "n" : "", n, i, creall(x[i]),\
                cabsl(x[i]) / (-2 * creall(x[i])), cabsl(x[i]));
        }
    }
}

/* command-line parser */
int parse(int argc, char *argv[], filter *bf) {
    int rtnv = 0;
    bf->order = 0;
    for (int i = 1; i < argc; i++) {
        if (*argv[i] == '-') {
            switch (*++argv[i]) {
                case 'n':
                if (!strncmp(argv[i], "normalize", strlen("normalize"))) {
                    flags |= F_NORMALIZE;
                } else {
                    rtnv = E_OPTION;
                }
                break;
                case 'o':
                if (!strncmp(argv[i], "order=", strlen("order="))) {
                    sscanf(&argv[i][strlen("order=")], "%d", &bf->order);
                    if (bf->order <= 0) {
                        rtnv = E_DOMAIN;
                    } else if (bf->order > MAXN) {
                        rtnv = E_DOMAIN;
                    }
                } else {
                    rtnv = E_OPTION;
                }
                break;
            default:
                rtnv = E_OPTION;
                break;
            }
        }
    }
    if (bf->order == 0) {
        bf->order = MAXN;
    }
    return (rtnv);
}
int main(int argc, char *argv[]) {
    extern int           order;
    extern long double   coef[MAXN + 1];
    long double          a[MAXN + 1];      /* coefficients */
    long double complex  x[MAXN];          /* approximation values of roots */
    long double          fc;               /* cut-off frequency */
    filter               bf;
    int                  r;
    int                  n;
    #ifdef __FreeBSD__
        fpsetprec(FP_PE);
        fpsetround(FP_RN);
    #endif
    if (argc == 1) {
        printf("%s -- %s\n  %s  %s\n%s %s\n  Usage : %s [-normalize] -order=n | n = 0..17\n",
             myname, descr, version, copyright, myname, disclaim, myname);
        return E_HELP;
    }
    if ((argc == 1) || (r = parse(argc, argv, &bf))) {
        printf("%s -- %s\n  %s  %s\n  Usage : %s [-normalize] -order=n | n = 0..17\n", myname, descr, version, copyright, myname);
        return r;
    } else {
        order = n = bf.order;                                     /* set order */
        rBessCoef(a, n);                                          /* set reverse Bessel coefficients to a[0..n] */
        for (int i = 0; i <= n; i++) coef[i] = a[i];              /* copy array to global work variables to be used by ostl() */
        r = ostl(1.0L, 5.0L, magfc, LDBL_EPSILON, MAX_ITTR, &fc); /* solve for cut-off frequency */
        if (flags & F_NORMALIZE) normalizePoly(a, fc, n);         /* normalize the polynomial coefficients */
        printf("# fc = %.11Lg\n", flags & F_NORMALIZE ? 1 : fc);  /* print cut-off frequency */
        for (int i = n; i >= 0; i--) a[i] /= a[0];                /* make a[] monic */
        setAberth(x, a, n);                                       /* set initial values */
        solvePoly(x, a, n);                                       /* solve polynomial by DKA method */
        printResults(x, n);                                       /* print results */
    }
    return 0;
}

 bf_dka_xp_mini-1.21.tar.gz — C source code for a minimal version of bf_dka for educational purposes. Can be compiled under *BSD or Linux.

Conclusion

I tried to find some textbooks and practical books related to analogue filters, but... I haven't been able to find an "accurate" table of 8+ digit numbers with a basis for calculations. There are a number of excellent CAS (Computer Algebra Systems) available nowadays, but I tried again to see how far I could get with 64-bit integers and C99 double-precision complex or extended-precision complex numbers. Unlike CAS, when you write in C, you can keep all the algorithms and optimisations for the limits of the computer in your hands, so I'm glad to be able to use it for other purposes. The results of the above program using the Durand-Kerner-Aberth (DKA) method in x86 extended-precision were valid up to 14 digits at 12th order, 12 digits at 16th order and 8 digits at 17th order in comparison with the reference value which was calculated to 25 digits of precision in CAS (c.f. fig.5). For the double precision version not listed here, 11 digits were valid up to the 12th order, 9 digits up to the 16th order and 8 digits up to the 17th order (c.f. fig.6).

fig.5 Maximum error V.S. order (Ehrilich-Aberth and Durand-Kerner method, x86 extended precision)


fig.6 Maximum error V.S. order (Ehrilich-Aberth and Durand-Kerner method, double precision)

Note

*1 : By definition, the cut-off frequency is the frequency at which the transferred power is halved, since it is 1/√2 ≃ -3 dB (i.e. 10-3/20 / (1/√2) - 1 ≃ 0.0011865297 ). In radio engineering, an error of about 0.1 % is negligible, but in the case of numerial tables it is a different story.

*2 : The Newton-Raphson method, Ostrowski's method ^[4][5] and others are used for solving equations, DKA (Durand-Kerner-Aberth) method ^{[4] [6] [9]} and others are used for solving higher-order algebraic equations.

Appendix — Design table

Design table for normalised Bessel filters (2nd to 17th order)

 normalized-Bessel-filter-design-parameters-1.31.pdf

Appendix — Circuit example

Bessel filter simulation circuit example (3rd order, VCVS type, for LTspice)

 bessel3-ltspice.zip

References

1. G.N.Watson. A TREATISE ON THE THEORY OF BESSEL FUNCTIONS. SECOND EDITION. Cambridge University Press, 1966, 815p
2. Thomson, W.E. "Delay networks having maximally flat frequency characteristics", Proceedings of the IEE - Part III: Radio and Communication Engineering, 1949, 96, (44), p. 487-490
3. Akio Katsumada "About the automatic design of Bessel digital filters" Seismic Test Report Vol. 56, 1993, p. 17-34
4. Hayato Togawa "Scientific and Technical Computing Handbook", Science Publishing, Tokyo, 1992, ISBN 4-7819-0666-4
5. .M.Ostrowskis. SOLUTION OF EQUATIONS IN EUCLIDEAN AND BANACH SPACES. THIRD EDITION. Academic Press, New Yoark and London, 1973
6. Tetsuro Yamamoto, Utaro Furukane, Kumi Nokura "Durand-Kerner method and Aberth method for solving algebraic equations", Information Processing Vol.18 No.6 June 1977, p. 566-571
7. Tsuneo Shiraki, Takao Terano, Hideharu Nakamura, and Chizuko Kurihara, "Solution of transcendental equations using the DKA method and its application examples," Proceedings of the Japan Society of Civil Engineers, No. 416/1-13, 1990, p. 295-302
8. Edited by Ryo Natori, co-authored by Hidehiko Hasegawa, Tetsuya Sakurai, Sumiko Hiyama, et al., "Numerical Calculation Method", Ohmsha, Tokyo, Japan, 1998.9, p.109-113, ISBN 4-274-13153-X
9. Kazufumi Ozawa, "Simple method for setting efficient initial values ??of the Durand-Kerner method", Transactions of the Japan Society of Applied Mathematics Vol.3, No.4, 1993, p.173-186
10. Anatol I. Zverev, "Handbook of Filter Synthesis", John Wiley and Sons, Inc, New York London Sydney, 1967, p. 61-71

See Also

• Solve by Ostrovski's method

External Links

• BesselFilterModel — WOLFRAM
• Bessel analog filter design — MathWorks
• Gershgorin Circle Theorem — Wolfram MathWorld
• Bessel filter — Wikipedia
• Pade approximant — Wikipedia
• Spice — a circuit simulator / UC Berkeley
• LTspice — free SPICE simulator software / Analog Devices / Linear Technology
• gnuplot — graphing utility
• GCC — the GNU Compiler Collection