Let's take as an example an IIR filter using the 'sos' output, as this is what I use the most. To apply a bandpass filter, you can do:
from scipy.signal import butter, sosfilt fs = 512 # Sampling rate, Hz lowcut = 1 # Hz, cutting frequency highcut = 40 # Hz, currting frequency low = lowcut / (0.5 * fs) high = highcut / (0.5 * fs) sos = butter(2, [low, high], btype='band', output='sos', fs=fs) data = sosfilt(sos, data)
This is great. Now let's consider that our data is not a simple 1D example, but instead a 2D array of
N channels sampled during 1 second. Thus, with
data.shape = samples x channels the last line becomes:
data = sosfilt(sos, data, axis=0)
Now, what if we want to take into account the initial condition? Then, the
zi argument must be provided. The documentation says:
Initial conditions for the cascaded filter delays. It is a (at least 2D) vector of shape (n_sections, ..., 2, ...), where ..., 2, ... denotes the shape of x, but with x.shape[axis] replaced by 2. If zi is None or is not given then initial rest (i.e. all zeros) is assumed.
Thus, if we want to reproduce the 'None' behavior, using an initial rest, and returning the
zf (final filter decay value),
zi can be initialized with
np.zeros((sos.shape, 2, N)) for the data example above of shape
samples x channels.
zf: ndarray, optional
If zi is None, this is not returned, otherwise, zf holds the final filter delay values.
All this makes sense. But then what about a non-resting-state initial condition?
That is where the function
scipy.signal.sosfilt_zi() should come into play. But I don't understand what it does, and how it determines initial conditions from the
sos argument and not from the actual signal to filter.
Moreover, it returns a
zi of shape
(n_sections, 2) which is only good for the 1D case. What about a multidimensional array filtered along an axis? How should the
zi returned by
scipy.signal.sosfilt_zi() be used to create a
(n_sections, ..., 2, ...) array?