Doc mention cutoff frequency correction for signal.filtfilt

[getzze]

The filtfilt function runs two filters to cancel the phase shift. Hence to compensate for the twice -3dB attenuation, one should define a higher cutoff frequency equal to f_c = f_wanted / 0.8 . See explanations in https://nbviewer.jupyter.org/github/demotu/BMC/blob/master/notebooks/DataFiltering.ipynb, Cell[8]
This should be mentioned in the doc and the examples changed appropriately for filtfilt and lfilter

Simple example:

import numpy as np
import scipy
import scipy.signal
import matplotlib.pyplot as plt

f0 = 0.1
N = 2
rip = 0.05

x = np.linspace(0,100, 1001)
y = np.random.randn(len(x))

cutoffs = [None, f0, f0, f0/0.8, f0*0.8]
lss = ["-", "--", ":", ":", ":"]
labels = ["raw", "single pass", "dp: no corr", "dp: good corr", "dp: bad corr"]
Y = np.vstack((y, ys, yf, yff, yfff))

fig,ax = plt.subplots(figsize=(16,12))
single_pass=True
for cutoff, label, ls in zip(cutoffs, labels, lss):
    if cutoff is None:
        yf = y
    else:
        b, a = scipy.signal.cheby1(N, rip, cutoff)
        if single_pass:
            yf = scipy.signal.lfilter(b, a, y)
            single_pass = False
        else:
            yf = scipy.signal.filtfilt(b, a, y)    
    ax.plot(x, yf, label=label, linestyle=ls, linewidth=2)

ax.set_xlabel("time")
ax.legend()
ax.set_xlim(2,5)
ax.set_ylim(-1.5,1.5)

The amplitude of the single pass lfilter (yellow) is comparable to the filtfilt with the corrected cutoff f0/0.8 (red).

On python3.7, Scipy 1.1

[getzze]

Incidently, I think there is a mistake in the decimate function.

In the source, the low-pass filter is defined with the cutoff frequency f0*0.8 instead of f0/0.8 :
(line 3434 in v1.1.0) system = dlti(*cheby1(n, 0.05, 0.8 / q)).

Here is a minimal example demonstrating the difference:

import numpy as np
import scipy.signal
import matplotlib.pyplot as plt

x = np.linspace(0, 100, 1001)
fs = 1/(x[1] - x[0]) # sampling frequency
f1 = 0.001*fs
y = 2*np.sin(2*np.pi*f1 * x) # sine with 1000 points per cycle
y += 0.5*np.sin(2*np.pi*10*f1 * x) # sine with 100 points per cycle

## Minimal scipy.signal.decimate
def good_dec(x, q, n=8):
    filtered = scipy.signal.filtfilt(*scipy.signal.cheby1(n, 0.05, 1 / (0.8 * q)), x)
    return filtered[::q]
def bad_dec(x, q, n=8):
    filtered = scipy.signal.filtfilt(*scipy.signal.cheby1(n, 0.05, 0.8 / q), x)
    return filtered[::q]

## Decimate
fact = 100  # only the first sine will survive
xd = x[::fact]
yc = good_dec(y, fact)
yb = bad_dec(y, fact)

fig,ax = plt.subplots(figsize=(16,12))
ax.plot(x, y, label="raw", linestyle="-", linewidth=2)
ax.plot(xd, yc, label="f0/0.8: r={}".format(fact), linestyle=":", linewidth=3)
ax.plot(xd, yb, label="f0*0.8: r={}".format(fact), linestyle="--", linewidth=3)

ax.set_xlabel("time")
ax.legend()
ax.set_xlim(0,100)

The decimation process with the correct frequency cutoff seems to give a better result than the current implementation.

I could not find any reference for the choice of f0*0.8, but the Matlab implementation is the same.

[aweinstein]

Please notice that the 0.8 correction is valid only for a Butterworth filter. Other filters should use different correction values.

[larsoner]

I could not find any reference for the choice of f0*0.8, but the Matlab implementation is the same.

I see two possibilities. One is that it's a typo. The other is that it's not related to the filtfilt vs lfilter, but instead a slightly more conservative choice for the lowpass corner frequency to reduce aliasing. However, given that the cutoff freq in the lfilter/fir code path of decimate is 1. / q and the cutoff freq in the filtfilt/iir code path is 0.8 / q it does seem more like a typo. The weird part is that, as @aweinstein mentions, is that Chebychev filters -- like those used in decimate -- could require a different constant. So this might actually be two bugs, one that should be the reciprocal of the correction value, and that the correction value 0.8 might not be correct (someone would need to check).

One option to move forward would be to open a PR and see if any existing tests break. Then extend the existing demos here to show that using the reciprocal of the correct value works better.

While we're at it, I'm not sure why we bother going through dlti when we could just get the TF coefficients directly.

And really we should probably use sosfiltfilt for better stability. Maybe not all in the same PR, though :)

[aweinstein]

@larsoner You can find a reference to the 0.8 correction factor in page 69 of this book:

https://books.google.com.br/books?id=_bFHL08IWfwC&printsec=frontcover&source=gbs_ge_summary_r&cad=0#v=onepage&q&f=false

This is a screenshot of the relevant paragraph:

Section 2.6 of this paper also discuss the effects of filtfilt (how one should report its use, etc.). Perhaps the documentation should point to this paper.

[Xurotash]

This site mention the correction factor also (but use it like this c_butter = 1/C ). And reference to this paper.

All references I seen is for lowpass filters (f_c_low = f_wanted / C). Now my quest is what about highpass filters? I have tested and found that by inverting C I get an rather good response (f_c_high = f_wanted * C), but is this correct?

Thanks for det book reference @aweinstein, does it state anything about highpass filters?

[getzze]

For a Butterworth high-pass filter, the correction factor would be f_c_low = C*f_wanted .
Indeed, a Chebyshev filter is used for the decimate function. Given that the slope is very steep, setting C=1 should be ok.

For filtfilt documentation, it would be good to emphasize that the cutoff frequency could need to be corrected and that the effective order of the filter is doubled (because two filters are applied).

[Xurotash]

I dug some deeper into this cutoff correction, and realised that the correction factor should NOT be applied directly on the cutoff frequency, i.e. f_c_low = f_wanted / C is incorrect. I find this so interesting that i wrote an article on codeproject to highlight this issue.
Multi-Pass Filter Cutoff Correction

[getzze]

@Xurotash So the CORRECT corrected cutoff frequency is:
f_corr = fs/pi * arctan( tan(pi*f_c/fs) / C)
with C = (2**(1/n) - 1)**(1/4) , n number of passes, fs the sampling frequency and f_c the wanted cutoff frequency.
Tested with scipy.signal.butter(2, f_corr/(fs/2)) it gives me the same result as rewriting butter with a correction factor.
When f_c << fs/2, the correct correction is almost equal to the wrong correction (f_c/C), so it was not an error easy to spot!

[piecot]

This issue is related to #9402 on signal attenuation in signal.decimate, isn't it?

[larsoner]

@piecot that one has to do with the DC component being affected by the Chebychev ripple, whereas this issue has to do with a frequency correction term for low-pass filter design when you know you will do two-pass filtering.

So to summarize the bugs:

1. The correction term is always applied, regardless of single or two-pass filtering being used (zero_phase)
2. The correction term is applied the opposite of how it should be (divide rather than multiply)
3. The correction term is for Butterworth filters, we use cheby1
4. The Chebychev filter (at least for even orders) does not treat the DC term properly (Signal attenuation in signal.decimate for even-order IIR filters #9402)

If someone can find the correction term for Chebychev filters, we can at least take care of points 1-3.

[roryyorke]

The IIR filter parameters in decimate seem to have been chosen to match Matlab's decimate; the documentation there in turn references a 1979 DSP book. My guess is the 0.8 is to improve anti-aliasing, and not a mistaken correction factor derived from 2nd order Butterworth filters.

As noted in #5392, Matlab's decimate offers only filtfilt behaviour. We could tweak decimate so that the forward response has a somewhat similar response to the forward-backward response, but I don't know if it's worth it.

To the original bug report: to design a cheby1 filter for filtfilt, I would halve the ripple, and maybe the order too, rather than fiddling with the cut-off frequency. I don't think the docs for filtfilt need to be updated for this; they already say The combined filter has zero phase and a filter order twice that of the original..