The transfer function H(p) of a n 'th order Bessel filter is given by
Bn(0) H(p) = ------- Bn(p)
with the denominator polynomial
n n (2n - k)! p^k Bn(p) = sum c_k*p^k = sum ----------- * ------- (1) k=0 k=0 (n - k)!k! 2^(n-k)
and the numerator
(2n)! 1 Bn(0) = c_0 = ------- * ---- . (2) n! 2^n
Although the coefficients c_k are integer numbers, it is not advisable to use the polynomials in an unfactorized form because the coefficients are fast growing with order n (c_0 is approximately 0.3e24 and 0.8e59 for order n=20 and order n=40 respectively).
Therefore, the polynomial Bn(p) is factorized to first and second order polynomials with real coefficients corresponding to zeros and poles representation that is used in this library.
The function returns the coefficients which resulted from factorization of the normalized transfer function
H'(p') = H(p), p' = p/w0
as well as
alpha = 1/w0
the reciprocal of the cut of frequency w0 where the gain of the transfer function is decreased 3dB.
Both, coefficients and cut off frequency were calculated symbolically and were eventually evaluated with high precision calculation. The results were stored in this function as real numbers.
Equation
abs(H(j*w0)) = abs(Bn(0)/Bn(j*w0)) = 10^(-3/20)
which must be fulfilled for cut off frequency w = w0 leads to
[Re(Bn(j*w0))]^2 + [Im(Bn(j*w0))]^2 - (Bn(0)^2)*10^(3/10) = 0
which has exactly one real solution w0 for each order n. This solutions of w0 are calculated symbolically first and evaluated by using high precise values of the coefficients c_k calculated by following (1) and (2).
With w0, the coefficients of the factorized polynomial can be computed by calculating the zeros of the denominator polynomial
n Bn(p) = sum w0^k*c_k*(p/w0)^k k=0
of the normalized transfer function H'(p'). There exist n/2 of conjugate complex pairs of zeros (beta +-j*gamma) if n is even and one additional real zero (alpha) if n is odd. Finally, the coefficients a, b1_k, b2_k of the polynomials
a*p + 1, n is odd
and
b2_k*p^2 + b1_k*p + 1, k = 1,... div(n,2)
results from
a = -1/alpha
and
b2_k = 1/(beta_k^2 + gamma_k^2) b1_k = -2*beta_k/(beta_k^2 + gamma_k^2)
function BesselBaseCoefficients extends Modelica.Icons.Function; import Modelica.Utilities.Streams; input Integer order "Order of filter in the range 1..41"; output Real c1[mod(order, 2)] "[p] coefficients of Bessel denominator polynomials (a*p + 1)"; output Real c2[integer(order / 2), 2] "[p^2, p] coefficients of Bessel denominator polynomials (b2*p^2 + b1*p + 1)"; output Real alpha "Normalization factor"; end BesselBaseCoefficients;