Bessel filter

In electronics and signal processing, a Bessel filter is a type of analog linear filter with a maximally flat group delay (i.e., maximally linear phase response), which preserves the wave shape of filtered signals in the passband. ^[1] Bessel filters are often used in audio crossover systems.

The filter's name is a reference to German mathematician Friedrich Bessel (1784–1846), who developed the mathematical theory on which the filter is based. The filters are also called Bessel–Thomson filters in recognition of W. E. Thomson, who worked out how to apply Bessel functions to filter design in 1949. ^[2]

The Bessel filter is very similar to the Gaussian filter, and tends towards the same shape as filter order increases. ^{[3] [4]} While the time-domain step response of the Gaussian filter has zero overshoot, ^[5] the Bessel filter has a small amount of overshoot, ^{[6] [7]} but still much less than other common frequency-domain filters, such as Butterworth filters. It has been noted that the impulse response of Bessel–Thomson filters tends towards a Gaussian as the order of the filter is increased. ^[3]

Compared to finite-order approximations of the Gaussian filter, the Bessel filter has better shaping factor, flatter phase delay, and flatter group delay than a Gaussian of the same order, although the Gaussian has lower time delay and zero overshoot. ^[8]

The transfer function

A plot of the gain and group delay for a fourth-order low-pass Bessel filter. Note that the transition from the passband to the stopband is much slower than for other filters, but the group delay is practically constant in the passband. The Bessel filter maximizes the flatness of the group delay curve at zero frequency.

A Bessel low-pass filter is characterized by its transfer function: ^[9]

${\displaystyle H(s)={\frac {\theta _{n}(0)}{\theta _{n}(s/\omega _{0})}}\,}$

where ${\displaystyle \theta _{n}(s)}$ is a reverse Bessel polynomial from which the filter gets its name and ${\displaystyle \omega _{0}}$ is a frequency chosen to give the desired cut-off frequency. The filter has a low-frequency group delay of ${\displaystyle 1/\omega _{0}}$. Since ${\displaystyle \theta _{n}(0)}$ is indeterminate by the definition of reverse Bessel polynomials, but is a removable singularity, it is defined that ${\displaystyle \theta _{n}(0)=\lim _{x\rightarrow 0}\theta _{n}(x)}$.

Bessel polynomials

The roots of the third-order Bessel polynomial are the poles of the filter transfer function in the ${\displaystyle s}$ plane, here plotted as crosses.

The transfer function of the Bessel filter is a rational function whose denominator is a reverse Bessel polynomial, such as the following:

${\displaystyle n=1:\quad s+1}$

${\displaystyle n=2:\quad s^{2}+3s+3}$

${\displaystyle n=3:\quad s^{3}+6s^{2}+15s+15}$

${\displaystyle n=4:\quad s^{4}+10s^{3}+45s^{2}+105s+105}$

${\displaystyle n=5:\quad s^{5}+15s^{4}+105s^{3}+420s^{2}+945s+945}$

The reverse Bessel polynomials are given by: ^[9]

${\displaystyle \theta _{n}(s)=\sum _{k=0}^{n}a_{k}s^{k},}$

where

${\displaystyle a_{k}={\frac {(2n-k)!}{2^{n-k}k!(n-k)!}}\quad k=0,1,\ldots ,n.}$

Setting the cutoff attenuation

There is no standard set attenuation value for Bessel filters, however, −3.0103 dB is a common choice. Some applications may use a higher or lower attenuation, such as −1 dB or −20 dB. The attenuation is set by frequency scaling the denominator. This can be approximated by interpolating Bessel filter attenuation tables, or calculated precisely with Newton's method or the use of a root finding algorithm.

The denominator may be frequency scaled by using a scaling factor that will be referred to as ${\displaystyle \omega _{c}}$. The frequency-scaled denominator of the Bessel transfer function may be rewritten for the ${\displaystyle n=3}$ case above and the example below as follows:

${\displaystyle H(s)'=H(s)_{{\text{desired dB at }}\omega =1}={\frac {15}{(s\omega _{c})^{3}+6(s\omega _{c})^{2}+15s\omega _{c}+15}}}$

The task for modifying ${\displaystyle H(s)}$ is to find an ${\displaystyle \omega _{c}}$ that results in the desired attenuation at 1 r/s for the normalized transfer function. Newton's method can be summarized to do this by its basic definition:

${\displaystyle \omega _{a}=\omega _{a}+(|H(j\omega _{a})|-{\text{DesiredAttenuation}})/(d|H(j\omega _{a})|/dj\omega _{a})}$

through successive iterations of ${\displaystyle \omega _{a}}$ until ${\displaystyle \omega _{a}}$is the frequency that attenuates ${\displaystyle |H(j\omega _{a})|}$ to the desired attenuation.

While the above summary may be concise and easily understood, the mechanics of obtaining an accurate derivative of the magnitude function along the ${\displaystyle j\omega }$ axis may be problematic. Digital techniques may be used, but it is generally better to apply analytical techniques with continuous functions, if possible, so as to maximize accuracy. ${\displaystyle \omega _{a}}$ may be determined by either of the following methods.

Finding attenuation frequency with Newton's method

Newton's method requires a known magnitude value and derivative magnitude value for the for ${\displaystyle |H(j\omega a)|}$. To obtain these values, make the following alterations:

1. Multiply ${\displaystyle H(s)}$ by ${\displaystyle H(-s)}$ to obtain ${\displaystyle H(s)H(-s)}$. This will eliminate complex numerical results when evaluating ${\displaystyle |H(s)H(-s)|}$ by removing the odd order terms in the polynomial.
2. In ${\displaystyle H(s)H(-s)}$, negate all terms of ${\displaystyle s^{n}}$ when ${\displaystyle (n+2)}$ is divisible by ${\displaystyle 4}$. That would be ${\displaystyle s^{2}}$, ${\displaystyle s^{6}}$, ${\displaystyle s^{10}}$, and so on. We will call the modified function ${\displaystyle H_{2}(s)H_{2}(-s)}$, and this modification will allow the use of real numbers instead of complex numbers when evaluating the polynomial and its derivative. That is, we can now use the real ${\displaystyle \omega _{a}}$in place of the complex ${\displaystyle j\omega _{a}}$
3. Convert the desired attenuation in dB, ${\displaystyle A_{dB}}$, to an arithmetic attenuation value, ${\displaystyle A_{arith}}$, using ${\displaystyle A_{arith}=10^{A_{dB}/20}}$. For example, −3.010 dB is 0.7071, −1 dB is 0.8913 and so on. This simplifies the derivative evaluation.
4. Square the resulting arithmetic attenuation to correlate with the squaring of the ${\displaystyle H(s)}$ function.
5. Calculate the modified ${\displaystyle [H_{2}(\omega _{a})H_{2}(-\omega _{a})]}$ in Newton's method using the real value, ${\displaystyle \omega _{a}}$. Always take the absolute value.
6. Calculate the derivative the modified ${\displaystyle H_{2}(\omega _{a})H_{2}(-\omega _{a})}$ with respect to the real value, ${\displaystyle \omega _{a}}$. Negate the derivative computation to account for the effects due to the modifications made to create ${\displaystyle H_{2}(\omega _{a})H_{2}(-\omega _{a})}$. DO NOT take the absolute value of the derivative.

When steps 1) through 4) are complete, the expression involving Newton's method may be written:

${\displaystyle \omega _{a}=\omega _{a}+([H_{2}(\omega _{a})H_{2}(-\omega _{a})]-{\text{DesiredAttenuation}})/(d[H_{2}(\omega _{a})H_{2}(-\omega _{a})]/d\omega _{a})}$

using a real value for ${\displaystyle \omega _{a}}$with no complex arithmetic needed. The movement of ${\displaystyle \omega _{a}}$ should be limited to prevent it from going negative early in the iterations for increased reliability. When complete, ${\displaystyle \omega _{a}}$ can used for the ${\displaystyle \omega _{c}}$ that can be used to scale the original ${\displaystyle H(s)}$ transfer function denominator. The attenuation of the modified ${\displaystyle H(s)'}$ will then be virtually the exact desired value at 1 rad/sec. If performed properly, only a handful of iterations are needed to set the attenuation through a wide range of desired attenuation values for both small and very large order Bessel filters.

Finding attenuation frequency from the roots

Since ${\displaystyle |H(j\omega _{a})|}$ does not contain any phase information, directly factoring the transfer function will not produce usable results. However, the transfer function may be modified by multiplying it with ${\displaystyle H(-s)}$ to eliminate all odd powers of ${\displaystyle H(j\omega a)}$, which in turn forces ${\displaystyle H(j\omega a)}$ to be real at all frequencies, and then finding the frequency that result on the square of the desired attention.

1. Multiply ${\displaystyle H(s)}$ by ${\displaystyle H(-s)}$ to obtain ${\displaystyle H(s)H(-s)}$.
2. Convert the desired attenuation in dB, ${\displaystyle A_{dB}}$, to a squared arithmetic attenuation value, ${\displaystyle A_{arith}^{2}}$, by using ${\displaystyle A_{arith}^{2}=10^{A_{dB}/10}}$.
3. Find ${\displaystyle P(S)=H_{num}(S)H_{num}(-S)-A_{arith}^{2}H_{den}(S)H_{den}(-S)}$
4. Find the roots of P(S) using a root finding algorithm.
5. Of the set of roots from above, there will be one and only one positive imaginary root for a function such as Bessel filters that continuously attenuates as frequency increases. The imaginary value of this root will be frequency that results in the desired attenuation, ${\displaystyle \omega c}$.

For example, setting the example 3rd order Bessel filter to -20dB at a cutoff frequency of 1 rad/sec is performed as follows:

${\displaystyle {\begin{aligned}&A_{arith}^{2}=10^{-(-20)/10}=0.01\\&H(s)={\frac {15}{s^{3}+6s^{2}+15s+15}}\\&H(s)H(-s)={\frac {225}{-s^{6}+6s^{4}-45^{2}s+225}}\\&P(s)=225-0.01(-s^{6}+6s^{4}-45^{2}s+225)=0.01S^{6}--0.06S^{4}+0.45S+222.75\\&R_{positive,imaginary}=j5.0771344\\&\omega _{-20dB{\text{ atten}}}=\omega _{c}=5.0771344{\text{ rad/sec}}\\&H(s)'=H(s)_{-20dB{\text{ at }}\omega =1}={\frac {15}{(5.0771344^{3})s^{3}+(6*5.0771344^{2})s^{2}+(15*5.0771344)s+15}}\\&\\&{\text{Check:}}\\&|H(j)'|={\bigg |}{\frac {15}{(5.0771344^{3})j^{3}+(6*5.0771344^{2})j^{2}+(15*5.0771344)j+15}}{\bigg |}=0.1=-20dB\end{aligned}}}$

Example

Gain plot of the third-order Bessel low-pass filter, versus normalized frequency.

Group delay plot of the third-order Bessel low-pass filter, illustrating flat unit delay in the passband.

The transfer function for a third-order (three-pole) Bessel low-pass filter with ${\displaystyle \omega _{0}=1}$ is

${\displaystyle H(s)={\frac {15}{s^{3}+6s^{2}+15s+15}},}$

where the numerator has been chosen to give unity gain at zero frequency (${\displaystyle s=0}$).The roots of the denominator polynomial, the filter's poles, include a real pole at ${\displaystyle s=-2.3222}$, and a complex-conjugate pair of poles at ${\displaystyle s=-1.8389\pm j1.7544}$, plotted above.

The gain is then

${\displaystyle G(\omega )=|H(j\omega )|={\frac {15}{\sqrt {\omega ^{6}+6\omega ^{4}+45\omega ^{2}+225}}}.\,}$

The −3-dB point, where ${\displaystyle |H(j\omega )|={\frac {1}{\sqrt {2}}},\,}$ occurs at ${\displaystyle \omega =1.756}$. This is conventionally called the cut-off frequency.

The phase is

${\displaystyle \phi (\omega )=-\arg(H(j\omega ))=\arctan \left({\frac {15\omega -\omega ^{3}}{15-6\omega ^{2}}}\right).\,}$

The group delay is

${\displaystyle D(\omega )=-{\frac {d\phi }{d\omega }}={\frac {6\omega ^{4}+45\omega ^{2}+225}{\omega ^{6}+6\omega ^{4}+45\omega ^{2}+225}}.\,}$

The Taylor series expansion of the group delay is

${\displaystyle D(\omega )=1-{\frac {\omega ^{6}}{225}}+{\frac {\omega ^{8}}{1125}}+\cdots .}$

Note that the two terms in ${\displaystyle \omega ^{2}}$ and ${\displaystyle \omega ^{4}}$ are zero, resulting in a very flat group delay at ${\displaystyle \omega =0}$. This is the greatest number of terms that can be set to zero, since there are a total of four coefficients in the third-order Bessel polynomial, requiring four equations in order to be defined. One equation specifies that the gain be unity at ${\displaystyle \omega =0}$ and a second specifies that the gain be zero at ${\displaystyle \omega =\infty }$, leaving two equations to specify two terms in the series expansion to be zero. This is a general property of the group delay for a Bessel filter of order ${\displaystyle n}$: the first ${\displaystyle n-1}$ terms in the series expansion of the group delay will be zero, thus maximizing the flatness of the group delay at ${\displaystyle \omega =0}$.

Digital

Although the bilinear transform is used to convert continuous-time (analog) filters to discrete-time (digital) infinite impulse response (IIR) filters with comparable frequency response, IIR filters obtained by the bilinear transformation do not have constant group delay. ^[10] Since the important characteristic of a Bessel filter is its maximally-flat group delay, the bilinear transform is inappropriate for converting an analog Bessel filter into a digital form.

The digital equivalent is the Thiran filter, also an all-pole low-pass filter with maximally-flat group delay, ^{[11] [12]} which can also be transformed into an allpass filter, to implement fractional delays. ^{[13] [14]}

See also

• Bessel function
• Butterworth filter
• Chebyshev filter
• Comb filter
• Elliptic filter
• Group delay and phase delay

References

1. "Bessel Filter". 2013. Archived from the original on 2013-01-24. Retrieved 2022-05-14.
2. Thomson, W. E. (November 1949). "Delay networks having maximally flat frequency characteristics" (PDF). Proceedings of the IEE - Part III: Radio and Communication Engineering. 96 (44): 487–490. doi:10.1049/pi-3.1949.0101.
3. Roberts, Stephen (2001). "Transient Response and Transforms: 3.1 Bessel-Thomson filters" (PDF).
4. "comp.dsp | IIR Gaussian Transition filters". www.dsprelated.com. Retrieved 2022-05-14. An analog Bessel filter is an approximation to a Gaussian filter, and the approximation improves as the filter order increases.
5. "Gaussian Filters". www.nuhertz.com. Archived from the original on 2020-01-11. Retrieved 2022-05-14.
6. "How to choose a filter? (Butterworth, Chebyshev, Inverse Chebyshev, Bessel–Thomson)". www.etc.tuiasi.ro. Retrieved 2022-05-14.
7. "Free Analog Filter Program". www.kecktaylor.com. Retrieved 2022-05-14. the Bessel filter has a small overshoot and the Gaussian filter has no overshoot.
8. Paarmann, L. D. (2001). Design and Analysis of Analog Filters: A Signal Processing Perspective. Springer Science & Business Media. ISBN   9780792373735. the Bessel filter has slightly better Shaping Factor, flatter phase delay, and flatter group delay than that of a Gaussian filter of equal order. However, the Gaussian filter has less time delay, as noted by the unit impulse response peaks occurring sooner than they do for Bessel filters of equal order.
9. Bianchi, Giovanni; Sorrentino, Roberto (2007). Electronic filter simulation & design. McGraw–Hill Professional. pp. 31–43. ISBN   978-0-07-149467-0.
10. Zhang, Xi (2008-07-01). "Design of maximally flat IIR filters with flat group delay responses". Signal Processing. 88 (7): 1792–1800. doi:10.1016/j.sigpro.2008.01.016. ISSN   0165-1684.
11. Thiran, J.-P. (1971). "Recursive digital filters with maximally flat group delay". IEEE Transactions on Circuit Theory. 18 (6): 659–664. doi:10.1109/TCT.1971.1083363. ISSN   0018-9324.
12. Madisetti, Vijay (1997). "Section 11.3.2.2 Classical IIR Filter Types". The Digital Signal Processing Handbook. CRC Press. p. 11-32. ISBN   9780849385728.
13. Smith III, Julius O. (2015-05-22). "Thiran Allpass Interpolators". W3K Publishing. Retrieved 2022-05-14.
14. Välimäki, Vesa (1995). Discrete-time modeling of acoustic tubes using fractional delay filters (PDF) (Thesis). Helsinki University of Technology.

External links

• Bessel and Linear Phase Filters — Nuhertz
• Bessel Filter Constants — C.R. Bond
• Bessel Filters Polynomials, Poles and Circuit Elements — C.R. Bond
• Java source code to compute Bessel filter poles