Bessel filter analog to digital
30 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1448152/image.png)
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")
0 Commenti
Risposte (2)
Paul
il 2 Ago 2023
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")
8 Commenti
Paul
il 5 Ago 2023
Modificato: Paul
il 6 Ago 2023
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')
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
il 5 Ago 2023
Modificato: Bruno Luong
il 6 Ago 2023
@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.
Star Strider
il 2 Ago 2023
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.
25 Commenti
Bruno Luong
il 5 Ago 2023
Modificato: Bruno Luong
il 5 Ago 2023
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.
Vedere anche
Categorie
Scopri di più su Bessel functions in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!