procedure BesselAnalog(Order: Integer; z: TVec; p: TVec; out k: Double);
Design an analog Bessel (Thomson) lowpass prototype filter of order Order, optimised for maximally-flat group delay (linear phase) rather than a flat magnitude. Returns zero-pole-gain with all zeros at infinity (z empty); the denominator is the reverse Bessel polynomial theta_n(s): H(s)=theta_n(0)/theta_n(s), theta_n(s)=sum_(k=0)^n((2n-k)!)/(2^(n-k) k! (n-k)!) s^k, normalised so that the DC gain is 1. The phase response is approximately linear near DC; the magnitude rolls off more gently than Butterworth. Domain: 1 <= Order <= MaxIirOrder. Poles satisfy Re(p_k)<0 (stable). z and p must share precision or an exception is raised.
| # | Name | Type | Description |
|---|---|---|---|
| 1 | Order | Integer | |
| 2 | z | TVec | source TVec |
| 3 | p | TVec | source TVec |
| 4 | k | Double |
Result: stored in self (calling object)
Design analog Bessel prototype filter of order Order. Place the resulting transfer function in zero-pole form in Z (zeros), P (poles) and K (gain). The cutoff frequency of the prototype filter is preset to 1 rad/sec. The filter has all zeros in infinity. The transfer function is defined as([1], p. 230):
d0 H(s) = -------- Bn(s) (2*n)! n d0 = ------- , Bn(s) = Sum(d[k]s^k), k = 0,...,n 2^n*n! k=0 (2*n-k)! d[k] = -------------- , n = order of the filter 2^(n-k)(n-k)! Filter poles must be scaled with d0^(1/n)
Roots of the Bessel polynomial Bn(s) are found with the PolyRoots routine. Bessel lowpass filters are charachterized by the property that the group delay is maximally flat at the origing of the s-plane. ([1], p. 228).
References:
[1] Theory and application of digital signal processing, Lawrence R. Rabiner and Bernard Gold. Prentice-Hall, 1975.
uses MtxExpr, Math387, MtxVec, SignalUtils, MtxVecTee, MtxVecEdit, IirFilters,
LinearSystems;
procedure TForm1.Button1Click(Sender: TObject);
var z,p, num,den, FreqFr,Response: Vector;
Order: integer;
k,Wc,BW: double;
begin
Order := 5; //design a fifth order filter.
BesselAnalog(Order,z,p,k); //design analog protype
Wc := 1.1;
LowpassToLowpass(z,p,k,Wc); //frequency transformation in s-domain
ZeroPoleToTransferFun(num,den,z,p,k);
FreqFr.Length := 1000; //Define the frequency grid (logarithmic)
LogRamp(FreqFr,-1,1); //between 0.1 (=10^(-1)) and 10 (=10^1) rad/sec
FrequencyResponseS(num,den,FreqFr,Response); //Laplace
DrawIt(Response); //Y axis linear, X axis logarithmic;
end;