Repository: pyfar/pyfar (python)

Issue:
Feature request: easy Filter to Signal conversion like filter.to_ir()
As a regular pyfar user, I often find the filter classes to be inconvenient. The main problems I have with filter objects are:
- They cannot be directly plotted using pf.plot functions
- so far, FilterFIR always uses direct convolution, which is almost always painfully slow compared to FFT/overlap-add convolutions
- They don't work with utilities like dsp.linear_phase()

As a consequence, I find myself writing code like the following in pretty much every script I write with pyfar:
```py
dirac = pf.signals.impulse(ir_length, filter.sampling_rate)
filter_response = filter.process(dirac)
# the actual plotting or processing I want to do with the filter
```

or for FilterFIR:
```py
filter_response = pf.Signal(filter.coefficients, filter.sampling_rate)
```

I discussed this with @f-brinkmann and we concluded that instead of changing all the plotting and utility functions to support filter objects, a convenience function for filters that converts them into a Signal object containing their impulse response would help will all of these concerns.
The method could look like this:
`filter.to_ir(target_length)`, where filter is the Filter object to convert and target_length the number of samples that the impulse response should have. 
If we want more explicit naming, it could also be called `filter.to_impulse_response(len)` or similar, but this would make it less convenient to use in one-liners for plotting etc. 

Pull request:
Filter impulse response
### Which issue(s) are closed by this pull request?

Closes #749

### Changes proposed in this pull request:

- add `impulse_response` method to `pyfar.Filter` class
- add testing thereof


Patch:
diff --git a/pyfar/classes/filter.py b/pyfar/classes/filter.py
index 5bf09657..f1137b02 100644
--- a/pyfar/classes/filter.py
+++ b/pyfar/classes/filter.py
@@ -276,6 +276,41 @@ class Filter(object):
 
         return filtered_signal
 
+    def impulse_response(self, n_samples):
+        """
+        Compute or approximate the impulse response of the filter.
+
+        See `impulse_response` methods in derived classes for more details.
+
+        Parameters
+        ----------
+        n_samples : int
+            Length in samples for which the impulse response is computed.
+
+        Returns
+        -------
+        impulse_response : Signal
+            The impulse response of the filter.
+        """
+
+        # track the state (better than copying the entire filter)
+        if self.state is not None:
+            state = self.state.copy()
+            self._state = None
+            reset = True
+        else:
+            reset = False
+
+        # get impulse response
+        impulse_response = self.process(pf.signals.impulse(
+            n_samples, sampling_rate=self.sampling_rate), reset=reset)
+
+        # reset the state if required
+        if reset:
+            self._state = state
+
+        return impulse_response
+
     def reset(self):
         """Reset the filter state by filling it with zeros."""
         if self._state is not None:
@@ -399,6 +434,35 @@ class FilterFIR(Filter):
                 new_state[idx, ...] = spsignal.lfilter_zi(coeff[0], coeff[1])
         super().init_state(state=new_state)
 
+    def impulse_response(self, n_samples=None):
+        """
+        Compute the finite impulse response of the filter.
+
+        Parameters
+        ----------
+        n_samples : int, optional
+            Length in samples for which the impulse response is computed. The
+            default is ``None`` in which case the length of the impulse
+            response is determined from the filter :py:func:`~FilterFIR.order`.
+            A warning is returned if ``n_samples`` is to short to compute the
+            entire impulse response.
+
+        Returns
+        -------
+        impulse_response : Signal
+            The impulse response of the filter.
+        """
+
+        # set or check the impulse response length
+        if n_samples is None:
+            n_samples = self.coefficients.shape[-1]
+        elif self.coefficients.shape[-1] > n_samples:
+            warnings.warn(
+                ('n_samples should be at least as long as the filter, '
+                 f'which is {self.coefficients.shape[-1]}'), stacklevel=1)
+
+        return super().impulse_response(n_samples)
+
     @staticmethod
     def _process(coefficients, data, zi=None):
         """Process a single filter channel.
@@ -473,6 +537,27 @@ class FilterIIR(Filter):
                 new_state[idx, ...] = spsignal.lfilter_zi(coeff[0], coeff[1])
         return super().init_state(state=new_state)
 
+    def impulse_response(self, n_samples):
+        """
+        Approximate the infinite impulse response of the filter by a finite
+        impulse response.
+
+        Note that the number of samples must be sufficiently long for
+        `impulse_response` to be a good approximation of the theoretically
+        infinitely long impulse response of the filter.
+
+        Parameters
+        ----------
+        n_samples : int
+            Length in samples for which the impulse response is computed.
+
+        Returns
+        -------
+        impulse_response : Signal
+            The impulse response of the filter.
+        """
+        return super().impulse_response(n_samples)
+
     @staticmethod
     def _process(coefficients, data, zi=None):
         """Process a single filter channel.
@@ -569,6 +654,27 @@ class FilterSOS(Filter):
                 new_state[idx, ...] = spsignal.sosfilt_zi(coeff)
         return super().init_state(state=new_state)
 
+    def impulse_response(self, n_samples):
+        """
+        Approximate the infinite impulse response of the filter by a finite
+        impulse response.
+
+        Note that the number of samples must be sufficiently long for
+        `impulse_response` to be a good approximation of the theoretically
+        infinitely long impulse response of the filter.
+
+        Parameters
+        ----------
+        n_samples : int
+            Length in samples for which the impulse response is computed.
+
+        Returns
+        -------
+        impulse_response : Signal
+            The impulse response of the filter.
+        """
+        return super().impulse_response(n_samples)
+
     @staticmethod
     def _process(sos, data, zi=None):
         """Process a single filter channel.


Test patch:
diff --git a/tests/test_filter_classes.py b/tests/test_filter_classes.py
index c646af5f..0b45f0bd 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
 
 
@@ -43,6 +44,65 @@ def test_filter_comment():
         pf.Signal([1, 2, 3], 44100, comment=[1, 2, 3])
 
 
+@pytest.mark.parametrize('filter_object', [
+    # FIR filter, single-channel
+    fo.FilterFIR(np.array([[1, -1, .1, -.1]]), 48000),
+    # FIR filter, multi-channel
+    fo.FilterFIR(np.array([[1, -1, .1, -.1], [2, -2, .2, -.2]]), 48000),
+    # IIR filter, single-channel
+    fo.FilterIIR(np.array([[[1, -1, .1, -.1], [1, .1, -.1, 0]]]), 48000),
+    # IIR filter, multi-channel
+    fo.FilterIIR(np.array([[[1, -1, .1, -.1], [1, .1, -.1, 0]],
+                           [[1, -2, .2, -.2], [1, .2, -.2, 0]]]), 48000),
+    # SOS filter, single-channel
+    fo.FilterSOS(np.array([[[1, -1, .1, 1, .1, -.1]]]), 48000),
+    # SOS filter, multi-channel
+    fo.FilterSOS(np.array([[[1, -1, .1, 1, .1, -.1]],
+                           [[2, -2, .2, 1, .2, -.2]]]), 48000),
+])
+def test_filter_impulse_response(filter_object):
+    """Test 'impulse_response' class methods for all Filter types."""
+
+    # test signal
+    n_samples = 256
+    impulse = pf.signals.impulse(n_samples, sampling_rate=48000)
+
+    # filter without state
+    impulse_response = filter_object.impulse_response(n_samples)
+    # test output
+    assert isinstance(impulse_response, pf.Signal)
+    assert impulse_response.n_samples == n_samples
+    npt.assert_equal(impulse_response, filter_object.process(impulse))
+    assert filter_object.state is None
+
+    # filter with state
+    filter_object.init_state(impulse.cshape, 'zeros')
+    state_before = filter_object.state.copy()
+    impulse_response_2 = filter_object.impulse_response(n_samples)
+    state_after = filter_object.state.copy()
+    # test output
+    assert isinstance(impulse_response, pf.Signal)
+    assert impulse_response.n_samples == n_samples
+    npt.assert_equal(impulse_response_2, impulse_response)
+    assert filter_object.state is not None
+    npt.assert_equal(state_after, state_before)
+
+
+def test_filter_fir_impulse_response_n_samples_warning():
+    """Test warning if impulse response length is less then filter length."""
+    filter_object = fo.FilterFIR([[1, 0, 0, 1]], 48000)
+    match = "n_samples should be at least as long as the filter, which is 4"
+    with pytest.warns(UserWarning, match=re.escape(match)):
+        _ = filter_object.impulse_response(2)
+
+
+def test_filter_fir_impulse_response_n_samples_default():
+    """Test default value for impulse response length."""
+    filter_object = fo.FilterFIR([[1, 0, 0, 1]], 48000)
+    impulse_response = filter_object.impulse_response()
+    assert impulse_response.n_samples == 4
+
+
 @pytest.mark.parametrize('filter_object', [
     (fo.FilterFIR([[1, -1]], 48000)),
     (fo.FilterIIR([[1, -1], [1, -1]], 48000)),
@@ -68,7 +128,6 @@ def test_filter_init_state_default_and_error(filter_object):
         filter_object.init_state((1, ), 'random')
 
 
-
 def test_filter_iir_init():
     coeff = np.array([[1, 1/2, 0], [1, 0, 0]])
     filt = fo.FilterIIR(coeff, sampling_rate=2*np.pi)