Bessel filter analog to digital

Question

Maurizio 2023-8-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2003602-bessel-filter-analog-to-digital

编辑： Bruno Luong 2023-8-6

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")

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2023-8-2

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2003602-bessel-filter-analog-to-digital#answer_1282317

在 MATLAB Online 中打开

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 个评论
显示 6更早的评论隐藏 6更早的评论

Paul 2023-8-4

编辑：Paul 2023-8-5

在 MATLAB Online 中打开

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.

Paul 2023-8-5

编辑：Paul 2023-8-6

在 MATLAB Online 中打开

@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 2023-8-5

编辑：Bruno Luong 2023-8-6

@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 2023-8-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2003602-bessel-filter-analog-to-digital#answer_1282307

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 个评论
显示 23更早的评论隐藏 23更早的评论

Bruno Luong 2023-8-5

编辑：Bruno Luong 2023-8-5

在 MATLAB Online 中打开

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.

Bruno Luong 2023-8-5

编辑：Bruno Luong 2023-8-5

@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 2023-8-6

编辑：Paul 2023-8-6

@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.

请先登录，再进行评论。

Bessel filter analog to digital

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

8 个评论
显示 6更早的评论隐藏 6更早的评论

25 个评论
显示 23更早的评论隐藏 23更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Bessel filter analog to digital

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

8 个评论 显示 6更早的评论隐藏 6更早的评论

25 个评论 显示 23更早的评论隐藏 23更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论

25 个评论
显示 23更早的评论隐藏 23更早的评论