How can I normalize my filter?

46 views (last 30 days)
Amy Lg
Amy Lg on 30 Mar 2022
Commented: Star Strider on 7 Apr 2022
Hi,
I am applying a low pass filter to my signal and I should make sure the transfer function of my filter is one. How can I check it?
If the filter is not normalized, how can I normalize it?
I am learning about filters and I appreciate any help.
%lowpass filter
sig = MY SIGNAL;
sig_length = 5000001; % my signal length
fs = 1e13 % 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);
  5 Comments
Amy Lg
Amy Lg on 6 Apr 2022
Thank you so much for your help.
Star Strider
Star Strider on 7 Apr 2022
Right.
In spite of everything I wrote in my answer (now deleted) and the observation in the Wikipedia article (and every signal processing textbook I consulted on this), you are still going to try to use the bilinear function to discretise a Bessel filter!
That is not going to work, and no amount of wishful thinking on your (or anyone else’s) part is ever going to make it work!

Sign in to comment.

Accepted Answer

Paul
Paul on 3 Apr 2022
Why is fs defined and then multiplied by 1e9 in the call to bilinear? Which I guess might be o.k. but perhaps isn't what was intended?
Anyway, without that scaling the discrete filter output from bilinear() looks like it has just about the same response as the continuous filter. To show that, I modified the code to use
fs = 10000e9 % sampling rate
fs = 1.0000e+13
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); % just use fs
Now compare the analog and discrete filters (I'm used to using Control System Toolbox functions, I'm sure the same plots can be obtained with Signal Processing Toolbox functions).
hc = zpk(z,p,k);
hd = zpk(zd,pd,kd,1/fs);
bode(hc,hd);
Pretty good match until close to the Nyquist freqency pi*fs = pi*1e13.
As for the question about normalization, I'm not quite sure what "make sure the transfer function of my filter is one" means. Clearly, the tf can't be one at all frequencies. If just looking to ensure the dc gain is one, then we can check this with
format long e
freqresp(hc,0)
ans =
1
freqresp(hd,0)
ans =
9.999999999999665e-01
which may be close enough to one for all practical purposes.
  8 Comments
Paul
Paul on 5 Apr 2022
Edited: Paul on 5 Apr 2022
Filter design and implementation:
fs = 10000e9; % sampling rate
fc = 3e9; % cut 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); % just use fs
Compare the frequency responses of the analog and discrete filters (I'm used to using Control System Toolbox functions, I'm sure the same plots can be obtained with Signal Processing Toolbox functions).
sysc = zpk(z,p,k);
sysd = zpk(zd,pd,kd,1/fs);
bode(sysc,sysd);
xline(wo)
We see a net phase change of 180 degrees up to the cut-off frequency.
According to the doc page for besself(), an important feature of the Bessel filter is that the group delay is nearly constant up to the cut-off frequency. Group delay is the derivative of the phase wrt to frequency, so if the group delay is constant, we expect the filter to have linear phase.
First, use the grpdelay() function get the group delay.
[gd,f] = grpdelay(zp2sos(zd,pd,kd),1e6,fs);
Verify that linear phase appears to be the case for both continuous- and discrete-time filters up to the cut-off.
[~,pd] = bode(sysd,2*pi*f);
pd = pi/180*squeeze(pd);
[~,pc] = bode(sysc,2*pi*f);
pc = pi/180*squeeze(pc);
plot(f,pc,f,pd)
xlabel('Frequency (Hz)')
ylabel('Phase (rad)')
xlim([0 fc])
legend('Cont','Disc')
Both filters look like the phase has constant slope, but hard to say for sure.
Compare the computed group delay to simple approximations of the slope of the phase responses.
plot(f,-gradient(pc,2*pi*f),f,-gradient(pd,2*pi*f),f,gd/fs)
xlabel('Frequency (Hz)')
ylabel('GD (sec)')
xlim([0 fc])
legend('Cont','Disc','grpdelay')
It appears that the discrete-time filter using bilinear() may be a good approximation to the continuous-time filter, at least up to the cut-off frequency. I think this result obtains because the Nyquist frequency is three orders of magnitude larger than the cut-off frequency, and in the low frequency region below the cut-off the bilinear transformation does prettty well at preserving the phase. Of course, this result might very well be not true if (fs/2)/fc were much smaller.
is it wrong if I use bessel filter for my simulation?
I assume this question really means: is it wrong if I use the discrete-time approximation to the Besel filter for my simulation. I guess that's something you'll have to determine. I have no idea what your simulation is, nor what its accuracy requirements are, nor why it's using this Bessel filter in the first place, nor the types of signals that will be input to the filter. If you're using a tool like Simulink, you could, in principal, simulate the analog filter directly (with appropriate step size). If instead of simulation (which I assume to mean forward-in-time only), you really mean just analyziing some data, it might be more appropriate to directly design a discrete-time filter with the desired gain characteristic and then use filtfilt() for a zero-phase filter, as suggested by @Star Strider, or maybe use the lowpass() function, which "compensates for the delay introduced by the filter." Note that forward-only filtering with the Bessel filter is not zero-phase.
Also, the poles of the discrete-time filter are very near the unit circle, as expected for such a large fs, with potential for round-off errors in computations to cause problems. For implementation, you might want to consider alternatives to the zpk form if that becomes a problem.
Amy Lg
Amy Lg on 6 Apr 2022
Thank you very much for your detailed answer and helping me to understand this.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!