0
$\begingroup$

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[0], 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?

$\endgroup$
1
$\begingroup$

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.

From the documentation

Compute an initial state zi for the sosfilt function that corresponds to the steady state of the step response.

It assumes that the input signal is a unit step. That's useful if you input signal has a large DC offset (at least in the beginning). You can call sosfilt_zi(), multiply it with the DC offset and apply it as state.

What about a multidimensional array filtered along an axis?

You need to determine the initial DC offset or bias of each channel individually and cascade the results of sosfilt_zi()*bias for each channel.

$\endgroup$
5
  • $\begingroup$ Ok, so this is completely useless in the case of a signal which does not look like a unit test. In the case of a more general signal, I should keep using a zero-initialized zi? $\endgroup$ – Mathieu Jun 4 at 12:13
  • $\begingroup$ No. Initialize with the mean of of the first N samples, where N is roughly the length of the impulse response of the filter or where the bias is reasonably constant. $\endgroup$ – Hilmar Jun 4 at 16:12
  • $\begingroup$ @Mathieu In the case of a more general signal, then you should use that signal to generate zi, no? $\endgroup$ – endolith Jun 6 at 5:22
  • 1
    $\begingroup$ @endolith That is what I thought too. Hilmar gave me a second answer on a different related question where I show a raw signal with a very large DC bias. zi can be multiplied with the DC bias to get the correct zi for my signal. $\endgroup$ – Mathieu Jun 6 at 10:03
  • 1
    $\begingroup$ @endolith the trick here is that you need to make a guess what the signal was earlier, so you need to extrapolate backwards in time. If you don't specify zi the filter function will assume that the signal has been zero before the first sample. This may or may not be a good assumption. In case of a large DC bias it's clearly NOT a good assumption. $\endgroup$ – Hilmar Jun 6 at 23:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.