Introduction
Cutoff frequency for a zero-phase-shift (zero-lag) IIR filter should be corrected due to damping of the cutoff frequency caused by multiple filter passes. This also applies to MATLAB's filtfilt function but is not directly mentioned in the documentation.
It seems that few people are aware of this correction, while people who are aware often use an incorrect implementation. Therefore, I chose to publish my findings regarding this issue to make people aware of it and to point to the correct reference. I have used MATLAB for the examples, but it applies to all implementation.
Theory
As described by David A. Winter in his book, Biomechanics and Motor Control of Human Movement, 4th-edition (page 69), and by D. Gordon E. Robertson and colleagues in their book, Research Methods in Biomechanics, 2nd-edition (page 288), the cutoff frequency has to be corrected for a zero-phase-shift Butterworth filter. This is because to compensate for the phase-shift introduced by each uneven filter pass, the filter is run again in reverse time direction, and for each filter pass, the -3db cutoff frequency will be pushed lower.
Winter shows how to compensate a 2nd-order Butterworth filter for this damping of cutoff frequency by using a correction factor that is applied on the adjusted Angular cutoff frequency (not directly on Angular cutoff frequency). The variable name Winter used for the adjusted Angular cutoff frequency is a bit misleading, therefore I used variable names similar to the names Robertson et al. used (only Ω replaced by U).
C = (((2^(1/filter_passes))-1)^(1/4)); % David A. Winter Butterworth correction factor
Wc = 2*pi*cutoff_frequency; % angular cutoff frequency
Uc = tan(Wc/(2*sampling_frequency)); % adjusted angular cutoff frequency
Un = Uc / C; % David A. Winter correction
Please note that older book editions and papers from both of these authors might not use this updated correction method.
Background
If you search the internet for this cutoff adjustment, you will find that many have implemented this method incorrectly. Instead of applying the correction factor to the adjusted angular cutoff frequency (Uc), it has been applied on the cutoff frequency. For example, there are a couple of cases where the correction has been done like this for a dual pass filter: f_corrected = f_cutoff / 0.802. This will give the incorrect response, though it’s not very far off at lower frequencies.
To investigate how a correction factor would look like if applied directly on the cutoff frequency, the cutoff frequency can be derived from the corrected parameter Un (see Results section how that can be done). It turns out that it is frequency dependent. Let’s call this frequency dependent factor Cf = f_cutoff / f_corrected. If anyone knows how to get this factor directly, please let me know.
Results
The correct way, according to Winter, to get the 2nd-order low-pass Butterworth coefficients are:
C = (((2^(1/filter_passes))-1)^(1/4)); % David A. Winter butterworth correction factor
Wc = 2*pi*cutoff_frequency; % angular cutoff frequency
Uc = tan(Wc/(2*sampling_frequency)); % adjusted angular cutoff frequency
Un = Uc / C; % David A. Winter correction
k1 = sqrt(2) * Un;
k2 = Un^2;
a0 = k2 / (1 + k1 + k2);
a1 = 2 * a0;
a2 = a0;
k3 = a1 / k2;
b1 = k3 - a1;
b2 = 1 - a1 - k3;
b = [ a0 a1 a2 ];
a = [ 1 -b1 -b2 ];
data_filtered = myfilter(data,b,a,'low',filter_passes);
If you want to use an existing function that only takes cutoff frequency as input (like MATLAB butter(..) function), then you can transform the corrected parameter Un back to frequency.
filter_passes = 2; % filtfilt use 2 passes
C = (((2^(1/filter_passes))-1)^(1/4)); % David A. Winter butterworth correction factor
Wc = 2*pi*cutoff_frequency; % angular cutoff frequency
Uc = tan(Wc/(2*sampling_frequency)); % adjusted angular cutoff frequency
Un = Uc / C; % David A. Winter correction
f_corrected = atan(Un)*sampling_frequency/pi; % corrected cutoff frequency
[b, a] = butter(2, f_corrected / (sampling_frequency / 2), 'low');
data_filtered = filtfilt(b, a, data);
Now we can test the difference between a reference single-pass filter, an uncorrected multi-pass filter, the INCORRECT frequency corrected multi-pass filter, and Winter’s corrected multi-pass filter. To test what happens with multiple passes, a filter with number of filter-passes as input is used (see attached files). The graphs shown is the spectrum of filtered noise.
With two passes, the uncorrected cutoff becomes 134Hz instead of 150Hz, the wrongly corrected cutoff becomes too high, and the corrected cutoff that Winter presented is spot on as shown in the figure below:
With a higher number of passes, the difference becomes bigger (at 4 passes, the wrongly corrected cutoff reaches fs/2 already at fs/3).
The same principle also applies to high-pass filters.
To get the correct 2nd-order high-pass Butterworth coefficients, you can use the method presented by Robertson et al:
cutoff_frequency = sampling_frequency/2 - cutoff_frequency; % Robertson high-pass modification
C = (((2^(1/filter_passes))-1)^(1/4)); % David A. Winter butterworth correction factor
Wc = 2*pi*cutoff_frequency; % angular cutoff frequency
Uc = tan(Wc/(2*sampling_frequency)); % adjusted angular cutoff frequency
Un = Uc / C; % David A. Winter correction
k1 = sqrt(2) * Un;
k2 = Un^2;
a0 = k2 / (1 + k1 + k2);
a1 = 2 * a0;
a2 = a0;
k3 = a1 / k2;
b1 = k3 - a1;
b2 = 1 - a1 - k3;
% Robertson high-pass modification
a1 = -a1;
b1 = -b1;
b = [ a0 a1 a2 ];
a = [ 1 -b1 -b2 ];
data_filtered = myfilter(data,b,a,'high',filter_passes);
You can get the corrected high-pass cutoff frequency in a similar way as for low-pass by multiplying with the correction factor (instead of dividing with it).
filter_passes = 2; % filtfilt use 2 passes
C = (((2^(1/filter_passes))-1)^(1/4)); % David A. Winter butterworth correction factor
Wc = 2*pi*cutoff_frequency; % angular cutoff frequency
Uc = tan(Wc/(2*sampling_frequency)); % adjusted angular cutoff frequency
Un = Uc * C; % multiply by C for highpass correction
f_corrected = atan(Un)*sampling_frequency/pi; % corrected cutoff frequency
[b, a] = butter(2, f_corrected / (sampling_frequency / 2), 'high');
data_filtered = filtfilt(b, a, data);
If we look at the difference for high-pass filters, we see that the incorrect frequency correction is not as bad as for low-pass filter, but it gets worse at higher cutoff frequencies and more filter passes:
Discussion
It is important to note that this specific correction factor (C) only applies to 2nd-order Butterworth coefficients, not for higher orders. Higher order filters will not be affected to the same degree from this damping of cutoff frequency due to a sharper cutoff curve. It is, however, still affected and needs to be considered. If anyone knows a correction factor that takes into account the filter order, please let me know.
History
