How can I normalize my filter?

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

[Star Strider]

It is generally not possible to convert a Bessel filter to a discrete version.

From the documentation:

‘The besself function does not support the design of digital Bessel filters.’

However if you want to see its performance characteristics, use the freqz function.

Compare it with the transfer function version of the analogue filter with freqs.

.

[Amy Lg]

Thank you for your reply. You mean the code is not correct? So what can I do if I need a Bessel filter performance?

I found from one of the posts here that we can do conversion from analog to digital by bilinear function.

[Star Strider]

I do not advise using a Bessel filter as a discrete filter for the reasons I previously stated.

‘So, is it wrong if I use bessel filter for my simulation?’

Yes.

Use any other design and then use filtfilt to do the actual filtering. since it will make all discrete filter realisations have a phase-neutral response.

EDIT — See the Digital section of the Wikipedia article on the Bessel filter for an appropriate discussion.

[Amy Lg]

Thank you so much for your help.

[Star Strider]

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!

Answer 1

[Paul]

Hi @Amy Lg

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.

[Star Strider]

Note my earlier Comment.

[Paul]

I thought that statement in the doc was just saying that besself() can only be used for continuous-time filters, whereas other filter design functions, like ellip() and butter(), can be used for either continuous- or discrete-time filters.

My understanding of the issue with using the bilinear() on a Bessel filter is that the blinear transformation does not perfectly preserve the approximately constant group delay of the continuous-time filter below the cut-off frequency. Whether or not the group delay of the resulting discrete-time filter is good enough is definitely a question to consider.

[Star Strider]

The discrete implementation does not preserve the zero-phase behaviour of the Bessel filter.

[Paul]

What is the definition of zero-phase in this context?

[Star Strider]

No phase distortion across the passband.

[Amy Lg]

Hi @Paul

Thank you so much for your explanation and help.

Regarding fs, your assumption is correct. I corrected my scaling in the code.

and regarding transfer function I meant that the transfer function of my filter should peaks at one (normalized) to make sure that it does not increase the energy of input.

Concerning phase distortion across the passband, I guess the phase response is close to linear over the passband based on your figure. So, is it wrong if I use bessel filter for my simulation?

[Paul]

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]

Thank you very much for your detailed answer and helping me to understand this.