ベッセル・フィルタ (Bessel filter) の計算例と設計表


Takayuk HOSODA
Apr. 11, 2024

English edition is here.

前書き

今は西暦2024年の春、Generative Pre-trained Transformer も桜も花盛りだというのに、 技術雑誌の付録の手帳にアクティブフィルタの設計表が載っていた。 花も GPT も見るもので使うものじゃないというのは分かるけど、 その数表のベッセル・フィルタ^[1][2]の設計値が「ちょっと」違う気がして気になった。

そんなところにセレンディピティ的なのを発揮してもしょうがないけど、気になってしまって少し調べてみた。 数表の出所は辿れなかったけど、カットオフ周波数を √0.5 の代わりに、 本末転倒感のある^(*1)厳密な -3 dB で計算した値と比較すると掲載値と 5 桁程度合うみたい。 5 桁程っていうのも、日本の近代回路網理論の黎明期が 1950 年代だと考えると、 元々の計算が利用出来る計算環境的に単精度とかが理由だったりするのだろうか。

実際の回路設計ではベッセル・フィルタを採用することも少ないし、 仮にあったとしてもさまざまな設計用アプリケーションが使える今となっては、 わざわざ数表から設計することも稀だとは思うものの、 アナログ信号とディジタル信号処理の狭間や、扱う信号の高速化にともなうエラーレートや EMC の観点からも、 ベッセル・フィルタそのものではなくても、 数学的な理解や本質的にベッセル・フィルタを必要とする場合もある。 そんな理由で、これをきっかけとして改めてベッセル・フィルタの計算例と計算プログラムと設計表をまとめておこうと思う。

正規化ベッセル・フィルタの計算例

3次ベッセル・フィルタの伝達関数 TBF(s, 3) は、時間遅れの伝達関数

の連分数展開による近似 ^[2][3] (i.e. 0 : 3 次の Padé 近似) から次のようになります。

TBF(s, 3) の振幅の絶対値が 1 / √2 になるカットオフ周波数 f_c は次のように計算されます。

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

ここで、solve(f(s), s) は、 関数 f(s) を s について解く^(*2)ことを意味します。

f_c で正規化した 3次ベッセル・フィルタの正規化伝達関数 TBFN(s, 3) は近似的に次のようになります。

TBFN(s, 3) の極 p は分母を解いて次のように求められます。

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

p のひとつの極 p に対する 固有振動数 ω とクォリティファクタ Q は次のように定義されます。

従って、2段構成の 3次ベッセル・フィルターの正規化設計パラメータは次のように計算されます。

Appendix

正規化ベッセル・フィルタの振幅特性 (1〜12 次)

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

正規化ベッセル・フィルタの SPICE3 シミュレーション結果 (1〜12 次)

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.

ベッセル・フィルタの伝達関数 (1〜12 次)

ベッセル・フィルタのカットオフ周波数 (1〜12 次)

正規化ベッセル・フィルタの極 (2〜12 次)

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

脚注

*1 : 電力が 1/2 になる周波数がカットオフ周波数であって 1/√2 ≃ -3 dB (i.e. 10-3/20 / (1/√2) - 1 ≃ 0.0011865297 ). であるため。無線工学的には 0.1 % 程の誤差は無視できるが数表となると話は別である。

*2 : 方程式の解法にはニュートン・ラフソン法 (Newton-Raphson method) や オストロフスキーの方法 (Ostrowski's method) ^[4][5] 、高次代数方程式の解法には DKA (Durand-Kerner-Aberth) 法 ^{[4] [6] [9]} などが利用される。

ソースコード

bf_dka -- Durand-Kerner-Aberth 法による逆ベッセル多項式の求根

bf_dka.c ◆ 数式のクリックで展開／折りたたみ

/* bf_dka -- A solver for the reverse Bessel polynomials using the Durand-Kerner-Aberth method.
 * Rev. 1.13 (extended precision, minimal edition) (May 13, 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>
#include "ost.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.13 (extended precision, minimal edition) (May 13, 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);
}

/* 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) {
            fprintf(stdout, "# %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) {
            fprintf(stdout, "# %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                  err;
    int                  n;
    #ifdef __FreeBSD__
        fpsetprec(FP_PE);
        fpsetround(FP_RN);
    #endif
    if (argc == 1) {
        fprintf(stdout, "%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) || (err = parse(argc, argv, &bf))) {
        fprintf(stdout, "%s -- %s\n  %s  %s\n  Usage : %s [-normalize] -order=n | n = 0..17\n", myname, descr, version, copyright, myname);
        return err;
    } 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() */
        ostl(1.0L, 5.0L, magfc, LDBL_EPSILON * 2, MAX_ITTR, &fc); /* solve for cut-off frequency */
        if (flags & F_NORMALIZE) normalizePoly(a, fc, n);         /* normalize the polynomial coefficients */
        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.13.tar.gz — 教育目的のための bf_dka の最小バージョンの C ソース コードです。 *BSD または Linux 環境でコンパイルできます。

結論

アナログフィルタ関係の教科書や実務書をいくつか紐解いてみたが、 計算の根拠付きで "正確な" 8桁以上の数表が載っているものがうまく見つけられなかった。 今は優れた数式処理ソフトウェアがいくつもあるが、改めて C99 の 64-bit 整数と倍精度あるいは拡張倍精度の複素数でどこまでいけるかやってみた。 数式処理ソフトとかと違って C で書くとアルゴリズムとか計算機限界に合わせた最適化とか全てが手の内に残るので他にも活用出来て嬉しい。 x86 拡張倍精度で DKA (Durand-Kerner-Aberth) 法を用いた上記プログラムでの計算結果は、 数式処理ソフトウェアでの 25 桁の多倍長計算の基準値との比較では、 12 次までは 14 桁、16 次までは 12 桁有効で、17 次では 11 桁まで有効だった (fig.5 参照)。 ここには未掲載の倍精度版での計算結果は 12 次までは 11 桁、16 次までは 9 桁有効で、17 次では 8 桁まで有効だった (fig.6 参照)。

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


fig.6 Order V.S. Maximum error and Itteraions until convergence (Ehrilich-Aberth and Durand-Kerner method, double precision)

追補

正規化ベッセル・フィルタの設計表 (2〜17 次)

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

ベッセル・フィルタのシミュレーション回路例 (3 次, VCVS 型, LTspice 用)

 bessel3-ltspice.zip

参考文献

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. 勝間田 明男 "ベッセルディジタルフィルタの自動設計について" 験震時報第56巻, 1993, p. 17-34
4. 戸川隼人 "科学技術計算ハンドブック", サイエンス社, 東京, 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. 山本 哲朗・古金 卯太郎・野倉 久美 "代数方程式を解く Durand-Kerner 法と Aberth 法", 情報処理　Vol.18 No.6 June 1977, p. 566-571
7. 白木 恒雄・寺野 隆雄・中村 秀治・栗原 千鶴子 "DKA法 による超越方程式の解法とその適用例", 土木学会論文集 第416号/1-13, 1990, p. 295-302
8. 名取 亮 編 長谷川 秀彦・櫻井哲也・檜山澄子・他 共著, "数値計算法", オーム社, 東京, Japan, 1998.9, p.109-113, ISBN 4-274-13153-X
9. 小澤 一文, "Durand-Kerner 法の効率的な初期値の簡単な設定法", 日本応用数理学会論文誌 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
11. Dunster, T. & Gil, Amparo & Ruiz-Antolin, Diego & Segura, Javier. (2021). Computation of the reverse generalized Bessel polynomials and their zeros. Computational and Mathematical Methods. 3. 10.1002/cmm4.1198.
12. Srivastava, Hari. (2023). An Introductory Overview of Bessel Polynomials, the Generalized Bessel Polynomials and the q-Bessel Polynomials. Symmetry. 15. 822. 10.3390/sym15040822.

関連項目

• オストロフスキー法による求根 (Solve by Ostrovski's method)
• パデ近似 (Padé approximant)

外部リンク

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