using scipy.signal.sosfilt to filter a stereo signal with zi, zf

Question

I'm stuck trying to filter a two-channel signal using sosfilt() with the zi initial state. I need to use the zi parameter so I can carry the filter state over multiple buffers.

Omitting zi works as expected

import numpy as np
from scipy.signal import sosfilt

>>> x = np.arange(20, dtype=np.float32).reshape(2, -1)
>>> x
array([[ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.],
       [10., 11., 12., 13., 14., 15., 16., 17., 18., 19.]], dtype=float32)
>>> sos = np.array([2, 0, 0, 1, 0, 0], dtype=np.float32)  # a very simple filter :)
>>> sos
array([2., 0., 0., 1., 0., 0.], dtype=float32)
>>> y = sosfilt(sos, x)
>>> y
array([[ 0.,  2.,  4.,  6.,  8., 10., 12., 14., 16., 18.],
       [20., 22., 24., 26., 28., 30., 32., 34., 36., 38.]], dtype=float32)

... but what does sosfilt() need for zi?

The sosfilt() documentation for zi 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.

The filter has two sections, so n_sections = 2. I assume that x.shape[axis] = 2 (because I can't imagine that zi needs to know that the source is 10 frames long). So I tried a few things, but none of them worked:

>>> zi = np.zeros((2, 2, 2), dtype=np.float32)
>>> y, zf = sosfilt(sos, x, zi)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "...\_signaltools.py", line 4277, in sosfilt
    x_zi_shape[axis] = 2
TypeError: only integer scalar arrays can be converted to a scalar index
>>> zi = np.zeros((1, 2, 2), dtype=np.float32) # maybe it means one second-order section?
>>> y, zf = sosfilt(sos, x, zi)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "...\_signaltools.py", line 4277, in sosfilt
    x_zi_shape[axis] = 2
TypeError: only integer scalar arrays can be converted to a scalar index
>>> zi = np.zeros((2, 2), dtype=np.float32) # at least 2D?
>>> y, zf = sosfilt(sos, x, zi)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "...\_signaltools.py", line 4277, in sosfilt
    x_zi_shape[axis] = 2
TypeError: only integer scalar arrays can be converted to a scalar index

What am I missing?

Answer 1

After more fiddling, I figured it out:

1. need to specify axis=1
2. n_sections is the number of second-order sections, not filter order.

So if N is the number of input channels, the correct initialization for zi and the correct call to sosfilt() is:

zi = np.zeros((1, N, 2), dtype=np.float32)
y, zf = sosfilt(sos, x, zi=zi, axis=1)