Bessel filter analog to digital

I need to verify the transfer function of a 4th order bessel filter and modify it if needed. I can create the analog filter, convert to digital and look at the trasffer function iof the analog filter, but can't figure out how tio compare to the digital converted filter. the analòog filter TF plot is ok. With the digital filter I am not able to create a plot as for the analog filter, x axis seems to be normalyzed to 1 Hz end frequency, then the 2 TFs can't be compared too well. any suggestion is wellcome

thanks

fs = 330e9;  % sampling rate to have at least 10 samples per period
fc = 33e9;   % cutt off frequency
% Bessel filter is analog filter
w0 = 2*pi*fc;  % cut-off as rad/second
[zb,pb,kb] = besself(4, w0,'low'); % zero, pole and gain form bessel filter
[zd,pd,kd] = bilinear(zb,pb,kb,fs);  % Convert to digital filterz-domain zero/pole/gain
% Convert zero-pole-gain filter parameters to transfer function form for
% both A and D filters
[bb,ab] = zp2tf(zb,pb,kb); 
[bd,ad] = zp2tf(zd,pd,kd);
[hb,wb] = freqs(bb,ab,5000); % Frequency response and phase of analog filter
[hd,wd]  = freqz(bd,ad,5000;  % Frequency respèonse and phase of analog filter

[sosd,gd] = zp2sos(zd,pd,kd); % coefficients of digital filter to be used in T&M product SW filter
% first TF plot
mag1 = abs(hb);
phase1 = angle(hb);
phasedeg1 = phase1*180/pi;
subplot(2,1,1)
plot(wb,mag1)
grid on
xlabel('Frequency (rad/s)')
ylabel('Magnitude')
title ("Bessel 4th order analog filter")
subplot(2,1,2)
plot(wb,phasedeg1)
grid on
xlabel('Frequency (rad/s)')
ylabel('Phase (degrees)')

% second  TF plot
figure
mag2 = abs(hd);
phase2 = angle(hd);
phasedeg2 = phase2*180/pi;
subplot(2,1,1)
plot(wd,mag2)
grid on
xlabel('Frequency (rad/s)')
ylabel('Magnitude')
subplot(2,1,2)
plot(wd,phasedeg2)
grid on
xlabel('Frequency (rad/s)')
ylabel('Phase (degrees)')
title ("Bessel 4th order Digital filter")

Answer 1

[Paul]

Hi Maurizio,

freqs returns the frequency vector, wb, in units of rad/sec.

freqz, as used below, returns the frequency vector, wd, in units of rad/sample.

For comparing the results wd needs to be divided by the sampling period, or multiplied by the sampling frequency, to convert to rad/sec.

fs = 330e9;  % sampling rate to have at least 10 samples per period
fc = 33e9;   % cutt off frequency
% Bessel filter is analog filter
w0 = 2*pi*fc;  % cut-off as rad/second
[zb,pb,kb] = besself(4, w0,'low'); % zero, pole and gain form bessel filter
[zd,pd,kd] = bilinear(zb,pb,kb,fs);  % Convert to digital filterz-domain zero/pole/gain
% Convert zero-pole-gain filter parameters to transfer function form for
% both A and D filters
[bb,ab] = zp2tf(zb,pb,kb); 
[bd,ad] = zp2tf(zd,pd,kd);
[hb,wb] = freqs(bb,ab,5000); % Frequency response and phase of analog filter
[hd,wd]  = freqz(bd,ad,5000);  % Frequency respèonse and phase of analog filter
[sosd,gd] = zp2sos(zd,pd,kd); % coefficients of digital filter to be used in T&M product SW filter
% first TF plot
mag1 = abs(hb);
phase1 = angle(hb);
phasedeg1 = phase1*180/pi;
subplot(2,1,1)

Use semilogx instead of plot to recreate plot in the question

%plot(wb,mag1)
semilogx(wb,mag1)
grid on
xlabel('Frequency (rad/s)')
ylabel('Magnitude')
title ("Bessel 4th order analog filter")
subplot(2,1,2)
%plot(wb,phasedeg1)
semilogx(wb,phasedeg1)
grid on
xlabel('Frequency (rad/s)')
ylabel('Phase (degrees)')

% second  TF plot
figure
mag2 = abs(hd);
phase2 = angle(hd);
phasedeg2 = phase2*180/pi;
subplot(2,1,1)
%plot(wd,mag2)

Multiply wd by fs to convert to rad/sec (rad/sample * sample/sec). Set the xlimits to match the freqs results.

semilogx(wd*fs,mag2)
grid on
xlim([1e10 1e12])
xlabel('Frequency (rad/s)')
ylabel('Magnitude')
subplot(2,1,2)
%plot(wd,phasedeg2)
semilogx(wd*fs,phasedeg2)
grid on
xlim([1e10 1e12])
xlabel('Frequency (rad/s)')
ylabel('Phase (degrees)')
title ("Bessel 4th order Digital filter")

[Maurizio]

this helped..and looks much better if I set the horizontal axis in Hz instead of in rad/s. The 2 filters are not identical but pretty close.

thanks for the correctlion

[Paul]

Discretization of any analog filter (except a constant gain I suppose) will never yield the same frequency response as the analog filter. The question is whether or not the frequency response of the discretized filter is suitable for your application.

From a visual perspective, it would be better to overlay the responses on the same axes for comparison.

Good luck with your project.

[Maurizio]

the 2 filters are pretty similar if the sample rate is much higher than the cutoff frequency point....Overlapping the e plots, helps, but something does not look right. The -3dB ( 0.7 in amplitude) point is supposed to be at fc. I set the fc = 20GHz for semplicity, and the 0.7 is at 13.1GHz, not 20GHz, what am I doing wrong? at 20GHz the gain is 0.4 !!

[Maurizio]

Looks like I missunderstood the shape of a Bessel filter, and the purpose of this filter is not the gain, but phase liearity. The attenuation at the fc is not -3dB as for other filters. Only the first order bessel filter has-3dB attenuation at fc, higher order filters starts rolling off at lower frequency.

[Paul]

Apparently the definition of the cutoff frequency is non-standard (Reference).

Because you have already determined the 0.7 gain point is at 13.1 GHz, you can shift it to 20 GHz by scaling the coefficients of the denominator. Could also figure out the 13.1 GHz programatically.

fs = 330e9;  % sampling rate to have at least 10 samples per period
fc = 20e9;   % cut off frequency
% Bessel filter is analog filter
w0 = 2*pi*fc;  % cut-off as rad/second
[b,a] = besself(4, w0,'low');

Define scale factor

r = 13.1e9/20e9;

Scale the denominator

format short e
a1 = a.*r.^(4:-1:0)

a1 = 1×5

1.0e+00 * 1.8406e-01 1.1032e+11 2.9752e+22 4.1607e+33 2.4937e+44

b1 = b;

Rescale so that a1 has unity leading coefficient (not necessary)

b1 = b1/a1(1);
a1 = a1/a1(1);

Or you could twiddle the w0 input to besslef.

[b3,a3] = besself(4, w0/r,'low');


[h,wout] = freqs(b,a);
h1       = freqs(b1,a1,wout);
h3       = freqs(b3,a3,wout);

fout = wout/2/pi;

figure
subplot(211);
plot(fout,abs(h),fout,abs(h1),fout,abs(h3)),grid;
yline(0.7);xline(fc)
subplot(212);
plot(fout,unwrap(angle(h)),fout,unwrap(angle(h1)),fout,unwrap(angle(h3))),grid

BTW, it's probably best to keep everything in zpk form.

[Bruno Luong]

Note the the bilinear transform has non-linear - but one-to-one - relationship between analog and digital frequency

Omega_analog = 2/T tan(omega_digital/2).

So for that reason alone the frequency response cannot be matched at high frequency side. But the more T (digital sample period) is reduced the more linear part is extended father.

[Paul]

@Bruno Luong

Every discretization method that I'm aware of will distort the frequency reponse of the discretized approximation relative to that of the analog design. I gather from your response below that you prefer the impulse invariance method. Here is the comparison of all the techniques supported by the Control System Toolbox, except for Tustin with prewarp (I used a lower cutoff for readability).

fs = 1e3;  % sampling rate to have at least 10 samples per period
fc = 1e2;   % cut off frequency
% Bessel filter is analog filter
w0 = 2*pi*fc;  % cut-off as rad/second
[zc,pc,kc] = besself(4, w0,'low'); % zero, pole and gain form bessel filter
hc = zpk(zc,pc,kc);
method = {'zoh','foh','impulse','tustin','matched','least-squares'};
wc = logspace(0,4,1000);
[m,p] = bode(hc,wc);
f1 = figure;
ax1 = subplot(211);
semilogx(ax1,wc/2/pi,20*log10(m(:)));hold on
ax2 = subplot(212);
semilogx(ax2,wc/2/pi,180/pi*angle(exp(1j*p(:)*pi/180)));hold on

wd = wc(wc<pi*fs);

for ii = 1:numel(method)
    [m,p] = bode(c2d(hc,1/fs,method{ii}),wd);
    subplot(211);semilogx(ax1,wd/2/pi,20*log10(m(:)))
    subplot(212);semilogx(ax2,wd/2/pi,180/pi*angle(exp(1j*p(:)*pi/180)))
end

subplot(ax1),grid,xline(fs/2);ylim([-100 0]),ylabel('dB'),xline(fc)
legend(horzcat({'cont'},method),"Location",'SouthWest')
subplot(ax2),grid,xline(fs/2);ylabel('Degrees'),xline(fc)
xlabel('Frequency (Hz)')
copyobj(get(f1,'children'),figure)

h = get(gcf,'Children')

h = 
  3×1 graphics array:

  Axes
  Legend    (cont, zoh, foh, impulse, tustin, matched, least-squares)
  Axes

set(h([1 3]),'XLim',[1e2 1e3])
set(h(2),'Location','SouthWest')

Impulse invariance certainly looks pretty good for both gain and phase close to the Nyquist frequency.

One of the features of the Bessel filter is its nearly linear phase repsonse in the passband. The negative of the slope of the phase response is the group delay, which should be constant over frequencies where the phase is linear. Following is a comparison of the group delay for the analog filter and its discrete approximations.

fvec = linspace(0,500,1000);
[m,p] = bode(hc,fvec*2*pi);
gd = -gradient(p(:)*pi/180)./gradient(2*pi*fvec(:))*1000;  % msec
f3 = figure;
ax3 = gca;
plot(ax3,fvec,gd,'-o','MarkerIndices',1:50:numel(gd));hold on

for ii = 1:numel(method)
    hd = tf(c2d(hc,1/fs,method{ii}));
    b = hd.Num{:};
    a = hd.Den{:};
    gd = grpdelay(b,a,fvec,fs)/fs*1000;  % group delay in ms
    plot(ax3,fvec,gd)
end

grid,xline(fc);xlim([0 150])
legend(horzcat({'cont'},method),"Location",'SouthWest')
xlabel('Frequency (Hz)')
ylabel('Group Delay (msec)')

We see that the tustin transformation is the worst at maintaining constant group delay. The matched and zoh methods are pretty good, but (interestingly) introduce a group delay bias. The foh, impulse, and least squares methods seem to yiled a group delay basically identical to that of the analog filter.

[Bruno Luong]

@Paul Thanks. I don't have access to this toolbox so the result is precious for me.

Yes my own analysis (long ago) on Bessel filter makes me take the decision to use impulse invarirant. I'm glad your result confirms my finding.

EDIT: For my purpose, it is not critical that the group delay of various discretizations do not match that of the continuous one, the bias of zoh and matched GRD are OK.

I only care whereas the group delay stays constant before the cut-off. The bias in GRD implies the bias of the cut-off as well. So we could see zoh and matches as methods that do not provide very accurate cut-off of the original analog Bessel.

Answer 2

[Star Strider]

Simply stated, there is no reliable method of digitsing a Bessel filter, and even if you managed to get some sort of result using the Tustin approximation (the preferred method), a digital version loses the Bessel filter linear phase characteristics. (This is discussed in every digital signal processing text that I have ever read.) Bessel filters can only be realised in hardware.

If you want a linear phase response from any digital filter created in MATLAB, use the filtfilt function to do the actual filtering.

[Bruno Luong]

@Star Strider Can you please provide reference or argument for Tustin as prefered method of digital conversion of Bessel filter?

In one of my applications I use impulse invariant method and I would like to revisit and see whereas I select the wrong method.

Thanks

[Bruno Luong]

Based on the script made by @Paul here is the Bode's plot for comparison the Tustin and Impulse Invariant on Bessel filter

fs = 330e9;  % sampling rate to have at least 10 samples per period
fc = 33e9;   % cutt off frequency
% Bessel filter is analog filter
w0 = 2*pi*fc;  % cut-off as rad/second

[zb,pb,kb] = besself(4, w0,'low'); % zero, pole and gain form bessel filter

% Tustin digital converter
[zd,pd,kd] = bilinear(zb,pb,kb,fs);  % Convert to digital filterz-domain zero/pole/gain
% both A and D filters
[bb,ab] = zp2tf(zb,pb,kb); 
[bd,ad] = zp2tf(zd,pd,kd);

[bc,ac] = impinvar(bb,ab,fs);  % Convert to digital filterz-domain using impulse invariant method

[hb,wb] = freqs(bb,ab,5000); % Frequency response and phase of analog filter
[hd,wd]  = freqz(bd,ad,5000);  % Frequency respèonse and phase of bilinear digital filter
[hc,wc]  = freqz(bc,ac,5000);  % Frequency respèonse and phase of impinvar digital filter

% first TF plot
mag1 = abs(hb);
phase1 = unwrap(angle(hb));
phasedeg1 = rad2deg(phase1);

mag2 = abs(hd);
phase2 = unwrap(angle(hd));
phasedeg2 = rad2deg(phase2);

mag3 = abs(hc);
phase3 = unwrap(angle(hc));
phasedeg3 = rad2deg(phase3);

close all
subplot(2,1,1)
semilogx(wb,mag1,'LineWidth', 3)
hold on
semilogx(wd*fs,mag2,'LineWidth', 2)
semilogx(wc*fs,mag3, 'g', 'LineWidth', 1)
xlim([1e10 1e12])
grid on
xlabel('Frequency (rad/s)')
ylabel('Magnitude')
title ("Bessel 4th order filter amplitude")
legend('Analog', 'Digital Tustin ', 'Digital Impulse invariant', 'Location', 'best')

subplot(2,1,2)
semilogx(wb,phasedeg1,'LineWidth', 3)
hold on
semilogx(wd*fs,phasedeg2,'LineWidth', 2)
semilogx(wc*fs,phasedeg3, 'g', 'LineWidth', 1)
grid on
xlim([1e10 1e12])
xlabel('Frequency (rad/s)')
ylabel('Phase (degrees)')
title ("Bessel 4th order filter phase")
legend('Analog', 'Digital Tustin ', 'Digital Impulse invariant', 'Location', 'best')

It looks like impinvar is slighly better than bilinear (Tustin), at least for thiis specific case.

[Star Strider]

I generally prefer Tustin, however it likely does not make much difference. The point is that every digital signal processing text I have ever consulted (and the MATLAB documentation: ‘The besself function does not support the design of digital Bessel filters.’) discourages digitising Bessel filters simply because they do not retain their linear phase properties when digitised.

[Bruno Luong]

Yes, but the Bode plot shows the digital version does a quite decend job with fs >= 10*fcutoff. It is just a question of able to sample fast enough.

I don't understand why it is discourage this particular filter design.

Nowaday digital filter is used more than anlogue filter for obvious reason that many signal get digitalized as soon as the analog signal is acquired. The focus should be shifted to digital processing.

[Star Strider]

‘I don't understand why it is discourage this particular filter design.’

The purpose of using a Bessel filter is to take advantage of its phase-neutral characteristics, the reason that as an analog filter it is almost universally used as an anti-aliasing filter. In that application, it eliminates frequencies higher than the Nyquist frequency (of the following ADC) while not altering the phase relationships of the frequencies in the passband, so there is no phase distortion as there would be in other filter designs. Digital implementations of the Bessel design lose this characteristic. Also, as a digital filter, it is less computationally efficient than a comparable elliptic filter (likely the best design for a digital filter).

In MATLAB, all filters can be rendered essentially phase-neutral with the filtfilt function, so digitising a Bessel filter to take advantage of the linear-phase characteristic that it loses when digitised simply makes no sense.

.

[Bruno Luong]

filtfilt is not causal filter, in practice it is not desirable when data come in stream. I hope the self driving cars don't way the acquisition buffer to be filled before taking the decision.

I'm not argue bessel filter should not be implemented in analogue, I argue it can be applied in digital as well especially the data is never exacty sampling,it's already integrated over the sensor window. It's just one way to design a digital filter, close to Bessel analogue filter. Who cares if it's doesn't have linear phase above 100 fs, where the signal perhaps contains mostly HF noise?

MATLAB is the tool. They should not discourage people for using the tool. It's a job of the designer to decide how he uses the tool. The professor job is teaching how evaluate and quantify how good or bad a digital filter is compared to a analgue filter, not blindly telling the students never using digitalized bessel filter.

[Star Strider]

The advise against digitising Bessel filters is in every digital signal processing book I have ever read, for the reasons I have stated. The filtfilt function may not be perfect (it is a forward-backward filter approach that eliminates the phase distortion of functions such as filter), however it works for the purposes it is intended to address.

If you want to digitise a Bessel filter (or do anyting else, for that matter), you have my permission. Just don’t be surprised if the effort results in a complete failure.

EDIT — (5 Aug 2023 at 13:26)

The purpose of signal processing is to get a filtered signal ideally without distortion in the passband. Details such as the phase chracteristics are important, especially with physiological signals (that I usually work with).

[Bruno Luong]

"If you want to digitise a Bessel filter (or do anyting else, for that matter), you have my permission. Just don’t be surprised if the effort results in a complete failure. "

Sorry but I have also successully used the Bessel digital filter in industrial SW processing. I don't need your permisson.

So far you haven't give a solid logical argument againts digital bessel filter, and the simple Bode plots I gave above does show any big issue to me.

[Star Strider]

Of course you can digitise it. However digitising it will not preserve its analog linear phase characteristics.

Using any filter design with filtfilt will result in a linear phase characteristic, regardless of the design.

[Bruno Luong]

"it will not preserve its analog linear phase characteristics"

My goal: the phase response is close to linear in the frequency bandwidth I'm interested in. That is what matter for my app. An enough fast acquisition and Bessel digitalized does the job for me. I don't need a perfectly straight line phase over the entire seismic to gamma-ray frequency range, those ideal theoretical is fun limited case, but I don't need it in practice.

No question I use filtfilt, I need causal filter and I need by filter signal continuously streamed, and the output filter signal must be smooth and not appending two non-overlapped intervals. The goal of my filer design is different than yours.

[Star Strider]

Thanks for the clarification, however for my applications, ‘close’ is just not good enough since phase distortion can be a significant problem for my signals. I am generally interested in eliminating baselline variations and line (mains) frequency noise, so I either use bandpass filters with passbands of 1 to 50 Hz, or broader bandpass filters to eliminate baseline variations and high-frequency noise, with notch filters to eliminate the mains frequency noise.

My arguments (and my objections to digitising Bessel filters) nevertheless appliy generally.

[Bruno Luong]

I can understand that in physiology signal the HF content has an important information.

In my application case I'm only interested in detecting precisely when a specific signal occurs (and report it without minimum delay) and the bandwidth of those expected signal is relatively narrow. I work on the frequency about 1-10 MHz.

[Paul]

@Bruno Luong

I agree with you 100+% that the question is whether or not the discretized approximation to the Bessel filter is good enough for its intended application. Citing an accessible reference to such an application where the discretized Bessel filter is used would be a wonderful addition to this forum because, IMO, the idea that a Bessel filter might be useful for many different applications is rarely acknowledged. If you don't mind me asking, can you expand on your industrial software application? What is the software doing?

Also, I agree that filtfilt is not a viable solution in all applications, which is also rarely acknowledged in this forum.

My only addition here is that even though filtfilt is a useful tool, it does not perform zero-phase filtering despite the claim to the contrary on its doc page. I showed that in this Question filtfilt Does Not Perform Zero-Phase Filtering. Does It Matter? (in which you participated, but linking here for others who may be reading this thread).

[Star Strider]

Probably any filter design would work if you are only interested in detecting the presence of the signal, since its phase characteristics might not be important. An elliptic filter would be more computationally efficient, since elliptic filters are generally shorter.

[Bruno Luong]

I want detect very precisely when the signal occurs and the delay of my detection as well to correct when the causing event occurs.

@Paul: the application is to estimate in real time the vibration-induced mechanical stress of the blades in a turbine (when the engine is runing). We install 3-6 sensors around the engine casing and measure when the blades are passing. We can then know the amplitude and frequeny of blade vibrations and then derive the stress amplitude and natural modal frequency. This application applies for both avionic or energetic industry.

[Paul]

Very interesting. You use the discretized Bessel filter on the outputs of the sensors?

[Star Strider]

Bruno —

‘I want detect very precisely when the signal occurs and the delay of my detection as well to correct when the causing event occurs.’

Wouldn’t that be easier to do in hardware?

[Bruno Luong]

"Bessel filter on the outputs of the sensors"

Yes. The filter is only a small part of the processing but it must be done with care.

[Bruno Luong]

@Star Strider "Wouldn’t that be easier to do in hardware? "

No because the Bessel filter cutoff is adaptative to the engine speed and other thing. Digital filter (SW) is more flexible to cope with those dependencies.

Furthermore the noise after the HW filter, e.g. ADC stage, would not be removed.

[Star Strider]

That was not obvious from the initial description.

[Paul]

At the risk of stating the obvious, the phase of the analog, low pass, Bessel filter is not perfectly linear over the pass band. So is it not true that, even in the analog world, one has to make value judgements about the adequacy of the phase response of the filter for any specific application (and trade that off against complexity of the filter implementation)? Conceptually, I don't see how that's any different than making that same judgement about the phase response of the discretized filter.

[Star Strider]

It may not perfectly phase-linear throughout the entire passband, however it has significant advantages in that respect over other filter designs specifically because it significantly minimises phase dostortion in the passband. The discretised filter loses the linear phase characteristic, and other digital filter designs are more computationally efficient.

That is the last I am going to say about this.

[Bruno Luong]

That's true (EDIT: subjected to Paul's message, not Star Strider). The Bessel analogue filter is the "best" in the sense that the Taylor expansion of the group delay (GRD) at Omega = 0 has all 0 terms up to the filter order allows. But of course it can not achieve totally flat GRD.

One can decide the "best" (flatnest of GRD) criteria differently (there are many ways), instead of Taylor expansion.

I wonder if one could start to design from scratch a IIR digital filter at certain given order and solve for flatness GRD - Taylor sense - same as analogue filter?

PS: after written the above I just doing some google and it seems that some people start to do such optimization with Butterworth/Chebyshev/Elliptic filters, http://article.nadiapub.com/IJSIP/vol9_no7/21.pdf, none of the GRD they present look as nice as Bessel filters (analogue or discreted)

[Bruno Luong]

@Star Strider "The discretised filter loses the linear phase characteristic"

I don't see it - at least in any noticable degree - in the Bode's plot I made above to compare impinvar and bilinear and in Paul's even more comprehhensive methods of discretizations.

[Paul]

@Bruno Luong

I edited my comment above to clean it up a bit and add a plot comparing the group delay that results from the different discretization methods.

Unfortunately, I got "File Not Found (404)" error when I tried thta nadiapub.com link.