Repository: pyfar/pyfar

Issue:
Find IR delay - Polyfit Error
I was using the function `pyfar.dsp.find_impulse_response_delay()` and stumbled upon an error that I cannot seem to pinpoint exactly. The function works fine for impulse responses that are very similar to the one I attached to this issue, but for this one, the polynomial fit fails. Maybe you can understand why exactly this IR seems to raise an issue.

Minimal code example
------
```python
import pyfar as pf
import os

s = pf.io.read_audio(os.path.join(os.path.dirname(__file__), "ir_problem.flac"))
pf.dsp.find_impulse_response_delay(s)
```

Error message
---------
File "check_ir.py", line 5, in <module>
    pf.dsp.find_impulse_response_delay(s)
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^
  File "/opt/miniconda3/envs/playground/lib/python3.13/site-packages/pyfar/dsp/dsp.py", line 1325, in find_impulse_response_delay
    search_region_poly = np.polyfit(
        search_region_range[mask]-argmax, search_region[mask], N)
  File "/opt/miniconda3/envs/playground/lib/python3.13/site-packages/numpy/lib/_polynomial_impl.py", line 636, in polyfit
    raise TypeError("expected non-empty vector for x")
TypeError: expected non-empty vector for x

Environment
--------
Python: 3.13.1
pyfar: 0.7.1
Macbook Pro M1

[ir_problem.flac.zip](https://github.com/user-attachments/files/18966566/ir_problem.flac.zip)

Discussion:
@f-brinkmann: Hi, thanks for reporting. I had a weird intuition and tried `pf.dsp.find_impulse_response_delay(s * -1)` after seing that the largest peak has a negative value. This worked and could be used as quick fix, but it would be good to understand it in more detail.

Do the other IRs where you mentioned it worked also have larger negative peakes. @mberz if I remeber correctly you implemented this. Do you maybe have an idea. We could otherwise ask Anton to take a look.@hoyer-a: I actually ended up with the same quick solution when looking into this yesterday.
The implementation doesn't seem to like negative peaks in general (`warning for poorly conditioned Polyfit` or return `NaNs`), but i'm not familiar with the method used in the code, so maybe Marco can explain.@nico-franco-gomez: Oh ok, based on your comments, I think I found the error. Since the peak of the IR is phase-inverted, the cross-correlation with its minimum-phase equivalent will have a negative peak. This peak is successfully retrieved in line 1316 with `argmax = np.argmax(np.abs(correlation_analytic))`. Nonetheless, the analytic part has the opposite behavior of what is expected, leading in line 1322 to masking all values since the correlation is phase-inverted. A fix would be to multiply with -1 the analytic part of the search region in case the peak of the correlation function has been found to be a minimum.

My quick fix would be  add a line after 1318:
```python
if correlation_analytic[argmax].real < 0.0:
    search_region *= -1.
```

This seems to work just fine for my use cases. Should I write a PR or will you implement this?@f-brinkmann: Thanks for looking into this. I will create a pull and also add a test for this.@f-brinkmann: @nico-franco-gomez can you maybe check if you agree with the explanation of the fix: https://github.com/pyfar/pyfar/pull/776@nico-franco-gomez: Seems fine to me. Thank you for the fix!

Solution patch:
diff --git a/pyfar/classes/filter.py b/pyfar/classes/filter.py
index 5bf09657..50f12839 100644
--- a/pyfar/classes/filter.py
+++ b/pyfar/classes/filter.py
@@ -146,16 +146,8 @@ class Filter(object):
         else:
             self._coefficients = None
 
-        if state is not None:
-            if coefficients is None:
-                raise ValueError(
-                    "Cannot set a state without filter coefficients")
-            state = _atleast_3d_first_dim(state)
-            self._initialized = True
-        else:
-            self._initialized = False
+        self.state = state
 
-        self._state = state
         self._sampling_rate = sampling_rate
         self.comment = comment
 
@@ -208,12 +200,22 @@ class Filter(object):
 
     @property
     def state(self):
-        """
-        The current state of the filter as an array with dimensions
-        corresponding to the order of the filter and number of filter channels.
-        """
+        """Get or set the filter state."""
         return self._state
 
+    @state.setter
+    def state(self, state):
+        """Set state of the filter."""
+        if state is not None:
+            if self.coefficients is None:
+                raise ValueError(
+                    "Cannot set a state without filter coefficients")
+            self._initialized = True
+        else:
+            self._initialized = False
+
+        self._state = state
+
     @staticmethod
     def _process(coefficients, data, zi=None):
         raise NotImplementedError("Abstract class method.")
@@ -345,6 +347,13 @@ class FilterFIR(Filter):
 
     def __init__(self, coefficients, sampling_rate, state=None, comment=""):
 
+        if state is not None and np.asarray(state).ndim < 3:
+            state = _atleast_3d_first_dim(state)
+            warnings.warn(
+                "The state should be of shape (n_channels, *cshape, order). "
+                f"The state has been reshaped to shape=({state.shape}",
+                stacklevel=2)
+
         super().__init__(coefficients, sampling_rate, state, comment)
 
     @property
@@ -378,6 +387,48 @@ class FilterFIR(Filter):
         coeff = np.stack((b, a), axis=-2)
         self._coefficients = _atleast_3d_first_dim(coeff)
 
+    @property
+    def state(self):
+        """
+        Get or set the filter state.
+
+        The filter state sets the initial condition of the filter and can for
+        example be used for block-wise filtering.
+
+        Note that the state can also be initialized with the
+        :py:func:`~pyfar.classes.filter.FilterFIR.init_state` method.
+
+        Shape of the state must be  ``(n_filter_channels, *cshape, order)``,
+        with ``cshape`` being the channel shape of the
+        :py:class:`~pyfar.Signal` to be filtered.
+        """
+        return self._state
+
+    @state.setter
+    def state(self, state):
+        """
+        Set initial state of FIR filter.
+
+        Parameters
+        ----------
+        state : array, double
+            The state of the filter with dimensions
+            ``(n_channels, *cshape, order)``, where ``cshape`` is
+            the channel shape of the :py:class:`~pyfar.Signal`
+            to be filtered.
+        """
+        if state is not None:
+            state = np.asarray(state)
+            if (state.shape[-1] != self.order
+                or state.ndim < 3
+                or state.shape[0] != self.n_channels):
+                raise ValueError(
+                    "The state does not match the filter structure. Required "
+                    "shape for FilterFIR is (n_channels, *cshape, order).")
+
+        # sets filter state in parent class
+        Filter.state.fset(self, state)
+
     def init_state(self, cshape, state='zeros'):
         """Initialize the buffer elements to pre-defined initial conditions.
 
@@ -405,6 +456,11 @@ class FilterFIR(Filter):
         This is a hidden static method required for a shared processing
         function in the parent class.
         """
+        if zi is not None and zi.shape[0:-1] != data.shape[0:-1]:
+            raise ValueError("The initial state does not match the cshape of "
+                             "the signal. Required shape for `state` in "
+                             "FilterFIR is (n_channels, *cshape, order).")
+
         return spsignal.lfilter(coefficients[0], 1, data, zi=zi)
 
     def __repr__(self):
@@ -445,6 +501,13 @@ class FilterIIR(Filter):
 
     def __init__(self, coefficients, sampling_rate, state=None, comment=""):
 
+        if state is not None and np.asarray(state).ndim < 3:
+            state = _atleast_3d_first_dim(state)
+            warnings.warn(
+                "The state should be of shape (n_channels, *cshape, order). "
+                f"The state has been reshaped to shape=({state.shape}",
+                stacklevel=2)
+
         super().__init__(coefficients, sampling_rate, state, comment)
 
     @property
@@ -452,6 +515,48 @@ class FilterIIR(Filter):
         """The order of the filter."""
         return np.max(self._coefficients.shape[-2:]) - 1
 
+    @property
+    def state(self):
+        """
+        Get or set the filter state.
+
+        The filter state sets the initial condition of the filter and can for
+        example be used for block-wise filtering.
+
+        Note that the state can also be initialized with the
+        :py:func:`~pyfar.classes.filter.FilterIIR.init_state` method.
+
+        Shape of the state must be  ``(n_filter_channels, *cshape, order)``,
+        with ``cshape`` being the channel shape of the
+        :py:class:`~pyfar.Signal` to be filtered.
+        """
+        return self._state
+
+    @state.setter
+    def state(self, state):
+        """
+        Set initial state of IIR filter.
+
+        Parameters
+        ----------
+        state : array, double
+            The state of the filter with dimensions
+            ``(n_channels, *cshape, order)``, where ``cshape`` is
+            the channel shape of the :py:class:`~pyfar.Signal`
+            to be filtered.
+        """
+        if state is not None:
+            state = np.asarray(state)
+            if (state.shape[-1] != self.order
+                or state.ndim < 3
+                or state.shape[0] != self.n_channels):
+                raise ValueError(
+                    "The state does not match the filter structure. Required "
+                    "shape for FilterIIR is (n_channels, *cshape, order).")
+
+        # sets filter state in parent class
+        Filter.state.fset(self, state)
+
     def init_state(self, cshape, state='zeros'):
         """Initialize the buffer elements to pre-defined initial conditions.
 
@@ -479,6 +584,10 @@ class FilterIIR(Filter):
         This is a hidden static method required for a shared processing
         function in the parent class.
         """
+        if zi is not None and zi.shape[0:-1] != data.shape[0:-1]:
+            raise ValueError("The initial state does not match the cshape of "
+                             "the signal. Required shape for `state` in "
+                             "FilterIIR is (n_channels, *cshape, order).")
         return spsignal.lfilter(coefficients[0], coefficients[1], data, zi=zi)
 
     def __repr__(self):
@@ -519,8 +628,12 @@ class FilterSOS(Filter):
 
     def __init__(self, coefficients, sampling_rate, state=None, comment=""):
 
-        if state is not None:
+        if state is not None and np.asarray(state).ndim < 4:
             state = _atleast_4d_first_dim(state)
+            warnings.warn(
+                "The state should be of shape (n_channels, *cshape, n_sections"
+                f", 2). The state has been reshaped to shape=({state.shape}",
+                stacklevel=2)
 
         super().__init__(coefficients, sampling_rate, state, comment)
 
@@ -548,6 +661,51 @@ class FilterSOS(Filter):
         """The number of sections."""
         return self._coefficients.shape[-2]
 
+    @property
+    def state(self):
+        """
+        Get or set the filter state.
+
+        The filter state sets the initial condition of the filter and can for
+        example be used for block-wise filtering.
+
+        Note that the state can also be initialized with the
+        :py:func:`~pyfar.classes.filter.FilterSOS.init_state` method.
+
+        Shape of the state must be
+        ``(n_filter_channels, *cshape, n_sections, 2))``,
+        with ``cshape`` being the channel shape of the
+        :py:class:`~pyfar.Signal` to be filtered.
+        """
+        return self._state
+
+    @state.setter
+    def state(self, state):
+        """
+        Set initial state of SOS filter.
+
+        Parameters
+        ----------
+        state : array, double
+            The state of the filter with dimensions
+            ``(n_channels, *cshape, n_sections, 2)``, where ``cshape`` is
+            the channel shape of the :py:class:`~pyfar.Signal`
+            to be filtered.
+        """
+        if state is not None:
+            state = np.asarray(state)
+            if (state.shape[0] != self.n_channels
+                or state.shape[-1] != 2
+                or state.shape[-2] != self.n_sections
+                or state.ndim < 4):
+                raise ValueError(
+                    "The state does not match the filter structure. Required "
+                    "shape for FilterSOS is (n_channels, *cshape, "
+                    "n_sections, 2).")
+
+        # sets filter state in parent class
+        Filter.state.fset(self, state)
+
     def init_state(self, cshape, state='zeros'):
         """Initialize the buffer elements to pre-defined initial conditions.
 
@@ -575,6 +733,11 @@ class FilterSOS(Filter):
         This is a hidden static method required for a shared processing
         function in the parent class.
         """
+        if zi is not None and zi.shape[0:-2] != data.shape[0:-1]:
+            raise ValueError("The initial state does not match the cshape of "
+                             "the signal. Required shape for `state` in "
+                             "FilterSOS is (n_channels, *cshape, n_sections,"
+                             " 2).")
         if zi is not None:
             zi = zi.transpose(1, 0, 2)
         res = spsignal.sosfilt(sos, data, zi=zi, axis=-1)
diff --git a/pyfar/dsp/dsp.py b/pyfar/dsp/dsp.py
index 0aabd899..3ab8ee07 100644
--- a/pyfar/dsp/dsp.py
+++ b/pyfar/dsp/dsp.py
@@ -1292,7 +1292,8 @@ def find_impulse_response_delay(impulse_response, N=1):
             ir = impulse_response.time[ch]
             ir = np.real(ir) if mode == 'real' else np.imag(ir)
 
-            if np.max(ir) > 1e-16:
+            #Check absolute amximum, because peaks can be positive or negative.
+            if np.max(np.abs(ir)) > 1e-16:
                 # minimum phase warns if the input signal is not symmetric,
                 # which is not critical for this application
                 with warnings.catch_warnings():
@@ -1318,6 +1319,16 @@ def find_impulse_response_delay(impulse_response, N=1):
                 search_region = np.imag(
                     correlation_analytic[search_region_range])
 
+                # If this is true, it indicates that `correlation_analytic` has
+                # a negative peak, which can happen if the absolute maximum of
+                # `impulse_response` is negative. Changing the sign of the
+                # search region makes sure that the `mask` generated below
+                # works as intended. Fixing it this way is safer because it is
+                # theoretically possible that the absolute maximum of
+                # `impulse_response` is negative but
+                # `correlation_analytic[argmax].real` is positive.
+                search_region *= np.sign(correlation_analytic[argmax].real)
+
                 # mask values with a negative gradient
                 mask = np.gradient(search_region, search_region_range) > 0
 


Test patch:
diff --git a/tests/test_dsp.py b/tests/test_dsp.py
index e8a18f4d..98715054 100644
--- a/tests/test_dsp.py
+++ b/tests/test_dsp.py
@@ -764,7 +764,8 @@ def test_impulse_response_delay(is_complex):
     npt.assert_allclose(start_sample_est, start_sample, atol=1e-2)
 
 
-def test_impulse_response_delay_sinc():
+@pytest.mark.parametrize('sign', [1, -1])
+def test_impulse_response_delay_sinc(sign):
     """Test delay of a band-limited sinc function shifted by 1/2 samples."""
     sr = 44100
     n_samples = 128
@@ -774,6 +775,10 @@ def test_impulse_response_delay_sinc():
     sinc = np.sinc(samples - delay_samples)
     win = sgn.get_window('hann', n_samples, fftbins=False)
 
+    # Impulse responses with large negative peaks require special handling in
+    # `find_impulse_response_delay`. This makes tests the special case.
+    sinc *= sign
+
     ir = pf.Signal(sinc*win, sr)
     start_samples = pf.dsp.find_impulse_response_delay(ir)
     npt.assert_allclose(start_samples, delay_samples, atol=1e-3, rtol=1e-4)
diff --git a/tests/test_filter_classes.py b/tests/test_filter_classes.py
index c646af5f..4363b022 100644
--- a/tests/test_filter_classes.py
+++ b/tests/test_filter_classes.py
@@ -3,6 +3,7 @@ import numpy as np
 import numpy.testing as npt
 import pyfar.classes.filter as fo
 import pyfar as pf
+import re
 from scipy import signal as spsignal
 
 
@@ -34,6 +35,61 @@ def test_filter_init_with_state():
     npt.assert_array_equal(filt.state, state)
 
 
+def test_filter_state_setter():
+    coeff = np.array([[[1, 0, 0], [1, 0, 0]]])
+    state = np.array([[[1, 0]]])
+    filt = fo.Filter(coefficients=coeff)
+    filt.state = state
+    npt.assert_array_equal(filt.state, state)
+
+
+def test_filter_state_setter_errors():
+    """Test errors when setting a state with wrong shape."""
+    fs = 44100
+    coeff_iir = np.array([[[1, 0, 0], [1, 0, 0]]])
+    with pytest.raises(ValueError, match=re.escape(
+        "The state does not match the filter structure. Required "
+        "shape for FilterIIR is (n_channels, *cshape, order).")):
+        fo.FilterIIR(coeff_iir, fs, state=[[[1, 0, 0]]])
+
+    coeff_fir = np.array([1, 0, 0])
+    with pytest.raises(ValueError, match=re.escape(
+        "The state does not match the filter structure. Required "
+        "shape for FilterFIR is (n_channels, *cshape, order).")):
+        fo.FilterFIR(coeff_fir, fs, state=[[[1, 0, 0]]])
+
+    coeff_sos = [[[1, .5, .25, 1, .5, .25], [1, .5, .25, 1, .5, .25]]]
+    with pytest.raises(ValueError, match=re.escape(
+        "The state does not match the filter structure. Required shape for "
+        "FilterSOS is (n_channels, *cshape, n_sections, 2).")):
+        fo.FilterSOS(coeff_sos, fs, state=[[[[0], [1]]]])
+
+
+def test_filter_state_process_errors():
+    """Test errors when processing a signal with wrong state shape."""
+    fs = 44100
+    impulse = pf.signals.impulse(256, amplitude=np.ones((4, 1)))
+    coeff_iir = np.array([[[1, 0, 0], [1, 0, 0]]])
+    with pytest.raises(ValueError, match=re.escape(
+        "The initial state does not match the cshape of the signal. Required "
+        "shape for `state` in FilterIIR is (n_channels, *cshape, order).")):
+        fo.FilterIIR(coeff_iir, fs, state=[[[1, 0]]]).process(impulse)
+
+    coeff_fir = np.array([1, 0, 0])
+    with pytest.raises(ValueError, match=re.escape(
+        "The initial state does not match the cshape of the signal. Required "
+        "shape for `state` in FilterFIR is (n_channels, *cshape, order).")):
+        fo.FilterFIR(coeff_fir, fs, state=[[[1, 0]]]).process(impulse)
+
+    coeff_sos = [[[1, .5, .25, 1, .5, .25], [1, .5, .25, 1, .5, .25]]]
+    with pytest.raises(ValueError, match=re.escape(
+        "The initial state does not match the cshape of the signal. Required "
+        "shape for `state` in FilterSOS is (n_channels, *cshape, n_sections,"
+        " 2).")):
+        fo.FilterSOS(
+            coeff_sos, fs, state=[[[[0, 0], [1, 0]]]]).process(impulse)
+
+
 def test_filter_comment():
     filt = fo.Filter(coefficients=None, state=None, comment='Bla')
     assert filt.comment == 'Bla'