From 92426274f4b22669ff1e8056cf269afe613aba16 Mon Sep 17 00:00:00 2001 From: Galen Lynch Date: Sun, 11 Oct 2026 14:35:03 -0700 Subject: [PATCH 1/3] feat: add fir_filter step and FIR bandpass and highpass designs fir_filter(recording, taps) applies any FIR kernel, shared or per channel, by overlap-save. Its margin is the kernel's support, so its output does not depend on chunking. Channels are filtered in blocks sized to keep each block's buffers in the L2 cache. design_fir_bandpass and design_fir_highpass design Kaiser-window kernels from a stated spec: -6 dB at the cutoffs, a given transition width at each edge, and a given stopband attenuation. Keeping design separate from application lets kernels be inspected or combined before use. The defaults (175 Hz lower transition, 45 dB stopband) pass at most 0.03% of the probe's specified noise power below 150 Hz, and no more than the default Butterworth between 150 and 300 Hz, on every channel of a Neuropixels 2.0 recording. A 400 Hz lower transition passed 2.4 times the Butterworth's power between 150 and 300 Hz. check_chunk_invariance in core.testing compares chunked reads with one whole read, so any preprocessor's seam error can be tested directly. --- doc/api.rst | 3 + src/spikeinterface/core/testing.py | 29 +++ src/spikeinterface/preprocessing/__init__.py | 1 + .../preprocessing/filter_design.py | 116 ++++++++++ .../preprocessing/fir_filter.py | 218 ++++++++++++++++++ .../preprocessing/preprocessing_classes.py | 2 + .../preprocessing/tests/test_fir_filter.py | 207 +++++++++++++++++ .../preprocessing/tests/test_pipeline.py | 3 + 8 files changed, 579 insertions(+) create mode 100644 src/spikeinterface/preprocessing/filter_design.py create mode 100644 src/spikeinterface/preprocessing/fir_filter.py create mode 100644 src/spikeinterface/preprocessing/tests/test_fir_filter.py diff --git a/doc/api.rst b/doc/api.rst index c0ace59d3b..38bbfc0e98 100755 --- a/doc/api.rst +++ b/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 diff --git a/src/spikeinterface/core/testing.py b/src/spikeinterface/core/testing.py index 6ad7458c37..d19b2f14f3 100644 --- a/src/spikeinterface/core/testing.py +++ b/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, diff --git a/src/spikeinterface/preprocessing/__init__.py b/src/spikeinterface/preprocessing/__init__.py index 957b5a97cc..0ca4567f59 100644 --- a/src/spikeinterface/preprocessing/__init__.py +++ b/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, diff --git a/src/spikeinterface/preprocessing/filter_design.py b/src/spikeinterface/preprocessing/filter_design.py new file mode 100644 index 0000000000..2a2a1a6766 --- /dev/null +++ b/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 diff --git a/src/spikeinterface/preprocessing/fir_filter.py b/src/spikeinterface/preprocessing/fir_filter.py new file mode 100644 index 0000000000..b9c53c54c5 --- /dev/null +++ b/src/spikeinterface/preprocessing/fir_filter.py @@ -0,0 +1,218 @@ +from __future__ import annotations + +import math +import threading + +import numpy as np + +from spikeinterface.core.core_tools import define_function_handling_dict_from_class + +from .basepreprocessor import BasePreprocessor, BasePreprocessorSegment +from .filter import fix_dtype + + +def _detect_l2_bytes(): + """Per-core L2 size from Linux sysfs, else 1 MB, which fits in any current core's L2.""" + for index in range(5): + cache = f"/sys/devices/system/cpu/cpu0/cache/index{index}" + try: + with open(f"{cache}/level") as f: + level = int(f.read()) + with open(f"{cache}/type") as f: + cache_type = f.read().strip() + if level != 2 or cache_type not in ("Unified", "Data"): + continue + with open(f"{cache}/size") as f: + size = f.read().strip() + except (OSError, ValueError): + continue + scale = {"K": 1024, "M": 1024**2}.get(size[-1:], 1) + return int(size.rstrip("KM")) * scale + return 1_000_000 + + +_L2_BYTES = _detect_l2_bytes() + + +def _overlap_save_fft_length(signal_length, num_taps): + """FFT length minimising FFT work per output sample, as ``scipy.signal.oaconvolve`` chooses it. + + Unlike ``oaconvolve``, the length is rounded up to one that is fast for a real FFT: lengths with a factor of 7 + or 11 are fast for a complex FFT but up to 1.5 times slower for ``rfft``. + """ + import scipy.fft + from scipy.special import lambertw + + whole = int(scipy.fft.next_fast_len(signal_length + num_taps - 1, real=True)) + if num_taps == 1 or num_taps >= signal_length / 2: + return whole + overlap = num_taps - 1 + optimal = -overlap * lambertw(-1 / (2 * math.e * overlap), k=-1).real + return min(whole, int(scipy.fft.next_fast_len(math.ceil(optimal), real=True))) + + +def _channels_per_block(num_channels, fft_length, bytes_per_sample, pool_size=1, headroom=0.6): + """ + Channels per overlap-save block: as many as keep the block's buffers in L2, but no more than gives each of + `pool_size` workers a block. + """ + fit_l2 = max(1, int(headroom * _L2_BYTES / (bytes_per_sample * fft_length))) + return min(fit_l2, max(1, math.ceil(num_channels / pool_size))) + + +class FIRFilterRecording(BasePreprocessor): + """ + Applies a FIR filter to each channel, by overlap-save convolution. + + Output sample ``n`` is ``sum_j taps[j] * x[n + center - j]``, so with the default `center` a symmetric + kernel is zero-phase. Every output sample depends only on the ``len(taps) - 1`` input samples around it, + which the filter reads from the parent, so the output does not depend on how it is chunked. + + Parameters + ---------- + recording : Recording + The recording to filter. + taps : array-like + Kernel of shape ``(num_taps,)``, applied to every channel, or ``(num_kernels, num_taps)``, with each + channel's row given by `channel_groups`. Designed with, for example, ``design_fir_bandpass`` or + ``scipy.signal.firwin``. + channel_groups : array-like | None, default: None + For a 2D `taps`, the row of `taps` for each channel. None requires one row per channel. + center : int | None, default: None + Index of the tap aligned with the output sample. None means ``(num_taps - 1) // 2``, zero-phase for a + symmetric kernel; 0 makes the filter causal. + dtype : dtype | None, default: None + The dtype of the output. None uses the parent's. Computation is in float32, and integer output is + rounded. + add_reflect_padding : bool, default: False + If True, the input is reflected beyond the recording's edges, otherwise it is zero there. + + Returns + ------- + filtered_recording : FIRFilterRecording + The filtered recording. + """ + + def __init__(self, recording, taps, channel_groups=None, center=None, dtype=None, add_reflect_padding=False): + taps = np.asarray(taps, dtype=np.float64) + if taps.ndim not in (1, 2) or taps.shape[-1] == 0: + raise ValueError(f"taps must have shape (num_taps,) or (num_kernels, num_taps), got {taps.shape}") + kernels = np.atleast_2d(taps) + num_kernels, num_taps = kernels.shape + num_channels = recording.get_num_channels() + if channel_groups is None: + if num_kernels not in (1, num_channels): + raise ValueError(f"{num_kernels} kernels for {num_channels} channels need channel_groups") + groups = np.zeros(num_channels, dtype=np.intp) if num_kernels == 1 else np.arange(num_channels) + else: + groups = np.asarray(channel_groups, dtype=np.intp) + if groups.shape != (num_channels,) or groups.min() < 0 or groups.max() >= num_kernels: + raise ValueError(f"channel_groups must give a row of taps, 0 to {num_kernels - 1}, for each channel") + center = (num_taps - 1) // 2 if center is None else int(center) + if not 0 <= center < num_taps: + raise ValueError(f"center must be in [0, {num_taps}), got {center}") + 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(): + # A constant offset passes through scaled by each kernel's DC gain. + self.set_channel_offsets(self.get_channel_offsets() * kernels.sum(axis=1)[groups]) + + for parent_segment in recording.segments: + self.add_recording_segment( + FIRFilterRecordingSegment(parent_segment, kernels, groups, center, dtype, add_reflect_padding) + ) + + self._kwargs = dict( + recording=recording, + taps=taps.tolist(), + channel_groups=None if channel_groups is None else groups.tolist(), + center=center, + dtype=dtype.str, + add_reflect_padding=add_reflect_padding, + ) + + +class FIRFilterRecordingSegment(BasePreprocessorSegment): + def __init__(self, parent_recording_segment, kernels, channel_groups, center, dtype, add_reflect_padding): + BasePreprocessorSegment.__init__(self, parent_recording_segment) + self.kernels = kernels + self.num_taps = kernels.shape[1] + self.channel_groups = channel_groups + self.center = center + self.dtype = dtype + self.add_reflect_padding = add_reflect_padding + # Kernel spectra depend on the FFT length, which depends on the chunk length. + self._kernel_spectra = {} + self._kernel_spectra_lock = threading.Lock() + + def _get_kernel_spectra(self, fft_length): + import scipy.fft + + with self._kernel_spectra_lock: + spectra = self._kernel_spectra.get(fft_length) + if spectra is None: + spectra = scipy.fft.rfft(self.kernels.T, n=fft_length, axis=0).astype(np.complex64) + self._kernel_spectra[fft_length] = spectra + return spectra + + def _get_input(self, start_frame, end_frame, channel_indices): + """Parent traces for every input sample the kernel reaches, padded beyond the recording's edges.""" + left = self.num_taps - 1 - self.center + right = self.center + fetch_start = max(0, start_frame - left) + fetch_end = min(self.parent_recording_segment.get_num_samples(), end_frame + right) + traces = self.parent_recording_segment.get_traces(fetch_start, fetch_end, channel_indices) + pad = ((fetch_start - (start_frame - left), end_frame + right - fetch_end), (0, 0)) + if pad[0] != (0, 0): + traces = np.pad(traces, pad, mode="reflect" if self.add_reflect_padding else "constant") + return traces + + def get_traces(self, start_frame, end_frame, channel_indices): + import scipy.fft + + if start_frame is None: + start_frame = 0 + if end_frame is None: + end_frame = self.get_num_samples() + if channel_indices is None: + channel_indices = slice(None) + traces = self._get_input(start_frame, end_frame, channel_indices) + groups = self.channel_groups[channel_indices] + num_samples, num_channels = end_frame - start_frame, traces.shape[1] + num_taps = self.num_taps + + fft_length = _overlap_save_fft_length(num_samples + num_taps - 1, num_taps) + step = fft_length - num_taps + 1 + spectra = self._get_kernel_spectra(fft_length) + one_kernel = spectra.shape[1] == 1 + # The float32 block and its complex64 spectrum take 8 bytes per sample, and per-channel kernel spectra 4 more. + slab = _channels_per_block(num_channels, fft_length, bytes_per_sample=8 if one_kernel else 12) + + # Each block is rounded and cast as it is written, so no full float32 copy of the chunk is held. + out = np.empty((num_samples, num_channels), dtype=self.dtype) + round_to_integer = np.issubdtype(self.dtype, np.integer) + buffer = np.empty((fft_length, slab), dtype=np.float32) + for c0 in range(0, num_channels, slab): + c1 = min(num_channels, c0 + slab) + tile = buffer[:, : c1 - c0] + tile_spectra = spectra if one_kernel else spectra[:, groups[c0:c1]] + for t0 in range(0, num_samples, step): + num_out = min(step, num_samples - t0) + num_in = num_out + num_taps - 1 + tile[:num_in] = traces[t0 : t0 + num_in, c0:c1] + tile[num_in:] = 0 + spectrum = scipy.fft.rfft(tile, axis=0) + spectrum *= tile_spectra + # The first num_taps - 1 samples of each circular convolution wrap around and are discarded. + result = scipy.fft.irfft(spectrum, n=fft_length, axis=0, overwrite_x=True) + valid = result[num_taps - 1 : num_taps - 1 + num_out] + if round_to_integer: + np.rint(valid, out=valid) + out[t0 : t0 + num_out, c0:c1] = valid + + return out + + +fir_filter = define_function_handling_dict_from_class(source_class=FIRFilterRecording, name="fir_filter") diff --git a/src/spikeinterface/preprocessing/preprocessing_classes.py b/src/spikeinterface/preprocessing/preprocessing_classes.py index 406e239380..d1a7aed7a6 100644 --- a/src/spikeinterface/preprocessing/preprocessing_classes.py +++ b/src/spikeinterface/preprocessing/preprocessing_classes.py @@ -13,6 +13,7 @@ causal_filter, ) from .filter_gaussian import GaussianFilterRecording, gaussian_filter +from .fir_filter import FIRFilterRecording, fir_filter from .normalize_scale import ( NormalizeByQuantileRecording, normalize_by_quantile, @@ -60,6 +61,7 @@ HighpassFilterRecording: highpass_filter, NotchFilterRecording: notch_filter, GaussianFilterRecording: gaussian_filter, + FIRFilterRecording: fir_filter, # gain offset stuff NormalizeByQuantileRecording: normalize_by_quantile, ScaleRecording: scale, diff --git a/src/spikeinterface/preprocessing/tests/test_fir_filter.py b/src/spikeinterface/preprocessing/tests/test_fir_filter.py new file mode 100644 index 0000000000..a24b58d537 --- /dev/null +++ b/src/spikeinterface/preprocessing/tests/test_fir_filter.py @@ -0,0 +1,207 @@ +import numpy as np +import pytest +import scipy.signal + +from spikeinterface import NumpyRecording, load +from spikeinterface.core.testing import check_chunk_invariance +from spikeinterface.preprocessing import design_fir_bandpass, design_fir_highpass, fir_filter + +FS = 30_000.0 + + +@pytest.fixture(scope="module") +def recording(): + rng = np.random.default_rng(0) + traces = (rng.standard_normal((20_000, 6)) * 50).astype(np.float32) + return NumpyRecording(traces, FS) + + +def expected_output(traces, taps, center): + """``sum_j taps[j] * x[n + center - j]`` with zeros beyond the edges, by direct convolution.""" + full = np.stack([np.convolve(traces[:, c].astype(np.float64), taps) for c in range(traces.shape[1])], axis=1) + return full[center : center + traces.shape[0]] + + +def bandpass_taps(): + return design_fir_bandpass(300.0, 6000.0, FS) + + +@pytest.mark.parametrize("center", [None, 0, 7]) +def test_fir_filter_matches_direct_convolution(recording, center): + taps = np.random.default_rng(1).standard_normal(31) + traces = recording.get_traces() + out = fir_filter(recording, taps, center=center).get_traces() + expected = expected_output(traces, taps, 15 if center is None else center) + np.testing.assert_allclose(out, expected, atol=1e-5 * np.abs(expected).max()) + + +def test_fir_filter_causal(recording): + # With center=0 an impulse at sample n moves only samples n and later. + traces = np.zeros((1000, 1), dtype=np.float32) + traces[500] = 1.0 + taps = np.arange(1.0, 6.0) + out = fir_filter(NumpyRecording(traces, FS), taps, center=0).get_traces()[:, 0] + np.testing.assert_allclose(out[:500], 0.0, atol=1e-6) + np.testing.assert_allclose(out[500:505], taps, rtol=1e-6) + + +def test_fir_filter_chunk_invariance(recording): + rec = fir_filter(recording, bandpass_taps()) + scale = np.abs(rec.get_traces()).max() + for chunk_size in (997, 3000, 7001): + check_chunk_invariance(rec, chunk_size, atol=1e-5 * scale) + + +def test_fir_filter_per_channel_kernels(recording): + rng = np.random.default_rng(2) + kernels = rng.standard_normal((3, 21)) + groups = np.array([0, 1, 2, 0, 1, 2]) + traces = recording.get_traces() + out = fir_filter(recording, kernels, channel_groups=groups).get_traces() + for channel, group in enumerate(groups): + expected = expected_output(traces[:, [channel]], kernels[group], 10) + np.testing.assert_allclose(out[:, [channel]], expected, atol=1e-5 * np.abs(expected).max()) + # One row per channel needs no channel_groups. + one_per_channel = rng.standard_normal((6, 21)) + out = fir_filter(recording, one_per_channel).get_traces() + expected = expected_output(traces[:, [4]], one_per_channel[4], 10) + np.testing.assert_allclose(out[:, [4]], expected, atol=1e-5 * np.abs(expected).max()) + + +def test_fir_filter_channel_blocks(recording, monkeypatch): + import sys + + # Blocks of 4 split the 6 channels unevenly and across kernel groups. + monkeypatch.setattr( + sys.modules["spikeinterface.preprocessing.fir_filter"], "_channels_per_block", lambda *a, **k: 4 + ) + kernels = np.random.default_rng(3).standard_normal((2, 21)) + groups = np.array([0, 1, 1, 0, 1, 0]) + traces = recording.get_traces() + out = fir_filter(recording, kernels, channel_groups=groups).get_traces() + expected = np.concatenate([expected_output(traces[:, [c]], kernels[g], 10) for c, g in enumerate(groups)], axis=1) + np.testing.assert_allclose(out, expected, atol=1e-5 * np.abs(expected).max()) + + +def test_fir_filter_channel_subset(recording): + rng = np.random.default_rng(3) + kernels = rng.standard_normal((6, 41)) + rec = fir_filter(recording, kernels) + whole = rec.get_traces() + ids = recording.channel_ids[[4, 1]] + subset = rec.get_traces(start_frame=1234, end_frame=5678, channel_ids=ids) + np.testing.assert_allclose(subset, whole[1234:5678, [4, 1]], atol=1e-5 * np.abs(whole).max()) + + +def test_fir_filter_integer_output_rounds(recording): + int_rec = NumpyRecording(np.rint(recording.get_traces()).astype(np.int16), FS) + taps = bandpass_taps() + out = fir_filter(int_rec, taps).get_traces() + assert out.dtype == np.int16 + as_float = fir_filter(int_rec, taps, dtype="float32").get_traces() + np.testing.assert_array_equal(out, np.rint(as_float).astype(np.int16)) + + +def test_fir_filter_reflect_padding(): + # A constant signal reflected at the edges stays constant through a unit-gain lowpass. + traces = np.full((2000, 1), 10.0, dtype=np.float32) + taps = scipy.signal.firwin(101, 1000.0, fs=FS) + rec = NumpyRecording(traces, FS) + np.testing.assert_allclose(fir_filter(rec, taps, add_reflect_padding=True).get_traces(), traces, rtol=1e-4) + assert abs(fir_filter(rec, taps).get_traces()[0, 0] - 10.0) > 1.0 + + +def test_fir_filter_offsets_scale_by_dc_gain(recording): + rec = NumpyRecording(recording.get_traces(), FS) + rec.set_channel_gains(2.0) + rec.set_channel_offsets(5.0) + taps = bandpass_taps() + np.testing.assert_allclose(fir_filter(rec, taps).get_channel_offsets(), 5.0 * taps.sum()) + lowpass = scipy.signal.firwin(101, 1000.0, fs=FS) + np.testing.assert_allclose(fir_filter(rec, lowpass).get_channel_offsets(), 5.0 * lowpass.sum()) + + +def test_fir_filter_save_load(tmp_path): + from spikeinterface.core import generate_recording + + parent = generate_recording(num_channels=4, durations=[1.0], sampling_frequency=FS, seed=0) + rec = fir_filter(parent, np.random.default_rng(4).standard_normal((2, 15)), channel_groups=[0, 1, 0, 1]) + rec.dump_to_json(tmp_path / "rec.json") + np.testing.assert_array_equal(load(tmp_path / "rec.json").get_traces(), rec.get_traces()) + saved = rec.save(folder=tmp_path / "saved") + np.testing.assert_array_equal(saved.get_traces(), rec.get_traces()) + + +def test_fir_filter_invalid_arguments(recording): + with pytest.raises(ValueError, match="channel_groups"): + fir_filter(recording, np.ones((3, 5))) + with pytest.raises(ValueError, match="channel_groups"): + fir_filter(recording, np.ones((3, 5)), channel_groups=[0, 1, 2, 3, 0, 1]) + with pytest.raises(ValueError, match="center"): + fir_filter(recording, np.ones(5), center=5) + + +def response_db(taps, freqs): + return 20 * np.log10(np.abs(scipy.signal.freqz(taps, worN=np.asarray(freqs, dtype=float), fs=FS)[1])) + + +@pytest.mark.parametrize( + "freq_min, freq_max, transition_low_hz, transition_high_hz, stopband_db", + [(300.0, 6000.0, 175.0, 2000.0, 45.0), (300.0, 6000.0, 400.0, 2000.0, 60.0), (500.0, 3000.0, 600.0, 300.0, 40.0)], +) +def test_design_fir_bandpass_meets_spec(freq_min, freq_max, transition_low_hz, transition_high_hz, stopband_db): + taps = design_fir_bandpass(freq_min, freq_max, FS, transition_low_hz, transition_high_hz, stopband_db) + assert taps.size % 2 == 1 + np.testing.assert_allclose(taps, taps[::-1], atol=1e-15) + np.testing.assert_allclose(response_db(taps, [freq_min, freq_max]), -6.02, atol=0.1) + stop_low = np.linspace(0.0, freq_min - transition_low_hz / 2, 2000) + stop_high = np.linspace(freq_max + transition_high_hz / 2, FS / 2, 2000) + assert response_db(taps, np.r_[stop_low, stop_high]).max() <= -stopband_db + passband = np.linspace(freq_min + transition_low_hz / 2, freq_max - transition_high_hz / 2, 2000) + delta = 10 ** (-stopband_db / 20) + gain = 10 ** (response_db(taps, passband) / 20) + assert np.abs(gain - 1).max() <= 2 * delta + # Each edge keeps its own width: a quarter of the way into the wider transition, one window sized for the + # narrower edge would already be in its stopband. + width = max(transition_low_hz, transition_high_hz) + inside = freq_min - width / 4 if transition_low_hz >= transition_high_hz else freq_max + width / 4 + assert response_db(taps, [inside])[0] > -30 + + +@pytest.mark.parametrize("freq_min, transition_low_hz, stopband_db", [(300.0, 175.0, 45.0), (300.0, 400.0, 60.0)]) +def test_design_fir_highpass_meets_spec(freq_min, transition_low_hz, stopband_db): + taps = design_fir_highpass(freq_min, FS, transition_low_hz, stopband_db) + assert taps.size % 2 == 1 + np.testing.assert_allclose(taps, taps[::-1], atol=1e-15) + np.testing.assert_allclose(response_db(taps, [freq_min]), -6.02, atol=0.1) + assert response_db(taps, np.linspace(0.0, freq_min - transition_low_hz / 2, 2000)).max() <= -stopband_db + passband = np.linspace(freq_min + transition_low_hz / 2, FS / 2, 4000) + gain = 10 ** (response_db(taps, passband) / 20) + assert np.abs(gain - 1).max() <= 2 * 10 ** (-stopband_db / 20) + + +def test_design_fir_highpass_invalid_arguments(): + with pytest.raises(ValueError, match="transition band"): + design_fir_highpass(80.0, FS) + with pytest.raises(ValueError, match="positive"): + design_fir_highpass(300.0, FS, stopband_db=0.0) + + +def test_fir_filter_bandpass_attenuates_out_of_band(): + t = np.arange(60_000) / FS + tones = {50.0: 1.0, 1000.0: 1.0, 12_000.0: 1.0} + traces = sum(np.sin(2 * np.pi * f * t) for f in tones)[:, np.newaxis].astype(np.float32) + out = fir_filter(NumpyRecording(traces, FS), bandpass_taps()).get_traces()[10_000:-10_000, 0] + t = t[10_000:-10_000] + np.testing.assert_allclose(out, np.sin(2 * np.pi * 1000.0 * t), atol=4 * 10 ** (-45 / 20)) + + +def test_design_fir_bandpass_invalid_arguments(): + with pytest.raises(ValueError, match="transition bands"): + design_fir_bandpass(80.0, 6000.0, FS) + with pytest.raises(ValueError, match="transition bands"): + design_fir_bandpass(300.0, 14_500.0, FS) + with pytest.raises(ValueError, match="transition bands"): + design_fir_bandpass(6000.0, 300.0, FS) + with pytest.raises(ValueError, match="positive"): + design_fir_bandpass(300.0, 6000.0, FS, stopband_db=0.0) diff --git a/src/spikeinterface/preprocessing/tests/test_pipeline.py b/src/spikeinterface/preprocessing/tests/test_pipeline.py index 7b53b0ab92..52c7b92215 100644 --- a/src/spikeinterface/preprocessing/tests/test_pipeline.py +++ b/src/spikeinterface/preprocessing/tests/test_pipeline.py @@ -83,6 +83,9 @@ def test_pipeline_equiv_to_step(): elif pp_name == "remove_artifacts": pp_dict[pp_name] = {"list_triggers": [0.5]} pp_rec_from_class = pp_class(rec, list_triggers=[0.5]) + elif pp_name == "fir_filter": + pp_dict[pp_name] = {"taps": [0.25, 0.5, 0.25]} + pp_rec_from_class = pp_class(rec, taps=[0.25, 0.5, 0.25]) elif pp_name == "zero_channel_pad": pp_dict[pp_name] = {"num_channels": 2} pp_rec_from_class = pp_class(rec, num_channels=2) From f795fbf1d414fab152c192188c1f716460a3de67 Mon Sep 17 00:00:00 2001 From: Galen Lynch Date: Mon, 12 Oct 2026 15:08:10 -0700 Subject: [PATCH 2/3] test: import scipy inside the fir_filter tests that use it The core test jobs run without scipy and collect every test module, so a module-level scipy import failed collection there. --- src/spikeinterface/preprocessing/tests/test_fir_filter.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/spikeinterface/preprocessing/tests/test_fir_filter.py b/src/spikeinterface/preprocessing/tests/test_fir_filter.py index a24b58d537..3b9cfe7f78 100644 --- a/src/spikeinterface/preprocessing/tests/test_fir_filter.py +++ b/src/spikeinterface/preprocessing/tests/test_fir_filter.py @@ -1,6 +1,5 @@ import numpy as np import pytest -import scipy.signal from spikeinterface import NumpyRecording, load from spikeinterface.core.testing import check_chunk_invariance @@ -103,6 +102,8 @@ def test_fir_filter_integer_output_rounds(recording): def test_fir_filter_reflect_padding(): + import scipy.signal + # A constant signal reflected at the edges stays constant through a unit-gain lowpass. traces = np.full((2000, 1), 10.0, dtype=np.float32) taps = scipy.signal.firwin(101, 1000.0, fs=FS) @@ -112,6 +113,8 @@ def test_fir_filter_reflect_padding(): def test_fir_filter_offsets_scale_by_dc_gain(recording): + import scipy.signal + rec = NumpyRecording(recording.get_traces(), FS) rec.set_channel_gains(2.0) rec.set_channel_offsets(5.0) @@ -142,6 +145,8 @@ def test_fir_filter_invalid_arguments(recording): def response_db(taps, freqs): + import scipy.signal + return 20 * np.log10(np.abs(scipy.signal.freqz(taps, worN=np.asarray(freqs, dtype=float), fs=FS)[1])) From 578ff12e1051e1a9622579c21062e9ce77956ef0 Mon Sep 17 00:00:00 2001 From: Galen Lynch Date: Mon, 12 Oct 2026 16:58:14 -0700 Subject: [PATCH 3/3] docs: drop the measured slowdown from the rfft length docstring --- src/spikeinterface/preprocessing/fir_filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spikeinterface/preprocessing/fir_filter.py b/src/spikeinterface/preprocessing/fir_filter.py index b9c53c54c5..9e7cd2ab36 100644 --- a/src/spikeinterface/preprocessing/fir_filter.py +++ b/src/spikeinterface/preprocessing/fir_filter.py @@ -38,7 +38,7 @@ def _overlap_save_fft_length(signal_length, num_taps): """FFT length minimising FFT work per output sample, as ``scipy.signal.oaconvolve`` chooses it. Unlike ``oaconvolve``, the length is rounded up to one that is fast for a real FFT: lengths with a factor of 7 - or 11 are fast for a complex FFT but up to 1.5 times slower for ``rfft``. + or 11 are fast for a complex FFT but slower for ``rfft``. """ import scipy.fft from scipy.special import lambertw