Add method="fir" to highpass_filter and bandpass_filter

doc/api.rst

@@ -188,6 +188,8 @@ spikeinterface.preprocessing
     .. autofunction:: save_motion_info
     .. autofunction:: decimate
     .. autofunction:: depth_order
+    .. autofunction:: design_fir_bandpass
+    .. autofunction:: design_fir_highpass
     .. autofunction:: detect_bad_channels
     .. autofunction:: detect_and_interpolate_bad_channels
     .. autofunction:: detect_and_remove_bad_channels
@@ -197,6 +199,7 @@ spikeinterface.preprocessing
     .. autofunction:: detect_and_remove_artifacts
     .. autofunction:: directional_derivative
     .. autofunction:: filter
+    .. autofunction:: fir_filter
     .. autofunction:: gaussian_filter
     .. autofunction:: highpass_filter
     .. autofunction:: highpass_spatial_filter

src/spikeinterface/core/testing.py

@@ -91,6 +91,35 @@ def check_recordings_equal(
         check_extractor_properties_equal(RX1, RX2)
 
 
+def check_chunk_invariance(
+    recording: BaseRecording, chunk_size: int, atol: float = 0.0, segment_index: int = 0
+) -> float:
+    """
+    Check that traces read in chunks of `chunk_size` and concatenated equal traces read in one call.
+
+    A preprocessor whose margin covers its full support should pass with `atol` at float rounding; one that
+    truncates an unbounded response, such as an IIR filter, should pass with `atol` at its stated seam error.
+
+    Returns
+    -------
+    max_error : float
+        The largest absolute difference between the chunked and the whole traces.
+    """
+    whole = recording.get_traces(segment_index=segment_index).astype(np.float64)
+    num_samples = whole.shape[0]
+    chunked = np.concatenate(
+        [
+            recording.get_traces(
+                segment_index=segment_index, start_frame=start, end_frame=min(start + chunk_size, num_samples)
+            )
+            for start in range(0, num_samples, chunk_size)
+        ]
+    )
+    max_error = float(np.abs(chunked - whole).max())
+    assert max_error <= atol, f"chunked traces differ from whole traces by {max_error}, above {atol}"
+    return max_error
+
+
 def check_sortings_equal(
     SX1: BaseSorting,
     SX2: BaseSorting,

src/spikeinterface/preprocessing/__init__.py

@@ -12,6 +12,7 @@
 from .preprocessing_tools import get_spatial_interpolation_kernel
 from .detect_bad_channels import detect_bad_channels
 from .correct_lsb import correct_lsb
+from .filter_design import design_fir_bandpass, design_fir_highpass
 
 from .pipeline import (
     apply_preprocessing_pipeline,

src/spikeinterface/preprocessing/filter.py

@@ -157,6 +157,32 @@ def __init__(
             direction=direction,
         )
 
+    def _init_fir(self, recording, taps, dtype, filter_kwargs):
+        """Set this recording up to apply the zero-phase FIR `taps` by overlap-save, for ``method="fir"``."""
+        from .fir_filter import FIRFilterRecordingSegment
+
+        iir_defaults = dict(filter_order=5, ftype="butter", filter_mode="sos", coeff=None, direction="forward-backward")
+        for name, value in filter_kwargs.items():
+            if name not in iir_defaults and name != "add_reflect_padding":
+                raise TypeError(f"Unexpected argument {name!r}")
+            default = iir_defaults.get(name)
+            if name in iir_defaults and (value is not None if default is None else value != default):
+                raise ValueError(f"{name} applies to method='iir' only; method='fir' is set by its transition bands")
+        dtype = fix_dtype(recording, dtype)
+
+        BasePreprocessor.__init__(self, recording, dtype=dtype)
+        self.annotate(is_filtered=True)
+        if "offset_to_uV" in self.get_property_keys():
+            self.set_channel_offsets(0)
+        self.margin_samples = taps.size // 2
+        kernels = taps[np.newaxis, :]
+        groups = np.zeros(recording.get_num_channels(), dtype=np.intp)
+        add_reflect_padding = filter_kwargs.get("add_reflect_padding", False)
+        for parent_segment in recording.segments:
+            self.add_recording_segment(
+                FIRFilterRecordingSegment(parent_segment, kernels, groups, taps.size // 2, dtype, add_reflect_padding)
+            )
+
 
 class FilterRecordingSegment(BasePreprocessorSegment):
     def __init__(
@@ -243,6 +269,16 @@ class BandpassFilterRecording(FilterRecording):
         The dtype of the returned traces. If None, the dtype of the parent recording is used
     ignore_low_freq_error : bool, default: False
         If True, does not raise an error if freq_min is too low for the sampling frequency.
+    method : "iir" | "fir", default: "iir"
+        "iir" applies a Butterworth (or `ftype`) filter with `scipy.signal`. "fir" applies a linear-phase Kaiser
+        bandpass from ``design_fir_bandpass`` by overlap-save: it is -6 dB at `freq_min` and `freq_max`, at least
+        `stopband_db` down beyond the transition bands, and its output does not depend on chunking. With "fir",
+        `margin_ms` is ignored, and `filter_order`, `ftype`, `filter_mode`, `coeff` and `direction` must be left
+        at their defaults.
+    transition_low_hz, transition_high_hz : float, default: 175.0, 2000.0
+        For `method="fir"`, the width of the transition band centred on each cutoff, in Hz.
+    stopband_db : float, default: 45.0
+        For `method="fir"`, the least attenuation in both stopbands, in dB.
     {}
 
     Returns
@@ -260,19 +296,32 @@ def __init__(
         dtype=None,
         ignore_low_freq_error=False,
         _skip_margin_warning_for_old_version=False,
+        method="iir",
+        transition_low_hz=175.0,
+        transition_high_hz=2000.0,
+        stopband_db=45.0,
         **filter_kwargs,
     ):
-        if margin_ms == "auto":
-            margin_ms = adjust_margin_ms_for_highpass(freq_min)
-        highpass_check(
-            freq_min,
-            margin_ms,
-            ignore_low_freq_error=ignore_low_freq_error,
-            skip_warning=_skip_margin_warning_for_old_version,
-        )
-        FilterRecording.__init__(
-            self, recording, band=[freq_min, freq_max], margin_ms=margin_ms, dtype=dtype, **filter_kwargs
-        )
+        if method == "iir":
+            if margin_ms == "auto":
+                margin_ms = adjust_margin_ms_for_highpass(freq_min)
+            highpass_check(
+                freq_min,
+                margin_ms,
+                ignore_low_freq_error=ignore_low_freq_error,
+                skip_warning=_skip_margin_warning_for_old_version,
+            )
+            FilterRecording.__init__(
+                self, recording, band=[freq_min, freq_max], margin_ms=margin_ms, dtype=dtype, **filter_kwargs
+            )
+        elif method == "fir":
+            from .filter_design import design_fir_bandpass
+
+            fs = recording.get_sampling_frequency()
+            taps = design_fir_bandpass(freq_min, freq_max, fs, transition_low_hz, transition_high_hz, stopband_db)
+            self._init_fir(recording, taps, dtype, filter_kwargs)
+        else:
+            raise ValueError(f"method must be 'iir' or 'fir', got {method!r}")
         dtype = fix_dtype(recording, dtype)
         self._kwargs = dict(
             recording=recording,
@@ -281,6 +330,10 @@ def __init__(
             margin_ms=margin_ms,
             dtype=dtype.str,
             ignore_low_freq_error=ignore_low_freq_error,
+            method=method,
+            transition_low_hz=transition_low_hz,
+            transition_high_hz=transition_high_hz,
+            stopband_db=stopband_db,
         )
         self._kwargs.update(filter_kwargs)
 
@@ -314,6 +367,15 @@ class HighpassFilterRecording(FilterRecording):
         The dtype of the returned traces. If None, the dtype of the parent recording is used
     ignore_low_freq_error : bool, default: False
         If True, does not raise an error if freq_min is too low for the sampling frequency.
+    method : "iir" | "fir", default: "iir"
+        "iir" applies a Butterworth (or `ftype`) filter with `scipy.signal`. "fir" applies a linear-phase Kaiser
+        highpass from ``design_fir_highpass`` by overlap-save: it is -6 dB at `freq_min`, at least `stopband_db`
+        down below the transition band, and its output does not depend on chunking. With "fir", `margin_ms` is
+        ignored, and `filter_order`, `ftype`, `filter_mode`, `coeff` and `direction` must be left at their defaults.
+    transition_low_hz : float, default: 175.0
+        For `method="fir"`, the width of the transition band centred on `freq_min`, in Hz.
+    stopband_db : float, default: 45.0
+        For `method="fir"`, the least attenuation in the stopband, in dB.
     {}
 
     Returns
@@ -330,26 +392,40 @@ def __init__(
         dtype=None,
         ignore_low_freq_error=False,
         _skip_margin_warning_for_old_version=False,
+        method="iir",
+        transition_low_hz=175.0,
+        stopband_db=45.0,
         **filter_kwargs,
     ):
-        if margin_ms == "auto":
-            margin_ms = adjust_margin_ms_for_highpass(freq_min)
-        highpass_check(
-            freq_min,
-            margin_ms,
-            ignore_low_freq_error=ignore_low_freq_error,
-            skip_warning=_skip_margin_warning_for_old_version,
-        )
-        FilterRecording.__init__(
-            self, recording, band=freq_min, margin_ms=margin_ms, dtype=dtype, btype="highpass", **filter_kwargs
-        )
+        if method == "iir":
+            if margin_ms == "auto":
+                margin_ms = adjust_margin_ms_for_highpass(freq_min)
+            highpass_check(
+                freq_min,
+                margin_ms,
+                ignore_low_freq_error=ignore_low_freq_error,
+                skip_warning=_skip_margin_warning_for_old_version,
+            )
+            FilterRecording.__init__(
+                self, recording, band=freq_min, margin_ms=margin_ms, dtype=dtype, btype="highpass", **filter_kwargs
+            )
+        elif method == "fir":
+            from .filter_design import design_fir_highpass
+
+            taps = design_fir_highpass(freq_min, recording.get_sampling_frequency(), transition_low_hz, stopband_db)
+            self._init_fir(recording, taps, dtype, filter_kwargs)
+        else:
+            raise ValueError(f"method must be 'iir' or 'fir', got {method!r}")
         dtype = fix_dtype(recording, dtype)
         self._kwargs = dict(
             recording=recording,
             freq_min=freq_min,
             margin_ms=margin_ms,
             dtype=dtype.str,
             ignore_low_freq_error=ignore_low_freq_error,
+            method=method,
+            transition_low_hz=transition_low_hz,
+            stopband_db=stopband_db,
         )
         self._kwargs.update(filter_kwargs)
 

src/spikeinterface/preprocessing/filter_design.py

@@ -0,0 +1,116 @@
+"""Filter design, separate from filter application, so kernels can be inspected and composed before use."""
+
+from __future__ import annotations
+
+import numpy as np
+
+
+def _stopband_db(taps, fs, edges):
+    """Smallest attenuation, in dB, of `taps` over the frequency intervals `edges`, sampled 16 times per ripple."""
+    import scipy.signal
+
+    worst = 0.0
+    for low, high in edges:
+        freqs = np.linspace(low, high, max(64, int(16 * (high - low) / fs * taps.size)))
+        _, response = scipy.signal.freqz(taps, worN=freqs, fs=fs)
+        worst = max(worst, np.abs(response).max())
+    return -20 * np.log10(max(worst, 1e-300))
+
+
+def design_fir_bandpass(freq_min, freq_max, fs, transition_low_hz=175.0, transition_high_hz=2000.0, stopband_db=45.0):
+    """
+    Design a linear-phase bandpass FIR kernel by the Kaiser window method.
+
+    The response is -6 dB at `freq_min` and `freq_max`, and at least `stopband_db` down below
+    ``freq_min - transition_low_hz / 2`` and above ``freq_max + transition_high_hz / 2``. Passband ripple is at
+    most about ``2 * 10 ** (-stopband_db / 20)`` in amplitude. The kernel convolves a highpass and a lowpass,
+    each windowed to its own transition: one window would give both edges the narrower width, and a sharp
+    upper edge rings for longer.
+
+    Parameters
+    ----------
+    freq_min, freq_max : float
+        Lower and upper -6 dB cutoffs, in Hz.
+    fs : float
+        Sampling frequency, in Hz.
+    transition_low_hz, transition_high_hz : float, default: 175.0, 2000.0
+        Width of the transition band centred on each cutoff, in Hz.
+    stopband_db : float, default: 45.0
+        Least attenuation in both stopbands, in dB.
+
+    Returns
+    -------
+    taps : np.ndarray
+        Symmetric, odd-length float64 kernel, for ``fir_filter``.
+    """
+    import scipy.signal
+
+    stop_low = freq_min - transition_low_hz / 2
+    stop_high = freq_max + transition_high_hz / 2
+    if transition_low_hz <= 0 or transition_high_hz <= 0 or stopband_db <= 0:
+        raise ValueError("transition_low_hz, transition_high_hz and stopband_db must be positive")
+    if not 0 < stop_low < freq_min < freq_max < stop_high < fs / 2:
+        raise ValueError(
+            f"The transition bands, [{stop_low}, {freq_min + transition_low_hz / 2}] and "
+            f"[{freq_max - transition_high_hz / 2}, {stop_high}] Hz, must lie in (0, {fs / 2}) Hz in that order"
+        )
+
+    def edge(cutoff, width, pass_zero):
+        num_taps, beta = scipy.signal.kaiserord(stopband_db, width / (fs / 2))
+        return num_taps | 1, beta, cutoff, pass_zero
+
+    edges = [edge(freq_min, transition_low_hz, False), edge(freq_max, transition_high_hz, True)]
+    stopbands = [(0.0, stop_low), (stop_high, fs / 2)]
+    while True:
+        highpass, lowpass = (
+            scipy.signal.firwin(n, cutoff, window=("kaiser", beta), pass_zero=pass_zero, fs=fs)
+            for n, beta, cutoff, pass_zero in edges
+        )
+        taps = np.convolve(highpass, lowpass)
+        # The two kernels' ripples can add, so the window's nominal attenuation is checked and the failing edge
+        # lengthened.
+        misses = [_stopband_db(taps, fs, [band]) < stopband_db for band in stopbands]
+        if not any(misses):
+            return taps
+        edges = [(n + 2 * miss, beta, cutoff, pz) for (n, beta, cutoff, pz), miss in zip(edges, misses)]
+
+
+def design_fir_highpass(freq_min, fs, transition_low_hz=175.0, stopband_db=45.0):
+    """
+    Design a linear-phase highpass FIR kernel by the Kaiser window method.
+
+    The response is -6 dB at `freq_min` and at least `stopband_db` down below ``freq_min - transition_low_hz / 2``,
+    with passband ripple of at most about ``10 ** (-stopband_db / 20)`` in amplitude.
+
+    Parameters
+    ----------
+    freq_min : float
+        The -6 dB cutoff, in Hz.
+    fs : float
+        Sampling frequency, in Hz.
+    transition_low_hz : float, default: 175.0
+        Width of the transition band centred on the cutoff, in Hz.
+    stopband_db : float, default: 45.0
+        Least attenuation in the stopband, in dB.
+
+    Returns
+    -------
+    taps : np.ndarray
+        Symmetric, odd-length float64 kernel, for ``fir_filter``.
+    """
+    import scipy.signal
+
+    stop = freq_min - transition_low_hz / 2
+    if transition_low_hz <= 0 or stopband_db <= 0:
+        raise ValueError("transition_low_hz and stopband_db must be positive")
+    if not 0 < stop < freq_min + transition_low_hz / 2 < fs / 2:
+        raise ValueError(
+            f"The transition band, [{stop}, {freq_min + transition_low_hz / 2}] Hz, must lie in (0, {fs / 2}) Hz"
+        )
+    num_taps, beta = scipy.signal.kaiserord(stopband_db, transition_low_hz / (fs / 2))
+    num_taps |= 1  # a highpass needs a type I kernel
+    while True:
+        taps = scipy.signal.firwin(num_taps, freq_min, window=("kaiser", beta), pass_zero=False, fs=fs)
+        if _stopband_db(taps, fs, [(0.0, stop)]) >= stopband_db:
+            return taps
+        num_taps += 2