# Matlab: How to design digital equivalent for a lowpass Bessel filter (Thiran filter)?

In the process of applying a lowpass Bessel filter to my signal, I realized that besself function does not support the design of digital Bessel filters and the bilinear function can be used to convert an analog filter into a digital form, except for Bessel filters. The digital equivalent for Bessel filters is the Thiran filter. The only thing I know for my filter is that it should have less than 5GHz bandwidth ( let's say 3GHz bandwidth). I do appreciate if someone could help me write the code in Matlab. The code which I have so far is:

%lowpass filter
sig = MY SIGNAL;
sig_length = 5000001;                         % my signal length
fs  = 10000e9                                 % sampling rate
fc  = 3e9;                                    % cutt off frequency
order = 4;
wo = 2*pi*fc;
[z,p,k] = besself(order, wo,'low');           % zero, pole and gain form
% Convert to digital fileter
[zd,pd,kd] = bilinear(z,p,k,fs);          % z-domain zero/pole/gain
[sos,g] = zp2sos(zd,pd,kd);                   % convert to second order section
filteredSignal = filtfilt(sos, g, sig);


I don't have Matlab but the coefficients for the Thiran filter are given by the Gaussian hypergeometric function. If you have wxMaxima there already is a built-in function. If you'll run these two lines (the first one is only needed once), you'll get the denominator:

expand_hypergeometric : true$N : 4$
D : 1\$
hypergeometric([-N, 2*d], [2*D + N + 1], z);


N is the order and D the delay, in seconds. The above will show:

z^4/42-(4*z^3)/21+(9*z^2)/14-(8*z)/7+1


If Matlab doesn't have that function, you could implement it as:

$$H(z)=\dfrac{(2N)!}{N!}\dfrac{1}{\prod_{i=N+1}^{2N}{2D+i}}\dfrac{1}{\sum_{k=0}^N{(-1)^k\binom{N}{k}\left[\prod_{i=0}^N{\left(\dfrac{2D+i}{2D+k+i}z^{-k}\right)}\right]}} \tag{1}$$

Thenumerator will be the sum of the coefficients, for unity gain. It will be more numerically stable to split the denominator into 2nd order sections.

Also, as it is, the filter has very poor frequency attenuation, but if you add a numerator in $$z$$, you may get better results. See this for example, where the two zeroes can give a much better attenuation, at the cost of some extra delay. Or, if only a zero at Nyquist is fine, just add $$1+z^{-1}$$.

• Hi @Amy - Did this answer your question? If so please select it to close the question or comment on what remains to be answered. Thank you! Apr 16, 2022 at 12:58