From 5455c2ea5830a5d5a6e5fab61274863e7b08f7c3 Mon Sep 17 00:00:00 2001
From: Graham Findlay <graham.b.findlay@gmail.com>
Date: Wed, 17 Jun 2026 09:39:26 -0500
Subject: [PATCH 01/11] Add antialiased decimation and allow non-integer
 resampling rates

Adds antialised decimation to both `DecimateRecording` and
`ResampleRecording`. Opt-in for `DecimateRecording` with the `antialias`
parameter. For `ResampleRecording`, antialiasing is always applied.
---
 doc/api.rst                                   |   1 +
 .../preprocessing/_decimation_tools.py        | 175 ++++++++++++++++++
 src/spikeinterface/preprocessing/decimate.py  |  99 ++++++++--
 src/spikeinterface/preprocessing/resample.py  | 102 ++++++----
 .../preprocessing/tests/test_decimate.py      | 122 +++++++++++-
 5 files changed, 445 insertions(+), 54 deletions(-)
 create mode 100644 src/spikeinterface/preprocessing/_decimation_tools.py

diff --git a/doc/api.rst b/doc/api.rst
index ce850d1291..fd645e579c 100755
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -216,6 +216,7 @@ spikeinterface.preprocessing
     .. autofunction:: get_motion_parameters_preset
     .. autofunction:: load_motion_info
     .. autofunction:: save_motion_info
+    .. autofunction:: decimate
     .. autofunction:: depth_order
     .. autofunction:: detect_bad_channels
     .. autofunction:: detect_and_interpolate_bad_channels
diff --git a/src/spikeinterface/preprocessing/_decimation_tools.py b/src/spikeinterface/preprocessing/_decimation_tools.py
new file mode 100644
index 0000000000..fc5147c55e
--- /dev/null
+++ b/src/spikeinterface/preprocessing/_decimation_tools.py
@@ -0,0 +1,175 @@
+"""
+Helpers for splitting a (potentially large) integer decimation factor into several balanced
+sub-factors, so that anti-aliased decimation can be applied as multiple stable scipy.signal.decimate
+passes. Shared by DecimateRecording and ResampleRecording.
+"""
+
+import math
+import warnings
+
+import numpy as np
+
+from spikeinterface.core import get_chunk_with_margin
+
+# scipy.signal.decimate uses an order-8 Chebyshev type I IIR filter by default, and its
+# documentation recommends decimating in several balanced steps rather than a single step
+# for downsampling factors larger than this value.
+_MAX_SINGLE_PASS_DECIMATION = 13
+
+
+def _prime_factors(n):
+    """
+    Return the prime factors of a positive integer `n` (ascending, with multiplicity).
+
+    Examples
+    --------
+    >>> _prime_factors(60)
+    [2, 2, 3, 5]
+    >>> _prime_factors(17)
+    [17]
+    """
+    factors = []
+    divisor = 2
+    while divisor * divisor <= n:
+        while n % divisor == 0:
+            factors.append(divisor)
+            n //= divisor
+        divisor += 1
+    if n > 1:
+        factors.append(n)
+    return factors
+
+
+def _greedy_pack(primes_desc, num_bins):
+    """
+    Greedily pack `primes_desc` (largest first) into `num_bins` bins, keeping each bin's
+    product <= `_MAX_SINGLE_PASS_DECIMATION` and the bins as balanced as possible.
+
+    Returns the list of bin products, or None if some prime cannot be placed (i.e. `num_bins`
+    is too small to keep every bin <= the single-pass limit).
+
+    Examples
+    --------
+    Pack the prime factors of 48 into two balanced bins (6 and 8):
+
+    >>> _greedy_pack([3, 2, 2, 2, 2], 2)
+    [6, 8]
+
+    Two bins cannot hold 2 ** 7 = 128 without a bin exceeding the single-pass limit of 13:
+
+    >>> _greedy_pack([2, 2, 2, 2, 2, 2, 2], 2) is None
+    True
+    """
+    bins = [1] * num_bins
+    for prime in primes_desc:
+        fitting = [i for i in range(num_bins) if bins[i] * prime <= _MAX_SINGLE_PASS_DECIMATION]
+        if not fitting:
+            return None
+        # Place into the smallest fitting bin (ties broken by index, for determinism).
+        target = min(fitting, key=lambda i: (bins[i], i))
+        bins[target] *= prime
+    return bins
+
+
+def get_balanced_decimation_factors(decimation_factor):
+    """
+    Split `decimation_factor` into sub-factors, each <= 13, as balanced as possible (so their
+    products are close), for stable multi-pass anti-aliased decimation.
+
+    scipy recommends decimating in several balanced steps rather than one large step when the
+    factor exceeds 13 (e.g. 48 -> [8, 6] rather than [12, 4]). The product of the returned
+    factors always equals `decimation_factor`.
+
+    If `decimation_factor` has a prime factor greater than 13 (e.g. a large prime such as 17),
+    no valid split exists and `[decimation_factor]` is returned; it is the caller's
+    responsibility to handle this (e.g., warn that a single, potentially unstable, pass will
+    be used).
+    """
+    if decimation_factor <= _MAX_SINGLE_PASS_DECIMATION:
+        return [decimation_factor]
+
+    primes = _prime_factors(decimation_factor)
+    if max(primes) > _MAX_SINGLE_PASS_DECIMATION:
+        # If a prime factor > 13 cannot be split into sub-13 factors...
+        return [decimation_factor]
+
+    primes_desc = sorted(primes, reverse=True)
+    # Minimum number of passes so that, ideally, each pass decimates by <= 13.
+    num_passes = max(1, math.ceil(math.log(decimation_factor) / math.log(_MAX_SINGLE_PASS_DECIMATION)))
+    while num_passes <= len(primes_desc):
+        bins = _greedy_pack(primes_desc, num_passes)
+        if bins is not None:
+            return sorted(bins, reverse=True)
+        num_passes += 1
+    # Fallback: one prime per pass (always valid since every prime is <= 13).
+    return primes_desc
+
+
+def get_antialiased_decimated_traces(
+    parent_segment,
+    start_frame,
+    end_frame,
+    channel_indices,
+    decimation_factor,
+    decimation_factors,
+    margin,
+    dtype,
+    decimation_offset=0,
+):
+    """
+    Fetch a margined chunk from `parent_segment` and decimate it by `decimation_factor`, applied
+    as a cascade of the balanced `decimation_factors` passes of ``scipy.signal.decimate``.
+
+    The margin is rounded up to a multiple of the total `decimation_factor` so that
+    ``left_margin // decimation_factor`` is exact; combined with scipy's default
+    ``zero_phase=True`` (output sample i maps to filtered input sample i * factor), this keeps the
+    downsampled traces aligned across chunks (a chunked read matches a full read). Exactly
+    ``end_frame - start_frame`` decimated samples are returned.
+
+    Parameters
+    ----------
+    parent_segment : BaseRecordingSegment
+        The parent segment to read (full-rate) traces from.
+    start_frame, end_frame : int
+        Output (decimated) frame range to return.
+    channel_indices : slice | list | np.ndarray | None
+        Channels to read, forwarded to the parent segment.
+    decimation_factor : int
+        The total decimation factor (the product of `decimation_factors`).
+    decimation_factors : list[int]
+        The per-pass sub-factors (each <= 13), e.g. from `get_balanced_decimation_factors`.
+    margin : int
+        Margin in parent samples used to limit anti-aliasing filter edge effects. Rounded up
+        internally to a multiple of `decimation_factor`.
+    dtype : np.dtype | str
+        Output dtype. The decimation runs in float32 and the result is cast to `dtype`.
+    decimation_offset : int, default: 0
+        Index of the first parent frame, applied to the first output sample only.
+    """
+    from scipy import signal
+
+    q = decimation_factor
+    parent_start_frame = decimation_offset + start_frame * q
+    parent_end_frame = parent_start_frame + (end_frame - start_frame) * q
+    # Round the margin up to a multiple of q so that left_margin // q is exact.
+    margin = int(np.ceil(margin / q) * q)
+    parent_traces, left_margin, right_margin = get_chunk_with_margin(
+        parent_segment,
+        parent_start_frame,
+        parent_end_frame,
+        channel_indices,
+        margin,
+        add_reflect_padding=True,
+        dtype=np.float32,
+    )
+    decimated_traces = parent_traces
+    for sub_q in decimation_factors:
+        decimated_traces = signal.decimate(decimated_traces, q=sub_q, axis=0)
+    if np.any(np.isnan(decimated_traces)):
+        warnings.warn(
+            f"`scipy.signal.decimate` produced NaNs while decimating by {q}. "
+            f"Consider a different decimation factor."
+        )
+    start_drop = left_margin // q
+    n_out = end_frame - start_frame
+    return decimated_traces[start_drop : start_drop + n_out].astype(dtype)
diff --git a/src/spikeinterface/preprocessing/decimate.py b/src/spikeinterface/preprocessing/decimate.py
index 716cc46f70..348462ed6f 100644
--- a/src/spikeinterface/preprocessing/decimate.py
+++ b/src/spikeinterface/preprocessing/decimate.py
@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 from spikeinterface.core.core_tools import (
     define_function_handling_dict_from_class,
@@ -5,16 +7,23 @@
 
 from .basepreprocessor import BasePreprocessor
 from .filter import fix_dtype
+from ._decimation_tools import (
+    _MAX_SINGLE_PASS_DECIMATION,
+    get_balanced_decimation_factors,
+    get_antialiased_decimated_traces,
+)
 from spikeinterface.core import BaseRecordingSegment
 
 
 class DecimateRecording(BasePreprocessor):
     """
-    Decimate the recording extractor traces using array slicing
+    Decimate the recording extractor traces.
 
-    Important: This uses simple array slicing for decimation rather than eg scipy.decimate.
-    This might introduce aliasing, or skip across signal of interest.
-    Consider  spikeinterface.preprocessing.ResampleRecording for safe resampling.
+    By default this uses simple array slicing
+    (``<parent_traces>[<decimation_offset>::<decimation_factor>]``), which is fast but applies no
+    anti-aliasing filter and so might introduce aliasing, or skip across signal of interest. Set
+    `antialias=True` to low-pass filter before downsampling using ``scipy.signal.decimate`` (the
+    same anti-aliased decimation used by ``spikeinterface.preprocessing.ResampleRecording``).
 
     Parameters
     ----------
@@ -29,12 +38,25 @@ class DecimateRecording(BasePreprocessor):
         to ensure that the decimated recording has at least one frame. Consider combining DecimateRecording
         with FrameSliceRecording for fine control on the recording start and end frames.
         The same decimation offset is applied to all segments from the parent recording.
+    antialias : bool, default: False
+        If True, apply an anti-aliasing low-pass filter before downsampling, using
+        ``scipy.signal.decimate``. When `decimation_factor` exceeds 13, the decimation is
+        automatically performed in several balanced sub-13 passes (e.g. a factor of 48 is applied
+        as 8 then 6), as scipy recommends, to keep the IIR anti-aliasing filter stable. If False
+        (the default), traces are downsampled by plain array slicing with no filtering, and
+        `margin_ms` is ignored.
+    margin_ms : float, default: 100.0
+        Margin in ms used on each side of every chunk to limit edge effects of the anti-aliasing
+        filter. Only used when `antialias=True`. The margin is internally rounded up to a whole
+        number of output samples so the filtered, downsampled traces stay aligned across chunks.
+    dtype : dtype or None, default: None
+        The dtype of the returned traces. If None, the dtype of the parent recording is used.
 
     Returns
     -------
     decimate_recording: DecimateRecording
-        The decimated recording extractor object. The full traces of the child recording segment
-        correspond to the traces of the parent segment as follows:
+        The decimated recording extractor object. With `antialias=False` the full traces of the
+        child recording segment correspond to the traces of the parent segment as follows:
             ```<decimated_traces> = <parent_traces>[<decimation_offset>::<decimation_factor>]```
 
     """
@@ -44,6 +66,9 @@ def __init__(
         recording,
         decimation_factor,
         decimation_offset=0,
+        antialias=False,
+        margin_ms=100.0,
+        dtype=None,
     ):
         # Original sampling frequency
         self._orig_samp_freq = recording.get_sampling_frequency()
@@ -63,7 +88,21 @@ def __init__(
         self._decimation_offset = decimation_offset
         decimated_sampling_frequency = self._orig_samp_freq / self._decimation_factor
 
-        BasePreprocessor.__init__(self, recording, sampling_frequency=decimated_sampling_frequency)
+        # fix_dtype doesn't always returns the str, make sure it does
+        dtype = fix_dtype(recording, dtype).str
+
+        antialias_factors = get_balanced_decimation_factors(decimation_factor)
+        if antialias and decimation_factor > _MAX_SINGLE_PASS_DECIMATION and antialias_factors == [decimation_factor]:
+            warnings.warn(
+                f"`decimation_factor`={decimation_factor} cannot be split into anti-aliasing passes of <= 13 "
+                f"(it has a prime factor > 13). A single `scipy.signal.decimate` pass will be used, which may be "
+                f"unstable. Consider a `decimation_factor` without large prime factors."
+            )
+
+        # Margin (in parent samples) to limit anti-aliasing filter edge effects.
+        margin = int(margin_ms * self._orig_samp_freq / 1000)
+
+        BasePreprocessor.__init__(self, recording, sampling_frequency=decimated_sampling_frequency, dtype=dtype)
 
         for parent_segment in recording.segments:
             self.add_recording_segment(
@@ -74,6 +113,9 @@ def __init__(
                     decimation_factor,
                     decimation_offset,
                     self._dtype,
+                    antialias,
+                    margin,
+                    antialias_factors,
                 )
             )
 
@@ -81,6 +123,9 @@ def __init__(
             recording=recording,
             decimation_factor=decimation_factor,
             decimation_offset=decimation_offset,
+            antialias=antialias,
+            margin_ms=margin_ms,
+            dtype=dtype,
         )
 
 
@@ -93,6 +138,9 @@ def __init__(
         decimation_factor,
         decimation_offset,
         dtype,
+        antialias=False,
+        margin=0,
+        antialias_factors=None,
     ):
         if parent_recording_segment._time_vector is not None:
             time_vector = parent_recording_segment._time_vector[decimation_offset::decimation_factor]
@@ -113,6 +161,9 @@ def __init__(
         self._decimation_factor = decimation_factor
         self._decimation_offset = decimation_offset
         self._dtype = dtype
+        self._antialias = antialias
+        self._margin = margin
+        self._antialias_factors = antialias_factors if antialias_factors is not None else [decimation_factor]
 
     def get_num_samples(self):
         parent_n_samp = self._parent_segment.get_num_samples()
@@ -120,18 +171,30 @@ def get_num_samples(self):
         return int(np.ceil((parent_n_samp - self._decimation_offset) / self._decimation_factor))
 
     def get_traces(self, start_frame, end_frame, channel_indices):
-        # Account for offset and end when querying parent traces
-        parent_start_frame = self._decimation_offset + start_frame * self._decimation_factor
-        parent_end_frame = parent_start_frame + (end_frame - start_frame) * self._decimation_factor
-
-        # And now we can decimate without offsetting
-        return self._parent_segment.get_traces(
-            parent_start_frame,
-            parent_end_frame,
+        if not self._antialias:
+            # Simple array slicing, no anti-aliasing filter.
+            parent_start_frame = self._decimation_offset + start_frame * self._decimation_factor
+            parent_end_frame = parent_start_frame + (end_frame - start_frame) * self._decimation_factor
+            return self._parent_segment.get_traces(
+                parent_start_frame,
+                parent_end_frame,
+                channel_indices,
+            )[
+                :: self._decimation_factor
+            ].astype(self._dtype)
+
+        # Anti-aliased decimation as a cascade of balanced scipy.signal.decimate passes.
+        return get_antialiased_decimated_traces(
+            self._parent_segment,
+            start_frame,
+            end_frame,
             channel_indices,
-        )[
-            :: self._decimation_factor
-        ].astype(self._dtype)
+            self._decimation_factor,
+            self._antialias_factors,
+            self._margin,
+            self._dtype,
+            decimation_offset=self._decimation_offset,
+        )
 
 
 decimate = define_function_handling_dict_from_class(source_class=DecimateRecording, name="decimate")
diff --git a/src/spikeinterface/preprocessing/resample.py b/src/spikeinterface/preprocessing/resample.py
index a801c45eff..d6bf527287 100644
--- a/src/spikeinterface/preprocessing/resample.py
+++ b/src/spikeinterface/preprocessing/resample.py
@@ -8,6 +8,11 @@
 
 from .basepreprocessor import BasePreprocessor
 from .filter import fix_dtype
+from ._decimation_tools import (
+    _MAX_SINGLE_PASS_DECIMATION,
+    get_balanced_decimation_factors,
+    get_antialiased_decimated_traces,
+)
 from spikeinterface.core import get_chunk_with_margin, BaseRecordingSegment
 
 
@@ -15,17 +20,19 @@ class ResampleRecording(BasePreprocessor):
     """
     Resample the recording extractor traces.
 
-    If the original sampling rate is multiple of the resample_rate, it will use
-    the signal.decimate method from scipy. In other cases, it uses signal.resample. In the
-    later case, the resulting signal can have issues on the edges, mainly on the
-    rightmost.
+    If the parent sampling rate is an exact integer multiple of `resample_rate`, the
+    ``signal.decimate`` method from scipy is used (anti-aliased decimation). In other cases
+    ``signal.resample`` is used, in which case the resulting signal can have issues on the edges,
+    mainly on the rightmost. See Notes for a caveat on how the integer multiple is detected.
 
     Parameters
     ----------
     recording : Recording
         The recording extractor to be re-referenced
-    resample_rate : int
-        The resampling frequency
+    resample_rate : int | float
+        The resampling frequency. Integer ratios (parent_rate / resample_rate) use
+        ``scipy.signal.decimate``; non-integer ratios use ``scipy.signal.resample`` (FFT-based),
+        which can have edge effects, mainly on the rightmost samples.
     gap_tolerance_ms : float | None, default: None
         Maximum acceptable gap size in milliseconds for automatic segmentation.
 
@@ -60,6 +67,19 @@ class ResampleRecording(BasePreprocessor):
     resample_recording : ResampleRecording
         The resampled recording extractor object.
 
+    Notes
+    -----
+    The (anti-aliased) decimation path is selected by an exact check,
+    ``parent_rate % resample_rate == 0``. This only detects an integer downsampling factor when
+    both rates make that modulo exactly zero. If either the parent rate or `resample_rate` is a
+    non-integer float (i.e. ``float(int(x)) != float(x)``), a conceptually integer ratio can go
+    undetected and silently fall back to the FFT-based ``scipy.signal.resample`` path. For example,
+    decimating a 625 Hz recording by a factor of 6 means a target of 104.1666... Hz, and
+    ``625 % 104.1666... != 0``, so the integer-decimation path is not taken (whereas a factor of 5,
+    i.e. a 125 Hz target, is detected since ``625 % 125 == 0``). To force anti-aliased integer
+    decimation by a known factor regardless of the rates, use
+    ``spikeinterface.preprocessing.DecimateRecording`` with ``antialias=True``.
+
     """
 
     def __init__(
@@ -71,13 +91,24 @@ def __init__(
         dtype=None,
         skip_checks=False,
     ):
-        # Floating point resampling rates can lead to unexpected results, avoid actively
-        msg = "Non integer resampling rates can lead to unexpected results."
-        assert isinstance(resample_rate, (int, np.integer)), msg
-        # Original sampling frequency
         self._orig_samp_freq = recording.get_sampling_frequency()
         self._resample_rate = resample_rate
         self._sampling_frequency = resample_rate
+        # When the parent rate is an exact integer multiple of resample_rate, get_traces uses
+        # anti-aliased decimation. Large factors are split into several balanced sub-13 passes
+        # (see _decimation_tools); only an unsplittable factor (a prime > 13) falls back to a
+        # single, potentially unstable, pass and warns.
+        if self._orig_samp_freq % resample_rate == 0:
+            decimation_factor = int(self._orig_samp_freq / resample_rate)
+            decimation_factors = get_balanced_decimation_factors(decimation_factor)
+            if decimation_factors == [decimation_factor] and decimation_factor > _MAX_SINGLE_PASS_DECIMATION:
+                warnings.warn(
+                    f"Resampling by an integer factor of {decimation_factor} cannot be split into "
+                    f"anti-aliasing passes of <= 13 (it has a prime factor > 13); a single "
+                    f"`scipy.signal.decimate` pass will be used, which may be unstable."
+                )
+        else:
+            decimation_factors = None
         # fix_dtype not always returns the str, make sure it does
         dtype = fix_dtype(recording, dtype).str
         # Ensure that the requested resample rate is doable:
@@ -97,6 +128,7 @@ def __init__(
                     margin,
                     dtype,
                     gap_tolerance_ms,
+                    decimation_factors,
                 )
             )
 
@@ -119,6 +151,7 @@ def __init__(
         margin,
         dtype,
         gap_tolerance_ms=None,
+        decimation_factors=None,
     ):
         self._resample_rate = resample_rate
         self._parent_segment = parent_recording_segment
@@ -126,6 +159,9 @@ def __init__(
         self._margin = margin
         self._dtype = dtype
         self._has_gaps = False
+        # Per-pass integer decimation factors when the ratio is an exact integer, else None
+        # (non-integer ratio -> FFT-based scipy.signal.resample).
+        self._decimation_factors = decimation_factors
 
         # Compute time_vector or t_start, following the pattern from DecimateRecordingSegment.
         # Do not use BasePreprocessorSegment because we have to reset the sampling rate!
@@ -248,8 +284,23 @@ def get_traces(self, start_frame, end_frame, channel_indices):
         if self._has_gaps:
             return self._get_traces_gapped(start_frame, end_frame, channel_indices)
 
-        # Original code path: no gaps (or no time_vector)
-        # get parent traces with margin
+        if self._decimation_factors is not None:
+            # Integer ratio, no gaps: anti-aliased multi-pass decimation.
+            decimation_factor = int(self._parent_rate / self._resample_rate)
+            return get_antialiased_decimated_traces(
+                self._parent_segment,
+                start_frame,
+                end_frame,
+                channel_indices,
+                decimation_factor,
+                self._decimation_factors,
+                self._margin,
+                self._dtype,
+            )
+
+        # Non-integer ratio, no gaps: FFT-based resampling with proportional margins.
+        from scipy import signal
+
         parent_start_frame, parent_end_frame = [
             int((frame / self._resample_rate) * self._parent_rate) for frame in [start_frame, end_frame]
         ]
@@ -267,23 +318,9 @@ def get_traces(self, start_frame, end_frame, channel_indices):
             int((margin / self._parent_rate) * self._resample_rate) for margin in [left_margin, right_margin]
         ]
 
-        # get the size for the resampled traces in case of resample:
+        # get the size for the resampled traces
         num = int((end_frame + right_margin_rs) - (start_frame - left_margin_rs))
-
-        # Decimate can misbehave on some cases, while resample always looks nice enough.
-        # Check which method to use:
-        from scipy.signal import decimate, resample
-
-        if np.mod(self._parent_rate, self._resample_rate) == 0:
-            # Ratio between sampling frequencies
-            q = int(self._parent_rate / self._resample_rate)
-            # Decimate can have issues for some cases, returning NaNs
-            resampled_traces = decimate(parent_traces, q=q, axis=0)
-            # If that's the case, use signal.resample
-            if np.any(np.isnan(resampled_traces)):
-                resampled_traces = resample(parent_traces, num, axis=0)
-        else:
-            resampled_traces = resample(parent_traces, num, axis=0)
+        resampled_traces = signal.resample(parent_traces, num, axis=0)
 
         # now take care of the edges
         resampled_traces = resampled_traces[left_margin_rs : num - right_margin_rs]
@@ -368,10 +405,13 @@ def _get_traces_gapped(self, start_frame, end_frame, channel_indices):
             chunk_len = int(local_out_end - local_out_start)
             num = chunk_len + left_margin_rs + right_margin_rs
 
-            # Resample this section
+            # Resample this section. Integer ratios use anti-aliased multi-pass decimation
+            # (the balanced sub-13 passes computed in ResampleRecording.__init__), applied
+            # within the section only so filtering never crosses a gap.
             if is_integer_ratio:
-                q = int(self._parent_rate / self._resample_rate)
-                resampled = decimate(parent_traces, q=q, axis=0)
+                resampled = parent_traces
+                for sub_q in self._decimation_factors:
+                    resampled = decimate(resampled, q=sub_q, axis=0)
                 if np.any(np.isnan(resampled)):
                     resampled = resample(parent_traces, num, axis=0)
             else:
diff --git a/src/spikeinterface/preprocessing/tests/test_decimate.py b/src/spikeinterface/preprocessing/tests/test_decimate.py
index 93f70ea9bd..bdc880887d 100644
--- a/src/spikeinterface/preprocessing/tests/test_decimate.py
+++ b/src/spikeinterface/preprocessing/tests/test_decimate.py
@@ -2,8 +2,9 @@
 
 
 from spikeinterface import NumpyRecording
-from spikeinterface.core import generate_recording
-from spikeinterface.preprocessing.decimate import DecimateRecording
+from spikeinterface.core import generate_recording, load
+from spikeinterface.preprocessing.decimate import DecimateRecording, decimate, get_balanced_decimation_factors
+from spikeinterface.preprocessing.tests.test_resample import create_sinusoidal_traces
 import numpy as np
 
 
@@ -45,7 +46,8 @@ def test_decimate(num_segments, decimation_offset, decimation_factor):
         )
 
 
-def test_decimate_with_times():
+@pytest.mark.parametrize("antialias", [False, True])
+def test_decimate_with_times(antialias):
     rec = generate_recording(durations=[5, 10])
 
     # test with times
@@ -55,7 +57,7 @@ def test_decimate_with_times():
 
     decimation_factor = 2
     decimation_offset = 1
-    decimated_rec = DecimateRecording(rec, decimation_factor, decimation_offset=decimation_offset)
+    decimated_rec = DecimateRecording(rec, decimation_factor, decimation_offset=decimation_offset, antialias=antialias)
 
     for segment_index in range(rec.get_num_segments()):
         assert np.allclose(
@@ -68,7 +70,7 @@ def test_decimate_with_times():
     t_starts = [10, 20]
     for t_start, rec_segment in zip(t_starts, rec.segments):
         rec_segment._t_start = t_start
-    decimated_rec = DecimateRecording(rec, decimation_factor, decimation_offset=decimation_offset)
+    decimated_rec = DecimateRecording(rec, decimation_factor, decimation_offset=decimation_offset, antialias=antialias)
     for segment_index in range(rec.get_num_segments()):
         assert np.allclose(
             decimated_rec.get_times(segment_index),
@@ -76,5 +78,115 @@ def test_decimate_with_times():
         )
 
 
+@pytest.mark.parametrize(
+    "decimation_factor, expected",
+    [
+        (1, [1]),
+        (7, [7]),
+        (13, [13]),
+        (48, [8, 6]),
+        (50, [10, 5]),
+        (60, [10, 6]),
+        (100, [10, 10]),
+        (17, [17]),  # prime > 13: cannot be split
+        (23, [23]),  # prime > 13: cannot be split
+    ],
+)
+def test_balanced_decimation_factors(decimation_factor, expected):
+    factors = get_balanced_decimation_factors(decimation_factor)
+    assert factors == expected
+    # The product of the sub-factors always reconstructs the requested factor.
+    assert int(np.prod(factors)) == decimation_factor
+    # Every pass is <= 13 unless the factor is an unsplittable prime > 13.
+    if len(factors) > 1:
+        assert all(f <= 13 for f in factors)
+
+
+@pytest.mark.parametrize("decimation_factor", [6, 10, 48])
+def test_decimate_antialias_by_chunks(decimation_factor):
+    # Mirror test_resample_by_chunks: chunked reads must match a full read once the
+    # anti-aliasing margins are accounted for. Factor 48 exercises the internal multi-pass.
+    sampling_frequency = int(3e4)
+    duration = 30
+    traces, _ = create_sinusoidal_traces(sampling_frequency, duration, freqs_n=10, max_freq=1000, dtype=np.float32)
+    parent_rec = NumpyRecording(traces, sampling_frequency)
+    rms = np.sqrt(np.mean(parent_rec.get_traces() ** 2))
+    decimated_rate = sampling_frequency / decimation_factor
+
+    for margin_ms in [100, 1000]:
+        rec2 = DecimateRecording(parent_rec, decimation_factor, antialias=True, margin_ms=margin_ms)
+        chunk_size = int(decimated_rate * 2)  # ~2 seconds of the decimated signal
+        rec3 = rec2.save(format="memory", chunk_size=chunk_size, n_jobs=1, progress_bar=False)
+
+        traces2 = rec2.get_traces()
+        traces3 = rec3.get_traces()
+
+        # Drop the first and last chunk before comparing (as in test_resample_by_chunks).
+        sl = slice(chunk_size, -chunk_size)
+        error_mean = np.sqrt(np.mean((traces2[sl] - traces3[sl]) ** 2))
+        error_max = np.sqrt(np.max((traces2[sl] - traces3[sl]) ** 2))
+
+        assert error_mean / rms < 0.01
+        assert error_max / rms < 0.05
+
+
+@pytest.mark.parametrize("decimation_factor", [6, 10])
+@pytest.mark.parametrize("decimation_offset", [0, 1, 5])
+def test_decimate_antialias_with_offset(decimation_factor, decimation_offset):
+    sampling_frequency = 30000
+    # max_freq below every tested Nyquist, so anti-aliasing barely changes the signal.
+    traces, _ = create_sinusoidal_traces(sampling_frequency, duration=5, freqs_n=6, max_freq=500, dtype=np.float32)
+    parent_rec = NumpyRecording(traces, sampling_frequency)
+
+    dec_aa = DecimateRecording(
+        parent_rec, decimation_factor, decimation_offset=decimation_offset, antialias=True, dtype="float32"
+    )
+    dec_plain = DecimateRecording(
+        parent_rec, decimation_factor, decimation_offset=decimation_offset, antialias=False, dtype="float32"
+    )
+
+    # The anti-aliasing path returns the same number of samples as plain slicing.
+    parent_n = parent_rec.get_num_samples()
+    expected_n = int(np.ceil((parent_n - decimation_offset) / decimation_factor))
+    assert dec_aa.get_num_samples() == expected_n
+    assert dec_aa.get_num_samples() == dec_plain.get_num_samples()
+
+    # With only sub-Nyquist content, anti-aliased and plain-sliced traces stay aligned.
+    corr = np.corrcoef(dec_aa.get_traces().ravel(), dec_plain.get_traces().ravel())[0, 1]
+    assert corr > 0.95
+
+
+def test_decimate_antialias_multipass():
+    sampling_frequency = 30000
+    decimation_factor = 48
+    traces, _ = create_sinusoidal_traces(sampling_frequency, duration=10, freqs_n=8, max_freq=200, dtype=np.float32)
+    parent_rec = NumpyRecording(traces, sampling_frequency)
+
+    dec = decimate(parent_rec, decimation_factor, antialias=True)
+
+    # Multi-pass happens internally: a single DecimateRecording carries the full factor.
+    assert isinstance(dec, DecimateRecording)
+    assert dec._kwargs["decimation_factor"] == decimation_factor
+
+    segment = dec.segments[0]
+    assert int(np.prod(segment._antialias_factors)) == decimation_factor
+    assert all(f <= 13 for f in segment._antialias_factors)
+
+    parent_n = parent_rec.get_num_samples()
+    assert dec.get_num_samples() == int(np.ceil(parent_n / decimation_factor))
+
+    # Provenance round-trips and reproduces the traces.
+    dec_loaded = load(dec.to_dict())
+    np.testing.assert_allclose(dec_loaded.get_traces(), dec.get_traces())
+
+
+def test_decimate_antialias_large_prime_warns():
+    rec = generate_recording(durations=[2.0], num_channels=2)
+    with pytest.warns(UserWarning, match="prime factor > 13"):
+        dec = DecimateRecording(rec, 17, antialias=True)
+    # The unsplittable factor falls back to a single pass.
+    assert dec.segments[0]._antialias_factors == [17]
+
+
 if __name__ == "__main__":
     test_decimate()

From bfe7f537cb9c7e8f25d038208df28ea9b1d49727 Mon Sep 17 00:00:00 2001
From: Graham Findlay <graham.b.findlay@gmail.com>
Date: Fri, 7 Aug 2026 19:17:05 -0500
Subject: [PATCH 02/11] Round decimated traces before casting, per #4653

---
 src/spikeinterface/preprocessing/_decimation_tools.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/spikeinterface/preprocessing/_decimation_tools.py b/src/spikeinterface/preprocessing/_decimation_tools.py
index fc5147c55e..d6da9dd340 100644
--- a/src/spikeinterface/preprocessing/_decimation_tools.py
+++ b/src/spikeinterface/preprocessing/_decimation_tools.py
@@ -172,4 +172,7 @@ def get_antialiased_decimated_traces(
         )
     start_drop = left_margin // q
     n_out = end_frame - start_frame
-    return decimated_traces[start_drop : start_drop + n_out].astype(dtype)
+    decimated_traces = decimated_traces[start_drop : start_drop + n_out]
+    if np.issubdtype(np.dtype(dtype), np.integer):
+        np.round(decimated_traces, out=decimated_traces)  # Don't truncate towards zero
+    return decimated_traces.astype(dtype, copy=False)

From 2acd3a0ce19c22972acc315d01ab98d1860dfe75 Mon Sep 17 00:00:00 2001
From: Graham Findlay <graham.b.findlay@gmail.com>
Date: Thu, 10 Sep 2026 10:46:47 -0500
Subject: [PATCH 03/11] Share decimation helpers and fix sectionwise sample
 alignment

---
 .../preprocessing/_decimation_tools.py        | 123 ++++--------------
 src/spikeinterface/preprocessing/decimate.py  |  97 +++++++++++---
 src/spikeinterface/preprocessing/resample.py  |  64 ++++-----
 .../preprocessing/tests/test_decimate.py      |  61 ++++++---
 .../preprocessing/tests/test_resample.py      |  21 ++-
 5 files changed, 197 insertions(+), 169 deletions(-)

diff --git a/src/spikeinterface/preprocessing/_decimation_tools.py b/src/spikeinterface/preprocessing/_decimation_tools.py
index d6da9dd340..4b4bdeb9e3 100644
--- a/src/spikeinterface/preprocessing/_decimation_tools.py
+++ b/src/spikeinterface/preprocessing/_decimation_tools.py
@@ -7,12 +7,8 @@
 import math
 import warnings
 
-import numpy as np
-
-from spikeinterface.core import get_chunk_with_margin
-
 # scipy.signal.decimate uses an order-8 Chebyshev type I IIR filter by default, and its
-# documentation recommends decimating in several balanced steps rather than a single step
+# documentation recommends decimating in several steps rather than a single step
 # for downsampling factors larger than this value.
 _MAX_SINGLE_PASS_DECIMATION = 13
 
@@ -40,10 +36,10 @@ def _prime_factors(n):
     return factors
 
 
-def _greedy_pack(primes_desc, num_bins):
+def _greedy_pack(primes_desc, num_bins, max_factor=_MAX_SINGLE_PASS_DECIMATION):
     """
     Greedily pack `primes_desc` (largest first) into `num_bins` bins, keeping each bin's
-    product <= `_MAX_SINGLE_PASS_DECIMATION` and the bins as balanced as possible.
+    product <= `max_factor` and the bins as balanced as possible.
 
     Returns the list of bin products, or None if some prime cannot be placed (i.e. `num_bins`
     is too small to keep every bin <= the single-pass limit).
@@ -62,7 +58,7 @@ def _greedy_pack(primes_desc, num_bins):
     """
     bins = [1] * num_bins
     for prime in primes_desc:
-        fitting = [i for i in range(num_bins) if bins[i] * prime <= _MAX_SINGLE_PASS_DECIMATION]
+        fitting = [i for i in range(num_bins) if bins[i] * prime <= max_factor]
         if not fitting:
             return None
         # Place into the smallest fitting bin (ties broken by index, for determinism).
@@ -71,108 +67,41 @@ def _greedy_pack(primes_desc, num_bins):
     return bins
 
 
-def get_balanced_decimation_factors(decimation_factor):
+def get_balanced_decimation_factors(decimation_factor, max_factor=_MAX_SINGLE_PASS_DECIMATION):
     """
-    Split `decimation_factor` into sub-factors, each <= 13, as balanced as possible (so their
-    products are close), for stable multi-pass anti-aliased decimation.
+    Split `decimation_factor` into balanced sub-factors no greater than `max_factor`.
 
-    scipy recommends decimating in several balanced steps rather than one large step when the
-    factor exceeds 13 (e.g. 48 -> [8, 6] rather than [12, 4]). The product of the returned
-    factors always equals `decimation_factor`.
+    SciPy recommends multiple IIR decimation passes for factors above 13.
+    Balancing the factors (e.g. 48 -> [8, 6] rather than [12, 4]) further aids stability.
+    The product of the returned factors always equals `decimation_factor`.
 
     If `decimation_factor` has a prime factor greater than 13 (e.g. a large prime such as 17),
-    no valid split exists and `[decimation_factor]` is returned; it is the caller's
-    responsibility to handle this (e.g., warn that a single, potentially unstable, pass will
-    be used).
+    no valid split exists, a warning is issued, and `[decimation_factor]` is returned;
+    it is the caller's responsibility to handle this (e.g., warn that a single,
+    potentially unstable, pass will be used).
     """
-    if decimation_factor <= _MAX_SINGLE_PASS_DECIMATION:
+    if not isinstance(max_factor, int) or max_factor < 2:
+        raise ValueError("max_factor must be an integer greater than one")
+    if decimation_factor <= max_factor:
         return [decimation_factor]
 
     primes = _prime_factors(decimation_factor)
-    if max(primes) > _MAX_SINGLE_PASS_DECIMATION:
-        # If a prime factor > 13 cannot be split into sub-13 factors...
+    if max(primes) > max_factor:
+        warnings.warn(
+            f"`decimation_factor`={decimation_factor} cannot be split into anti-aliasing passes of <= {max_factor} "
+            f"(it has a prime factor > {max_factor}). A single `scipy.signal.decimate` pass will be used, "
+            f"which may be unstable. Consider a `decimation_factor` without large prime factors.",
+            stacklevel=2,
+        )
         return [decimation_factor]
 
     primes_desc = sorted(primes, reverse=True)
-    # Minimum number of passes so that, ideally, each pass decimates by <= 13.
-    num_passes = max(1, math.ceil(math.log(decimation_factor) / math.log(_MAX_SINGLE_PASS_DECIMATION)))
+    # Minimum number of passes so that, ideally, each pass decimates by <= max_factor.
+    num_passes = max(1, math.ceil(math.log(decimation_factor) / math.log(max_factor)))
     while num_passes <= len(primes_desc):
-        bins = _greedy_pack(primes_desc, num_passes)
+        bins = _greedy_pack(primes_desc, num_passes, max_factor)
         if bins is not None:
             return sorted(bins, reverse=True)
         num_passes += 1
-    # Fallback: one prime per pass (always valid since every prime is <= 13).
+    # Fallback: one prime per pass (always valid since every prime is <= max_factor).
     return primes_desc
-
-
-def get_antialiased_decimated_traces(
-    parent_segment,
-    start_frame,
-    end_frame,
-    channel_indices,
-    decimation_factor,
-    decimation_factors,
-    margin,
-    dtype,
-    decimation_offset=0,
-):
-    """
-    Fetch a margined chunk from `parent_segment` and decimate it by `decimation_factor`, applied
-    as a cascade of the balanced `decimation_factors` passes of ``scipy.signal.decimate``.
-
-    The margin is rounded up to a multiple of the total `decimation_factor` so that
-    ``left_margin // decimation_factor`` is exact; combined with scipy's default
-    ``zero_phase=True`` (output sample i maps to filtered input sample i * factor), this keeps the
-    downsampled traces aligned across chunks (a chunked read matches a full read). Exactly
-    ``end_frame - start_frame`` decimated samples are returned.
-
-    Parameters
-    ----------
-    parent_segment : BaseRecordingSegment
-        The parent segment to read (full-rate) traces from.
-    start_frame, end_frame : int
-        Output (decimated) frame range to return.
-    channel_indices : slice | list | np.ndarray | None
-        Channels to read, forwarded to the parent segment.
-    decimation_factor : int
-        The total decimation factor (the product of `decimation_factors`).
-    decimation_factors : list[int]
-        The per-pass sub-factors (each <= 13), e.g. from `get_balanced_decimation_factors`.
-    margin : int
-        Margin in parent samples used to limit anti-aliasing filter edge effects. Rounded up
-        internally to a multiple of `decimation_factor`.
-    dtype : np.dtype | str
-        Output dtype. The decimation runs in float32 and the result is cast to `dtype`.
-    decimation_offset : int, default: 0
-        Index of the first parent frame, applied to the first output sample only.
-    """
-    from scipy import signal
-
-    q = decimation_factor
-    parent_start_frame = decimation_offset + start_frame * q
-    parent_end_frame = parent_start_frame + (end_frame - start_frame) * q
-    # Round the margin up to a multiple of q so that left_margin // q is exact.
-    margin = int(np.ceil(margin / q) * q)
-    parent_traces, left_margin, right_margin = get_chunk_with_margin(
-        parent_segment,
-        parent_start_frame,
-        parent_end_frame,
-        channel_indices,
-        margin,
-        add_reflect_padding=True,
-        dtype=np.float32,
-    )
-    decimated_traces = parent_traces
-    for sub_q in decimation_factors:
-        decimated_traces = signal.decimate(decimated_traces, q=sub_q, axis=0)
-    if np.any(np.isnan(decimated_traces)):
-        warnings.warn(
-            f"`scipy.signal.decimate` produced NaNs while decimating by {q}. "
-            f"Consider a different decimation factor."
-        )
-    start_drop = left_margin // q
-    n_out = end_frame - start_frame
-    decimated_traces = decimated_traces[start_drop : start_drop + n_out]
-    if np.issubdtype(np.dtype(dtype), np.integer):
-        np.round(decimated_traces, out=decimated_traces)  # Don't truncate towards zero
-    return decimated_traces.astype(dtype, copy=False)
diff --git a/src/spikeinterface/preprocessing/decimate.py b/src/spikeinterface/preprocessing/decimate.py
index 348462ed6f..208f30acf8 100644
--- a/src/spikeinterface/preprocessing/decimate.py
+++ b/src/spikeinterface/preprocessing/decimate.py
@@ -7,12 +7,8 @@
 
 from .basepreprocessor import BasePreprocessor
 from .filter import fix_dtype
-from ._decimation_tools import (
-    _MAX_SINGLE_PASS_DECIMATION,
-    get_balanced_decimation_factors,
-    get_antialiased_decimated_traces,
-)
-from spikeinterface.core import BaseRecordingSegment
+from ._decimation_tools import get_balanced_decimation_factors
+from spikeinterface.core import BaseRecordingSegment, get_chunk_with_margin
 
 
 class DecimateRecording(BasePreprocessor):
@@ -91,13 +87,7 @@ def __init__(
         # fix_dtype doesn't always returns the str, make sure it does
         dtype = fix_dtype(recording, dtype).str
 
-        antialias_factors = get_balanced_decimation_factors(decimation_factor)
-        if antialias and decimation_factor > _MAX_SINGLE_PASS_DECIMATION and antialias_factors == [decimation_factor]:
-            warnings.warn(
-                f"`decimation_factor`={decimation_factor} cannot be split into anti-aliasing passes of <= 13 "
-                f"(it has a prime factor > 13). A single `scipy.signal.decimate` pass will be used, which may be "
-                f"unstable. Consider a `decimation_factor` without large prime factors."
-            )
+        decimation_factors = get_balanced_decimation_factors(decimation_factor) if antialias else None
 
         # Margin (in parent samples) to limit anti-aliasing filter edge effects.
         margin = int(margin_ms * self._orig_samp_freq / 1000)
@@ -115,7 +105,7 @@ def __init__(
                     self._dtype,
                     antialias,
                     margin,
-                    antialias_factors,
+                    decimation_factors,
                 )
             )
 
@@ -140,7 +130,7 @@ def __init__(
         dtype,
         antialias=False,
         margin=0,
-        antialias_factors=None,
+        decimation_factors=None,
     ):
         if parent_recording_segment._time_vector is not None:
             time_vector = parent_recording_segment._time_vector[decimation_offset::decimation_factor]
@@ -163,7 +153,7 @@ def __init__(
         self._dtype = dtype
         self._antialias = antialias
         self._margin = margin
-        self._antialias_factors = antialias_factors if antialias_factors is not None else [decimation_factor]
+        self._decimation_factors = decimation_factors if decimation_factors is not None else [decimation_factor]
 
     def get_num_samples(self):
         parent_n_samp = self._parent_segment.get_num_samples()
@@ -190,11 +180,84 @@ def get_traces(self, start_frame, end_frame, channel_indices):
             end_frame,
             channel_indices,
             self._decimation_factor,
-            self._antialias_factors,
+            self._decimation_factors,
             self._margin,
             self._dtype,
             decimation_offset=self._decimation_offset,
         )
 
 
+def get_antialiased_decimated_traces(
+    parent_segment,
+    start_frame,
+    end_frame,
+    channel_indices,
+    decimation_factor,
+    decimation_factors,
+    margin,
+    dtype,
+    decimation_offset=0,
+):
+    """
+    Fetch a margined chunk from `parent_segment` and decimate it by `decimation_factor`, applied
+    as a cascade of the balanced `decimation_factors` passes of ``scipy.signal.decimate``.
+
+    The margin is rounded up to a multiple of the total `decimation_factor` so that
+    ``left_margin // decimation_factor`` is exact; combined with scipy's default
+    ``zero_phase=True`` (output sample i maps to filtered input sample i * factor), this keeps the
+    downsampled traces aligned across chunks (a chunked read matches a full read). Exactly
+    ``end_frame - start_frame`` decimated samples are returned.
+
+    Parameters
+    ----------
+    parent_segment : BaseRecordingSegment
+        The parent segment to read (full-rate) traces from.
+    start_frame, end_frame : int
+        Output (decimated) frame range to return.
+    channel_indices : slice | list | np.ndarray | None
+        Channels to read, forwarded to the parent segment.
+    decimation_factor : int
+        The total decimation factor (the product of `decimation_factors`).
+    decimation_factors : list[int]
+        The per-pass sub-factors (each <= 13), e.g. from `get_balanced_decimation_factors`.
+    margin : int
+        Margin in parent samples used to limit anti-aliasing filter edge effects. Rounded up
+        internally to a multiple of `decimation_factor`.
+    dtype : np.dtype | str
+        Output dtype. The decimation runs in float32 and the result is cast to `dtype`.
+    decimation_offset : int, default: 0
+        Index of the first parent frame, applied to the first output sample only.
+    """
+    from scipy import signal
+
+    q = decimation_factor
+    parent_start_frame = decimation_offset + start_frame * q
+    parent_end_frame = parent_start_frame + (end_frame - start_frame) * q
+    # Round the margin up to a multiple of q so that left_margin // q is exact.
+    margin = int(np.ceil(margin / q) * q)
+    parent_traces, left_margin, right_margin = get_chunk_with_margin(
+        parent_segment,
+        parent_start_frame,
+        parent_end_frame,
+        channel_indices,
+        margin,
+        add_reflect_padding=True,
+        dtype=np.float32,
+    )
+    decimated_traces = parent_traces
+    for sub_q in decimation_factors:
+        decimated_traces = signal.decimate(decimated_traces, q=sub_q, axis=0)
+    if np.any(np.isnan(decimated_traces)):
+        warnings.warn(
+            f"`scipy.signal.decimate` produced NaNs while decimating by {q}. "
+            f"Consider a different decimation factor."
+        )
+    start_drop = left_margin // q
+    n_out = end_frame - start_frame
+    decimated_traces = decimated_traces[start_drop : start_drop + n_out]
+    if np.issubdtype(np.dtype(dtype), np.integer):
+        np.round(decimated_traces, out=decimated_traces)  # Don't truncate towards zero
+    return decimated_traces.astype(dtype, copy=False)
+
+
 decimate = define_function_handling_dict_from_class(source_class=DecimateRecording, name="decimate")
diff --git a/src/spikeinterface/preprocessing/resample.py b/src/spikeinterface/preprocessing/resample.py
index d6bf527287..46e883a685 100644
--- a/src/spikeinterface/preprocessing/resample.py
+++ b/src/spikeinterface/preprocessing/resample.py
@@ -8,11 +8,8 @@
 
 from .basepreprocessor import BasePreprocessor
 from .filter import fix_dtype
-from ._decimation_tools import (
-    _MAX_SINGLE_PASS_DECIMATION,
-    get_balanced_decimation_factors,
-    get_antialiased_decimated_traces,
-)
+from ._decimation_tools import get_balanced_decimation_factors
+from .decimate import get_antialiased_decimated_traces
 from spikeinterface.core import get_chunk_with_margin, BaseRecordingSegment
 
 
@@ -56,7 +53,7 @@ class ResampleRecording(BasePreprocessor):
         - 1.0: Tolerate gaps up to 1 ms, split on larger gaps
         - 100.0: Only major pauses (>100 ms) create sections
     margin_ms : float, default: 100.0
-        Margin in ms for computations, will be used to decrease edge effects.
+        Margin in ms for computations, used to decrease edge effects.
     dtype : dtype or None, default: None
         The dtype of the returned traces. If None, the dtype of the parent recording is used.
     skip_checks : bool, default: False
@@ -94,19 +91,10 @@ def __init__(
         self._orig_samp_freq = recording.get_sampling_frequency()
         self._resample_rate = resample_rate
         self._sampling_frequency = resample_rate
-        # When the parent rate is an exact integer multiple of resample_rate, get_traces uses
-        # anti-aliased decimation. Large factors are split into several balanced sub-13 passes
-        # (see _decimation_tools); only an unsplittable factor (a prime > 13) falls back to a
-        # single, potentially unstable, pass and warns.
+        # Exact integer-factor downsampling uses one or more antialiased decimation passes.
         if self._orig_samp_freq % resample_rate == 0:
             decimation_factor = int(self._orig_samp_freq / resample_rate)
             decimation_factors = get_balanced_decimation_factors(decimation_factor)
-            if decimation_factors == [decimation_factor] and decimation_factor > _MAX_SINGLE_PASS_DECIMATION:
-                warnings.warn(
-                    f"Resampling by an integer factor of {decimation_factor} cannot be split into "
-                    f"anti-aliasing passes of <= 13 (it has a prime factor > 13); a single "
-                    f"`scipy.signal.decimate` pass will be used, which may be unstable."
-                )
         else:
             decimation_factors = None
         # fix_dtype not always returns the str, make sure it does
@@ -285,7 +273,7 @@ def get_traces(self, start_frame, end_frame, channel_indices):
             return self._get_traces_gapped(start_frame, end_frame, channel_indices)
 
         if self._decimation_factors is not None:
-            # Integer ratio, no gaps: anti-aliased multi-pass decimation.
+            # Integer-factor downsampling, no gaps: one or more antialiased decimation passes.
             decimation_factor = int(self._parent_rate / self._resample_rate)
             return get_antialiased_decimated_traces(
                 self._parent_segment,
@@ -298,7 +286,7 @@ def get_traces(self, start_frame, end_frame, channel_indices):
                 self._dtype,
             )
 
-        # Non-integer ratio, no gaps: FFT-based resampling with proportional margins.
+        # Non-integer-ratio downsampling or upsampling, no gaps: FFT with proportional margins.
         from scipy import signal
 
         parent_start_frame, parent_end_frame = [
@@ -350,7 +338,13 @@ def _get_traces_gapped(self, start_frame, end_frame, channel_indices):
         first_sec = max(first_sec, 0)
         last_sec = min(last_sec, len(self._sec_n_out) - 1)
 
-        is_integer_ratio = (self._parent_rate % self._resample_rate) == 0
+        is_integer_ratio = self._decimation_factors is not None
+        margin = self._margin
+        if is_integer_ratio:
+            # Integer-factor decimation rounds the margin up to a multiple of the total
+            # factor, in parent samples.
+            q = int(self._parent_rate / self._resample_rate)
+            margin = ((margin + q - 1) // q) * q
 
         pos = 0
         for k in range(first_sec, last_sec + 1):
@@ -371,15 +365,19 @@ def _get_traces_gapped(self, start_frame, end_frame, channel_indices):
             if local_out_end <= local_out_start:
                 continue
 
-            # Map within-section output frames to within-section parent frames
-            local_par_start = int((local_out_start / self._resample_rate) * self._parent_rate)
-            local_par_end = int((local_out_end / self._resample_rate) * self._parent_rate)
+            # Map within-section output frames to within-section parent frames.
+            if is_integer_ratio:
+                local_par_start = local_out_start * q
+                local_par_end = local_out_end * q
+            else:
+                local_par_start = int((local_out_start / self._resample_rate) * self._parent_rate)
+                local_par_end = int((local_out_end / self._resample_rate) * self._parent_rate)
             local_par_start = max(0, min(local_par_start, sec_n_parent))
             local_par_end = max(0, min(local_par_end, sec_n_parent))
 
             # Apply margin within section boundaries only (do not cross gaps)
-            left_margin = min(self._margin, local_par_start)
-            right_margin = min(self._margin, sec_n_parent - local_par_end)
+            left_margin = min(margin, local_par_start)
+            right_margin = min(margin, sec_n_parent - local_par_end)
 
             par_fetch_start = par_start_k + local_par_start - left_margin
             par_fetch_end = par_start_k + local_par_end + right_margin
@@ -390,16 +388,20 @@ def _get_traces_gapped(self, start_frame, end_frame, channel_indices):
             )
 
             # Apply reflect padding if margin was truncated at section edge
-            pad_left = self._margin - left_margin
-            pad_right = self._margin - right_margin
+            pad_left = margin - left_margin
+            pad_right = margin - right_margin
             if pad_left > 0 or pad_right > 0:
                 parent_traces = np.pad(parent_traces, [(pad_left, pad_right), (0, 0)], mode="reflect")
-                left_margin = self._margin
-                right_margin = self._margin
+                left_margin = margin
+                right_margin = margin
 
-            # Compute resampled margins
-            left_margin_rs = int((left_margin / self._parent_rate) * self._resample_rate)
-            right_margin_rs = int((right_margin / self._parent_rate) * self._resample_rate)
+            # Compute resampled margins on the output grid.
+            if is_integer_ratio:
+                left_margin_rs = left_margin // q
+                right_margin_rs = right_margin // q
+            else:
+                left_margin_rs = int((left_margin / self._parent_rate) * self._resample_rate)
+                right_margin_rs = int((right_margin / self._parent_rate) * self._resample_rate)
 
             # Total output samples including margins
             chunk_len = int(local_out_end - local_out_start)
diff --git a/src/spikeinterface/preprocessing/tests/test_decimate.py b/src/spikeinterface/preprocessing/tests/test_decimate.py
index bdc880887d..b17b616a31 100644
--- a/src/spikeinterface/preprocessing/tests/test_decimate.py
+++ b/src/spikeinterface/preprocessing/tests/test_decimate.py
@@ -1,9 +1,12 @@
+import warnings
+
 import pytest
 
 
 from spikeinterface import NumpyRecording
 from spikeinterface.core import generate_recording, load
 from spikeinterface.preprocessing.decimate import DecimateRecording, decimate, get_balanced_decimation_factors
+from spikeinterface.preprocessing.resample import ResampleRecording
 from spikeinterface.preprocessing.tests.test_resample import create_sinusoidal_traces
 import numpy as np
 
@@ -79,27 +82,34 @@ def test_decimate_with_times(antialias):
 
 
 @pytest.mark.parametrize(
-    "decimation_factor, expected",
+    "decimation_factor, max_factor, expected",
     [
-        (1, [1]),
-        (7, [7]),
-        (13, [13]),
-        (48, [8, 6]),
-        (50, [10, 5]),
-        (60, [10, 6]),
-        (100, [10, 10]),
-        (17, [17]),  # prime > 13: cannot be split
-        (23, [23]),  # prime > 13: cannot be split
+        (1, 13, [1]),
+        (7, 13, [7]),
+        (13, 13, [13]),
+        (48, 13, [8, 6]),
+        (50, 13, [10, 5]),
+        (60, 13, [10, 6]),
+        (100, 13, [10, 10]),
+        (17, 13, [17]),  # prime > 13: cannot be split
+        (23, 13, [23]),
+        (48, 8, [8, 6]),
+        (48, 4, [4, 4, 3]),
+        (10, 4, [10]),  # prime > the custom limit: cannot be split
     ],
 )
-def test_balanced_decimation_factors(decimation_factor, expected):
-    factors = get_balanced_decimation_factors(decimation_factor)
+def test_balanced_decimation_factors(decimation_factor, max_factor, expected):
+    if expected == [decimation_factor] and decimation_factor > max_factor:
+        with pytest.warns(UserWarning, match=f"prime factor > {max_factor}"):
+            factors = get_balanced_decimation_factors(decimation_factor, max_factor=max_factor)
+    else:
+        factors = get_balanced_decimation_factors(decimation_factor, max_factor=max_factor)
     assert factors == expected
     # The product of the sub-factors always reconstructs the requested factor.
     assert int(np.prod(factors)) == decimation_factor
-    # Every pass is <= 13 unless the factor is an unsplittable prime > 13.
+    # Every pass respects the limit unless no valid split exists.
     if len(factors) > 1:
-        assert all(f <= 13 for f in factors)
+        assert all(f <= max_factor for f in factors)
 
 
 @pytest.mark.parametrize("decimation_factor", [6, 10, 48])
@@ -169,8 +179,8 @@ def test_decimate_antialias_multipass():
     assert dec._kwargs["decimation_factor"] == decimation_factor
 
     segment = dec.segments[0]
-    assert int(np.prod(segment._antialias_factors)) == decimation_factor
-    assert all(f <= 13 for f in segment._antialias_factors)
+    assert int(np.prod(segment._decimation_factors)) == decimation_factor
+    assert all(f <= 13 for f in segment._decimation_factors)
 
     parent_n = parent_rec.get_num_samples()
     assert dec.get_num_samples() == int(np.ceil(parent_n / decimation_factor))
@@ -181,11 +191,20 @@ def test_decimate_antialias_multipass():
 
 
 def test_decimate_antialias_large_prime_warns():
-    rec = generate_recording(durations=[2.0], num_channels=2)
-    with pytest.warns(UserWarning, match="prime factor > 13"):
-        dec = DecimateRecording(rec, 17, antialias=True)
-    # The unsplittable factor falls back to a single pass.
-    assert dec.segments[0]._antialias_factors == [17]
+    rec = generate_recording(durations=[2.0], num_channels=2, sampling_frequency=34000)
+    for preprocess in [
+        lambda: DecimateRecording(rec, 17, antialias=True),
+        lambda: ResampleRecording(rec, rec.get_sampling_frequency() / 17),
+    ]:
+        with pytest.warns(UserWarning, match="prime factor > 13") as caught:
+            dec = preprocess()
+        assert len(caught) == 1
+        assert dec.segments[0]._decimation_factors == [17]
+
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always")
+        DecimateRecording(rec, 17, antialias=False)
+    assert not caught
 
 
 if __name__ == "__main__":
diff --git a/src/spikeinterface/preprocessing/tests/test_resample.py b/src/spikeinterface/preprocessing/tests/test_resample.py
index 089905aaee..1d9a59c88a 100644
--- a/src/spikeinterface/preprocessing/tests/test_resample.py
+++ b/src/spikeinterface/preprocessing/tests/test_resample.py
@@ -360,8 +360,9 @@ def test_resample_preserves_gaps_non_integer_ratio():
     [
         700,  # non-integer ratio (30000 / 700 ~= 42.857)
         500,  # integer ratio (30000 / 500 = 60)
+        625,  # factor 48 does not divide the default 3000-sample margin
     ],
-    ids=["non_integer_ratio", "integer_ratio"],
+    ids=["non_integer_ratio", "integer_ratio", "integer_ratio_unaligned_margin"],
 )
 def test_resample_traces_across_gap(resample_rate):
     """Section-wise resampling should match individually resampled sections.
@@ -436,11 +437,25 @@ def test_resample_traces_across_gap(resample_rate):
     assert gapped_s2.shape == ref_traces2.shape, f"Section 2 shape mismatch: {gapped_s2.shape} vs {ref_traces2.shape}"
     np.testing.assert_allclose(gapped_s2, ref_traces2, rtol=1e-5, atol=1e-5)
 
+    if sampling_frequency % resample_rate == 0:
+        # Exercise reads within a section and across the gap, including channel selection.
+        for start, end in [(123, 456), (n_out_1 - 123, n_out_1 + 456)]:
+            expected_pieces = []
+            for offset, section in [(0, resampled1), (n_out_1, resampled2)]:
+                local_start = max(0, start - offset)
+                local_end = min(section.get_num_samples(), end - offset)
+                if local_start < local_end:
+                    expected_pieces.append(
+                        section.get_traces(start_frame=local_start, end_frame=local_end, channel_ids=[1])
+                    )
+            actual = resampled.get_traces(start_frame=start, end_frame=end, channel_ids=[1])
+            np.testing.assert_allclose(actual, np.concatenate(expected_pieces), rtol=1e-5, atol=1e-5)
+
 
-def test_resample_gapped_chunked_consistency():
+@pytest.mark.parametrize("resample_rate", [700, 625])
+def test_resample_gapped_chunked_consistency(resample_rate):
     """Chunked .save() should match non-chunked for gapped recordings."""
     sampling_frequency = 30000
-    resample_rate = 700
     rec, _, _ = _make_gapped_recording(sampling_frequency=sampling_frequency, sec1_duration=2.0, sec2_duration=2.0)
 
     import warnings as _warnings

From 4e9506e6b1ac35395c1ede3ff60547354a5a6cac Mon Sep 17 00:00:00 2001
From: Graham Findlay <graham.b.findlay@gmail.com>
Date: Thu, 10 Sep 2026 11:15:44 -0500
Subject: [PATCH 04/11] Unify sectionwise resampling and raise on NaN output.

Reuses the same trace processing for gapped sections and standalone
recordings, removing duplicate padding logic and the sneaky silent FFT retry.
---
 src/spikeinterface/preprocessing/decimate.py |  44 +++++--
 src/spikeinterface/preprocessing/resample.py | 131 ++++---------------
 2 files changed, 54 insertions(+), 121 deletions(-)

diff --git a/src/spikeinterface/preprocessing/decimate.py b/src/spikeinterface/preprocessing/decimate.py
index 208f30acf8..7227f5cb82 100644
--- a/src/spikeinterface/preprocessing/decimate.py
+++ b/src/spikeinterface/preprocessing/decimate.py
@@ -1,5 +1,3 @@
-import warnings
-
 import numpy as np
 from spikeinterface.core.core_tools import (
     define_function_handling_dict_from_class,
@@ -224,7 +222,7 @@ def get_antialiased_decimated_traces(
         Margin in parent samples used to limit anti-aliasing filter edge effects. Rounded up
         internally to a multiple of `decimation_factor`.
     dtype : np.dtype | str
-        Output dtype. The decimation runs in float32 and the result is cast to `dtype`.
+        Output dtype. Integer output is rounded and clipped to its range.
     decimation_offset : int, default: 0
         Index of the first parent frame, applied to the first output sample only.
     """
@@ -234,7 +232,7 @@ def get_antialiased_decimated_traces(
     parent_start_frame = decimation_offset + start_frame * q
     parent_end_frame = parent_start_frame + (end_frame - start_frame) * q
     # Round the margin up to a multiple of q so that left_margin // q is exact.
-    margin = int(np.ceil(margin / q) * q)
+    margin = ((margin + q - 1) // q) * q
     parent_traces, left_margin, right_margin = get_chunk_with_margin(
         parent_segment,
         parent_start_frame,
@@ -242,22 +240,40 @@ def get_antialiased_decimated_traces(
         channel_indices,
         margin,
         add_reflect_padding=True,
-        dtype=np.float32,
     )
-    decimated_traces = parent_traces
+    working_dtype = np.result_type(parent_traces.dtype, dtype, np.float32)
+    decimated_traces = parent_traces.astype(working_dtype, copy=False)
     for sub_q in decimation_factors:
         decimated_traces = signal.decimate(decimated_traces, q=sub_q, axis=0)
-    if np.any(np.isnan(decimated_traces)):
-        warnings.warn(
-            f"`scipy.signal.decimate` produced NaNs while decimating by {q}. "
-            f"Consider a different decimation factor."
-        )
     start_drop = left_margin // q
     n_out = end_frame - start_frame
     decimated_traces = decimated_traces[start_drop : start_drop + n_out]
-    if np.issubdtype(np.dtype(dtype), np.integer):
-        np.round(decimated_traces, out=decimated_traces)  # Don't truncate towards zero
-    return decimated_traces.astype(dtype, copy=False)
+    return _cast_resampled_traces(decimated_traces, dtype)
+
+
+def _cast_resampled_traces(traces, dtype):
+    """Reject nonfinite output and round and saturate integer conversions."""
+    if not np.all(np.isfinite(traces)):
+        raise ValueError("Resampling produced nonfinite values. Check the input traces and resampling parameters.")
+
+    dtype = np.dtype(dtype)
+    if np.issubdtype(dtype, np.integer):
+        rounded = np.rint(traces)
+        limits = np.iinfo(dtype)
+        below = rounded <= limits.min
+        above = rounded >= limits.max
+
+        # Assign saturated endpoints after casting because apparently float64 can't 
+        # represent int64.max exactly. 
+        rounded[below | above] = 0
+        result = rounded.astype(dtype)
+        result[below] = limits.min
+        result[above] = limits.max
+        return result
+
+    if np.issubdtype(dtype, np.floating) and np.any(np.abs(traces) > np.finfo(dtype).max):
+        raise ValueError(f"Resampled values exceed the finite range of {dtype}.")
+    return traces.astype(dtype, copy=False)
 
 
 decimate = define_function_handling_dict_from_class(source_class=DecimateRecording, name="decimate")
diff --git a/src/spikeinterface/preprocessing/resample.py b/src/spikeinterface/preprocessing/resample.py
index 46e883a685..cf8c7b7024 100644
--- a/src/spikeinterface/preprocessing/resample.py
+++ b/src/spikeinterface/preprocessing/resample.py
@@ -9,8 +9,9 @@
 from .basepreprocessor import BasePreprocessor
 from .filter import fix_dtype
 from ._decimation_tools import get_balanced_decimation_factors
-from .decimate import get_antialiased_decimated_traces
+from .decimate import get_antialiased_decimated_traces, _cast_resampled_traces
 from spikeinterface.core import get_chunk_with_margin, BaseRecordingSegment
+from spikeinterface.core.frameslicerecording import FrameSliceRecordingSegment
 
 
 class ResampleRecording(BasePreprocessor):
@@ -56,6 +57,8 @@ class ResampleRecording(BasePreprocessor):
         Margin in ms for computations, used to decrease edge effects.
     dtype : dtype or None, default: None
         The dtype of the returned traces. If None, the dtype of the parent recording is used.
+        Integer output is rounded and clipped to the dtype range. Nonfinite resampled
+        output raises a ValueError before conversion.
     skip_checks : bool, default: False
         If True, checks on sampling frequencies and cutoff filter frequencies are skipped
 
@@ -272,11 +275,14 @@ def get_traces(self, start_frame, end_frame, channel_indices):
         if self._has_gaps:
             return self._get_traces_gapped(start_frame, end_frame, channel_indices)
 
+        return self._get_resampled_traces(self._parent_segment, start_frame, end_frame, channel_indices)
+
+    def _get_resampled_traces(self, parent_segment, start_frame, end_frame, channel_indices):
         if self._decimation_factors is not None:
             # Integer-factor downsampling, no gaps: one or more antialiased decimation passes.
             decimation_factor = int(self._parent_rate / self._resample_rate)
             return get_antialiased_decimated_traces(
-                self._parent_segment,
+                parent_segment,
                 start_frame,
                 end_frame,
                 channel_indices,
@@ -293,14 +299,15 @@ def get_traces(self, start_frame, end_frame, channel_indices):
             int((frame / self._resample_rate) * self._parent_rate) for frame in [start_frame, end_frame]
         ]
         parent_traces, left_margin, right_margin = get_chunk_with_margin(
-            self._parent_segment,
+            parent_segment,
             parent_start_frame,
             parent_end_frame,
             channel_indices,
             self._margin,
             add_reflect_padding=True,
-            dtype=np.float32,
         )
+        working_dtype = np.result_type(parent_traces.dtype, self._dtype, np.float32)
+        parent_traces = parent_traces.astype(working_dtype, copy=False)
         # get left and right margins for the resampled case
         left_margin_rs, right_margin_rs = [
             int((margin / self._parent_rate) * self._resample_rate) for margin in [left_margin, right_margin]
@@ -312,123 +319,33 @@ def get_traces(self, start_frame, end_frame, channel_indices):
 
         # now take care of the edges
         resampled_traces = resampled_traces[left_margin_rs : num - right_margin_rs]
-        return resampled_traces.astype(self._dtype)
+        return _cast_resampled_traces(resampled_traces, self._dtype)
 
     def _get_traces_gapped(self, start_frame, end_frame, channel_indices):
-        """Resample traces section-by-section, avoiding FFT processing across gaps."""
-        from scipy.signal import decimate, resample
-
-        # Determine the post-indexing channel count via a 1-sample parent fetch.
-        # channel_indices may be a slice, list, ndarray, or None, so we cannot
-        # simply use len(channel_indices).
+        """Resample each section with margins bounded by its own samples."""
         n_channels = self._parent_segment.get_traces(0, 1, channel_indices).shape[1]
-
-        # Pre-allocate the output buffer.
         result = np.empty((end_frame - start_frame, n_channels), dtype=self._dtype)
-
         if start_frame == end_frame:
             return result
 
-        # Find which sections overlap [start_frame, end_frame) in output space.
-        # _sec_boundaries_output[k] = [out_start_k, out_end_k)
-        sec_ends = self._sec_boundaries_output[:, 1]
         sec_starts = self._sec_boundaries_output[:, 0]
+        sec_ends = self._sec_boundaries_output[:, 1]
         first_sec = int(np.searchsorted(sec_ends, start_frame, side="right"))
-        last_sec = int(np.searchsorted(sec_starts, end_frame, side="left")) - 1
-        first_sec = max(first_sec, 0)
-        last_sec = min(last_sec, len(self._sec_n_out) - 1)
-
-        is_integer_ratio = self._decimation_factors is not None
-        margin = self._margin
-        if is_integer_ratio:
-            # Integer-factor decimation rounds the margin up to a multiple of the total
-            # factor, in parent samples.
-            q = int(self._parent_rate / self._resample_rate)
-            margin = ((margin + q - 1) // q) * q
-
-        pos = 0
-        for k in range(first_sec, last_sec + 1):
-            out_start_k = int(self._sec_boundaries_output[k, 0])
-            out_end_k = int(self._sec_boundaries_output[k, 1])
-            par_start_k = int(self._sec_boundaries_parent[k, 0])
-            par_end_k = int(self._sec_boundaries_parent[k, 1])
-            sec_n_parent = par_end_k - par_start_k
-            sec_n_output = int(self._sec_n_out[k])
-
-            if sec_n_output == 0:
-                continue
+        stop_sec = int(np.searchsorted(sec_starts, end_frame, side="left"))
 
-            # Clip the output range to the requested [start_frame, end_frame)
-            local_out_start = max(start_frame, out_start_k) - out_start_k
-            local_out_end = min(end_frame, out_end_k) - out_start_k
-
-            if local_out_end <= local_out_start:
+        for k in range(first_sec, stop_sec):
+            out_start = max(start_frame, int(sec_starts[k]))
+            out_end = min(end_frame, int(sec_ends[k]))
+            if out_start >= out_end:
                 continue
 
-            # Map within-section output frames to within-section parent frames.
-            if is_integer_ratio:
-                local_par_start = local_out_start * q
-                local_par_end = local_out_end * q
-            else:
-                local_par_start = int((local_out_start / self._resample_rate) * self._parent_rate)
-                local_par_end = int((local_out_end / self._resample_rate) * self._parent_rate)
-            local_par_start = max(0, min(local_par_start, sec_n_parent))
-            local_par_end = max(0, min(local_par_end, sec_n_parent))
-
-            # Apply margin within section boundaries only (do not cross gaps)
-            left_margin = min(margin, local_par_start)
-            right_margin = min(margin, sec_n_parent - local_par_end)
-
-            par_fetch_start = par_start_k + local_par_start - left_margin
-            par_fetch_end = par_start_k + local_par_end + right_margin
-
-            # Fetch parent traces for this section's sub-chunk
-            parent_traces = self._parent_segment.get_traces(par_fetch_start, par_fetch_end, channel_indices).astype(
-                np.float32
+            par_start, par_end = self._sec_boundaries_parent[k]
+            section = FrameSliceRecordingSegment(self._parent_segment, int(par_start), int(par_end))
+            result[out_start - start_frame : out_end - start_frame] = self._get_resampled_traces(
+                section, out_start - int(sec_starts[k]), out_end - int(sec_starts[k]), channel_indices
             )
 
-            # Apply reflect padding if margin was truncated at section edge
-            pad_left = margin - left_margin
-            pad_right = margin - right_margin
-            if pad_left > 0 or pad_right > 0:
-                parent_traces = np.pad(parent_traces, [(pad_left, pad_right), (0, 0)], mode="reflect")
-                left_margin = margin
-                right_margin = margin
-
-            # Compute resampled margins on the output grid.
-            if is_integer_ratio:
-                left_margin_rs = left_margin // q
-                right_margin_rs = right_margin // q
-            else:
-                left_margin_rs = int((left_margin / self._parent_rate) * self._resample_rate)
-                right_margin_rs = int((right_margin / self._parent_rate) * self._resample_rate)
-
-            # Total output samples including margins
-            chunk_len = int(local_out_end - local_out_start)
-            num = chunk_len + left_margin_rs + right_margin_rs
-
-            # Resample this section. Integer ratios use anti-aliased multi-pass decimation
-            # (the balanced sub-13 passes computed in ResampleRecording.__init__), applied
-            # within the section only so filtering never crosses a gap.
-            if is_integer_ratio:
-                resampled = parent_traces
-                for sub_q in self._decimation_factors:
-                    resampled = decimate(resampled, q=sub_q, axis=0)
-                if np.any(np.isnan(resampled)):
-                    resampled = resample(parent_traces, num, axis=0)
-            else:
-                resampled = resample(parent_traces, num, axis=0)
-
-            # Trim margins and write directly into the pre-allocated buffer.
-            # Clamp to the remaining space in case decimate's output length
-            # differs from `num` by a rounding sample.
-            trimmed = resampled[left_margin_rs : num - right_margin_rs]
-            write_len = min(len(trimmed), result.shape[0] - pos)
-            result[pos : pos + write_len] = trimmed[:write_len]
-            pos += write_len
-
-        # Return only the filled portion (normally equals end_frame - start_frame).
-        return result[:pos]
+        return result
 
 
 resample = define_function_handling_dict_from_class(source_class=ResampleRecording, name="resample")

From 160bfc8fe38b6cfb934c2c265c28c5d61961a16a Mon Sep 17 00:00:00 2001
From: Graham Findlay <graham.b.findlay@gmail.com>
Date: Thu, 10 Sep 2026 11:56:41 -0500
Subject: [PATCH 05/11] Warn of impending change to the default decimation
 antialiasing behavior

---
 doc/modules/preprocessing.rst                |  8 ++++++
 src/spikeinterface/preprocessing/decimate.py | 29 +++++++++++++++-----
 2 files changed, 30 insertions(+), 7 deletions(-)

diff --git a/doc/modules/preprocessing.rst b/doc/modules/preprocessing.rst
index e587ac45a5..577cef8588 100644
--- a/doc/modules/preprocessing.rst
+++ b/doc/modules/preprocessing.rst
@@ -53,6 +53,14 @@ CMR, and save it to a binary file in the "/path/to/preprocessed" folder. The :co
 
 **NOTE:** some sorters will automatically perform the saving operation internally.
 
+.. note::
+
+    ``decimate()`` currently uses slicing without antialiasing when ``antialias`` is omitted,
+    and emits a ``FutureWarning``. A future release will enable antialiasing by default.
+    Set ``antialias=True`` to filter before downsampling, or explicitly set ``antialias=False``
+    to retain slicing without the transition warning, for example when the recording has
+    already been sufficiently low-pass filtered for the output sampling rate.
+
 The Preprocessing Pipeline
 --------------------------
 
diff --git a/src/spikeinterface/preprocessing/decimate.py b/src/spikeinterface/preprocessing/decimate.py
index 7227f5cb82..217b347022 100644
--- a/src/spikeinterface/preprocessing/decimate.py
+++ b/src/spikeinterface/preprocessing/decimate.py
@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 from spikeinterface.core.core_tools import (
     define_function_handling_dict_from_class,
@@ -32,13 +34,15 @@ class DecimateRecording(BasePreprocessor):
         to ensure that the decimated recording has at least one frame. Consider combining DecimateRecording
         with FrameSliceRecording for fine control on the recording start and end frames.
         The same decimation offset is applied to all segments from the parent recording.
-    antialias : bool, default: False
+    antialias : bool | None, default: None
         If True, apply an anti-aliasing low-pass filter before downsampling, using
         ``scipy.signal.decimate``. When `decimation_factor` exceeds 13, the decimation is
         automatically performed in several balanced sub-13 passes (e.g. a factor of 48 is applied
-        as 8 then 6), as scipy recommends, to keep the IIR anti-aliasing filter stable. If False
-        (the default), traces are downsampled by plain array slicing with no filtering, and
-        `margin_ms` is ignored.
+        as 8 then 6), as scipy recommends, to keep the IIR anti-aliasing filter stable. If False,
+        traces are downsampled by plain array slicing with no filtering, and `margin_ms`
+        is ignored. If omitted or None, currently behaves as False and emits a FutureWarning:
+        a future release will enable antialiasing by default. Pass True or False explicitly
+        to select the behavior and silence the transition warning.
     margin_ms : float, default: 100.0
         Margin in ms used on each side of every chunk to limit edge effects of the anti-aliasing
         filter. Only used when `antialias=True`. The margin is internally rounded up to a whole
@@ -60,7 +64,7 @@ def __init__(
         recording,
         decimation_factor,
         decimation_offset=0,
-        antialias=False,
+        antialias=None,
         margin_ms=100.0,
         dtype=None,
     ):
@@ -85,6 +89,17 @@ def __init__(
         # fix_dtype doesn't always returns the str, make sure it does
         dtype = fix_dtype(recording, dtype).str
 
+        if antialias is None:
+            warnings.warn(
+                "The default for `antialias` will change to True in a future release. "
+                "Currently, decimation uses slicing without an anti-aliasing filter. "
+                "Pass antialias=True to enable the filter, "
+                "or antialias=False to explicitly retain slicing.",
+                FutureWarning,
+                stacklevel=2,
+            )
+            antialias = False
+
         decimation_factors = get_balanced_decimation_factors(decimation_factor) if antialias else None
 
         # Margin (in parent samples) to limit anti-aliasing filter edge effects.
@@ -263,8 +278,8 @@ def _cast_resampled_traces(traces, dtype):
         below = rounded <= limits.min
         above = rounded >= limits.max
 
-        # Assign saturated endpoints after casting because apparently float64 can't 
-        # represent int64.max exactly. 
+        # Assign saturated endpoints after casting because apparently float64 can't
+        # represent int64.max exactly.
         rounded[below | above] = 0
         result = rounded.astype(dtype)
         result[below] = limits.min

From f0e5bdae61465f0d12d4ff7418ca52295213b2c7 Mon Sep 17 00:00:00 2001
From: Graham Findlay <graham.b.findlay@gmail.com>
Date: Thu, 10 Sep 2026 12:20:04 -0500
Subject: [PATCH 06/11] Remove note in preprocessing docs about future
 antialiasing default change.

I decided that this doc is more of a tutorial than a reference, so this doesn't really belong here.
---
 doc/modules/preprocessing.rst | 7 -------
 1 file changed, 7 deletions(-)

diff --git a/doc/modules/preprocessing.rst b/doc/modules/preprocessing.rst
index 577cef8588..97958cd4d0 100644
--- a/doc/modules/preprocessing.rst
+++ b/doc/modules/preprocessing.rst
@@ -53,13 +53,6 @@ CMR, and save it to a binary file in the "/path/to/preprocessed" folder. The :co
 
 **NOTE:** some sorters will automatically perform the saving operation internally.
 
-.. note::
-
-    ``decimate()`` currently uses slicing without antialiasing when ``antialias`` is omitted,
-    and emits a ``FutureWarning``. A future release will enable antialiasing by default.
-    Set ``antialias=True`` to filter before downsampling, or explicitly set ``antialias=False``
-    to retain slicing without the transition warning, for example when the recording has
-    already been sufficiently low-pass filtered for the output sampling rate.
 
 The Preprocessing Pipeline
 --------------------------

From f987d85b338c0e07acc29a837bc13e41e3e653bc Mon Sep 17 00:00:00 2001
From: Graham Findlay <graham.b.findlay@gmail.com>
Date: Thu, 10 Sep 2026 12:20:51 -0500
Subject: [PATCH 07/11] Automatic decimation margin estimation

---
 .../preprocessing/_decimation_tools.py        | 38 +++++++++++++++++++
 src/spikeinterface/preprocessing/decimate.py  | 10 +++--
 src/spikeinterface/preprocessing/resample.py  | 13 ++++---
 .../preprocessing/tests/test_decimate.py      |  4 +-
 4 files changed, 54 insertions(+), 11 deletions(-)

diff --git a/src/spikeinterface/preprocessing/_decimation_tools.py b/src/spikeinterface/preprocessing/_decimation_tools.py
index 4b4bdeb9e3..e6b4aed6e6 100644
--- a/src/spikeinterface/preprocessing/_decimation_tools.py
+++ b/src/spikeinterface/preprocessing/_decimation_tools.py
@@ -13,6 +13,44 @@
 _MAX_SINGLE_PASS_DECIMATION = 13
 
 
+def get_resampling_margin(sampling_frequency, margin_ms, decimation_factors=None):
+    """Return the margin in input samples, using an automatic estimate when margin_ms is None.
+
+    The estimate is an extension of the pole-decay heuristic illustrated in SciPy's filtfilt
+    documentation (theirs is for a single stage, extended here to the cascade).
+
+    Automatic margins are at least 100 ms and aligned to the total decimation factor.
+    FFT resampling margins remain 100 ms because this estimate only applies to IIR decimation.
+    An explicit margin_ms overrides the estimate.
+    """
+    if margin_ms is not None:
+        if not math.isfinite(margin_ms) or margin_ms < 0:
+            raise ValueError("margin_ms must be finite and nonnegative, or None")
+        return int(margin_ms * sampling_frequency / 1000)
+
+    margin = math.ceil(0.1 * sampling_frequency)
+    if decimation_factors is None:
+        return margin
+
+    from scipy.signal import cheby1
+
+    # Estimation method: for each default Chebyshev IIR stage, 
+    # estimate settling as ceil(log(1e-6) / log(r)), where r is the largest pole magnitude. 
+    # Convert each stage's estimate to input samples and sum them.
+    cascade_margin = 0
+    input_stride = 1
+    for factor in decimation_factors:
+        _, poles, _ = cheby1(8, 0.05, 0.8 / factor, output="zpk")
+        radius = max(abs(poles))
+        if not 0 < radius < 1:
+            raise ValueError("Cannot estimate a stable decimation margin. Specify margin_ms explicitly.")
+        cascade_margin += input_stride * math.ceil(math.log(1e-6) / math.log(radius))
+        input_stride *= factor
+
+    margin = max(margin, cascade_margin)
+    return ((margin + input_stride - 1) // input_stride) * input_stride
+
+
 def _prime_factors(n):
     """
     Return the prime factors of a positive integer `n` (ascending, with multiplicity).
diff --git a/src/spikeinterface/preprocessing/decimate.py b/src/spikeinterface/preprocessing/decimate.py
index 217b347022..2c3f5838b6 100644
--- a/src/spikeinterface/preprocessing/decimate.py
+++ b/src/spikeinterface/preprocessing/decimate.py
@@ -7,7 +7,7 @@
 
 from .basepreprocessor import BasePreprocessor
 from .filter import fix_dtype
-from ._decimation_tools import get_balanced_decimation_factors
+from ._decimation_tools import get_balanced_decimation_factors, get_resampling_margin
 from spikeinterface.core import BaseRecordingSegment, get_chunk_with_margin
 
 
@@ -43,10 +43,12 @@ class DecimateRecording(BasePreprocessor):
         is ignored. If omitted or None, currently behaves as False and emits a FutureWarning:
         a future release will enable antialiasing by default. Pass True or False explicitly
         to select the behavior and silence the transition warning.
-    margin_ms : float, default: 100.0
+    margin_ms : float | None, default: None
         Margin in ms used on each side of every chunk to limit edge effects of the anti-aliasing
         filter. Only used when `antialias=True`. The margin is internally rounded up to a whole
         number of output samples so the filtered, downsampled traces stay aligned across chunks.
+        If None, a suitable margin estimate based on the filter properties is used, with
+        a minimum of 100 ms. A nonnegative value overrides the estimate.
     dtype : dtype or None, default: None
         The dtype of the returned traces. If None, the dtype of the parent recording is used.
 
@@ -65,7 +67,7 @@ def __init__(
         decimation_factor,
         decimation_offset=0,
         antialias=None,
-        margin_ms=100.0,
+        margin_ms=None,
         dtype=None,
     ):
         # Original sampling frequency
@@ -103,7 +105,7 @@ def __init__(
         decimation_factors = get_balanced_decimation_factors(decimation_factor) if antialias else None
 
         # Margin (in parent samples) to limit anti-aliasing filter edge effects.
-        margin = int(margin_ms * self._orig_samp_freq / 1000)
+        margin = get_resampling_margin(self._orig_samp_freq, margin_ms, decimation_factors) if antialias else 0
 
         BasePreprocessor.__init__(self, recording, sampling_frequency=decimated_sampling_frequency, dtype=dtype)
 
diff --git a/src/spikeinterface/preprocessing/resample.py b/src/spikeinterface/preprocessing/resample.py
index cf8c7b7024..96a1a7246c 100644
--- a/src/spikeinterface/preprocessing/resample.py
+++ b/src/spikeinterface/preprocessing/resample.py
@@ -8,7 +8,7 @@
 
 from .basepreprocessor import BasePreprocessor
 from .filter import fix_dtype
-from ._decimation_tools import get_balanced_decimation_factors
+from ._decimation_tools import get_balanced_decimation_factors, get_resampling_margin
 from .decimate import get_antialiased_decimated_traces, _cast_resampled_traces
 from spikeinterface.core import get_chunk_with_margin, BaseRecordingSegment
 from spikeinterface.core.frameslicerecording import FrameSliceRecordingSegment
@@ -53,8 +53,11 @@ class ResampleRecording(BasePreprocessor):
         - 0.0: Strict mode — split on any gap >= 1.5 sample periods
         - 1.0: Tolerate gaps up to 1 ms, split on larger gaps
         - 100.0: Only major pauses (>100 ms) create sections
-    margin_ms : float, default: 100.0
-        Margin in ms for computations, used to decrease edge effects.
+    margin_ms : float | None, default: None
+        Margin in ms for computations, used to decrease edge effects. If None, integer-factor
+        decimation estimates the cascade's settling margin from its filter poles, with a
+        minimum of 100 ms. FFT resampling uses 100 ms. A nonnegative value overrides the
+        estimate. The estimate is a heuristic, not a bound on output error.
     dtype : dtype or None, default: None
         The dtype of the returned traces. If None, the dtype of the parent recording is used.
         Integer output is rounded and clipped to the dtype range. Nonfinite resampled
@@ -87,7 +90,7 @@ def __init__(
         recording,
         resample_rate,
         gap_tolerance_ms=None,
-        margin_ms=100.0,
+        margin_ms=None,
         dtype=None,
         skip_checks=False,
     ):
@@ -107,7 +110,7 @@ def __init__(
             assert check_nyquist(recording, resample_rate), "The requested resample rate would induce errors!"
 
         # Get a margin to avoid issues later
-        margin = int(margin_ms * recording.get_sampling_frequency() / 1000)
+        margin = get_resampling_margin(self._orig_samp_freq, margin_ms, decimation_factors)
 
         BasePreprocessor.__init__(self, recording, sampling_frequency=resample_rate, dtype=dtype)
         for parent_segment in recording.segments:
diff --git a/src/spikeinterface/preprocessing/tests/test_decimate.py b/src/spikeinterface/preprocessing/tests/test_decimate.py
index b17b616a31..1ab7dbb502 100644
--- a/src/spikeinterface/preprocessing/tests/test_decimate.py
+++ b/src/spikeinterface/preprocessing/tests/test_decimate.py
@@ -112,7 +112,7 @@ def test_balanced_decimation_factors(decimation_factor, max_factor, expected):
         assert all(f <= max_factor for f in factors)
 
 
-@pytest.mark.parametrize("decimation_factor", [6, 10, 48])
+@pytest.mark.parametrize("decimation_factor", [6, 10, 48, 300])
 def test_decimate_antialias_by_chunks(decimation_factor):
     # Mirror test_resample_by_chunks: chunked reads must match a full read once the
     # anti-aliasing margins are accounted for. Factor 48 exercises the internal multi-pass.
@@ -123,7 +123,7 @@ def test_decimate_antialias_by_chunks(decimation_factor):
     rms = np.sqrt(np.mean(parent_rec.get_traces() ** 2))
     decimated_rate = sampling_frequency / decimation_factor
 
-    for margin_ms in [100, 1000]:
+    for margin_ms in [None, 100, 1000]:
         rec2 = DecimateRecording(parent_rec, decimation_factor, antialias=True, margin_ms=margin_ms)
         chunk_size = int(decimated_rate * 2)  # ~2 seconds of the decimated signal
         rec3 = rec2.save(format="memory", chunk_size=chunk_size, n_jobs=1, progress_bar=False)

From 55fb9b8347c69b4fb5f07ab2a48eca8daff24f6e Mon Sep 17 00:00:00 2001
From: Graham Findlay <graham.b.findlay@gmail.com>
Date: Thu, 10 Sep 2026 13:53:59 -0500
Subject: [PATCH 08/11] Reject negative decimation offsets

---
 src/spikeinterface/preprocessing/decimate.py            | 4 ++--
 src/spikeinterface/preprocessing/tests/test_decimate.py | 8 ++++++++
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/spikeinterface/preprocessing/decimate.py b/src/spikeinterface/preprocessing/decimate.py
index 2c3f5838b6..32c0110272 100644
--- a/src/spikeinterface/preprocessing/decimate.py
+++ b/src/spikeinterface/preprocessing/decimate.py
@@ -75,8 +75,8 @@ def __init__(
         if not isinstance(decimation_factor, int) or decimation_factor <= 0:
             raise ValueError(f"Expecting strictly positive integer for `decimation_factor` arg")
         self._decimation_factor = decimation_factor
-        if not isinstance(decimation_offset, int) or decimation_factor < 0:
-            raise ValueError(f"Expecting positive integer for `decimation_factor` arg")
+        if not isinstance(decimation_offset, int) or decimation_offset < 0:
+            raise ValueError("Expecting a nonnegative integer for `decimation_offset` arg")
         parent_min_n_samp = min(
             [recording.get_num_samples(segment_index) for segment_index in range(recording.get_num_segments())]
         )
diff --git a/src/spikeinterface/preprocessing/tests/test_decimate.py b/src/spikeinterface/preprocessing/tests/test_decimate.py
index 1ab7dbb502..c0f535d61f 100644
--- a/src/spikeinterface/preprocessing/tests/test_decimate.py
+++ b/src/spikeinterface/preprocessing/tests/test_decimate.py
@@ -11,6 +11,14 @@
 import numpy as np
 
 
+def test_decimate_rejects_invalid_offset():
+    rec = NumpyRecording(np.zeros((24, 1), dtype="float32"), 1000)
+    for offset in [-1, -12, 0.5]:
+        for antialias in [False, True]:
+            with pytest.raises(ValueError, match="nonnegative integer.*decimation_offset"):
+                decimate(rec, 12, decimation_offset=offset, antialias=antialias)
+
+
 @pytest.mark.parametrize("num_segments", [1, 2])
 @pytest.mark.parametrize("decimation_offset", [0, 1, 5, 21, 101])
 @pytest.mark.parametrize("decimation_factor", [1, 7, 50])

From 7d1f0538fca09907473892fa07e63f6c9cfb8646 Mon Sep 17 00:00:00 2001
From: Graham Findlay <graham.b.findlay@gmail.com>
Date: Thu, 10 Sep 2026 13:59:59 -0500
Subject: [PATCH 09/11] Remove obsolete resampling Nyquist checks

---
 src/spikeinterface/preprocessing/resample.py | 38 +-------------------
 1 file changed, 1 insertion(+), 37 deletions(-)

diff --git a/src/spikeinterface/preprocessing/resample.py b/src/spikeinterface/preprocessing/resample.py
index 96a1a7246c..f2294f123d 100644
--- a/src/spikeinterface/preprocessing/resample.py
+++ b/src/spikeinterface/preprocessing/resample.py
@@ -1,10 +1,7 @@
 import numpy as np
 import warnings
 
-from spikeinterface.core.core_tools import (
-    define_function_handling_dict_from_class,
-    recursive_key_finder,
-)
+from spikeinterface.core.core_tools import define_function_handling_dict_from_class
 
 from .basepreprocessor import BasePreprocessor
 from .filter import fix_dtype
@@ -62,8 +59,6 @@ class ResampleRecording(BasePreprocessor):
         The dtype of the returned traces. If None, the dtype of the parent recording is used.
         Integer output is rounded and clipped to the dtype range. Nonfinite resampled
         output raises a ValueError before conversion.
-    skip_checks : bool, default: False
-        If True, checks on sampling frequencies and cutoff filter frequencies are skipped
 
     Returns
     -------
@@ -92,7 +87,6 @@ def __init__(
         gap_tolerance_ms=None,
         margin_ms=None,
         dtype=None,
-        skip_checks=False,
     ):
         self._orig_samp_freq = recording.get_sampling_frequency()
         self._resample_rate = resample_rate
@@ -105,9 +99,6 @@ def __init__(
             decimation_factors = None
         # fix_dtype not always returns the str, make sure it does
         dtype = fix_dtype(recording, dtype).str
-        # Ensure that the requested resample rate is doable:
-        if skip_checks:
-            assert check_nyquist(recording, resample_rate), "The requested resample rate would induce errors!"
 
         # Get a margin to avoid issues later
         margin = get_resampling_margin(self._orig_samp_freq, margin_ms, decimation_factors)
@@ -132,7 +123,6 @@ def __init__(
             gap_tolerance_ms=gap_tolerance_ms,
             margin_ms=margin_ms,
             dtype=dtype,
-            skip_checks=skip_checks,
         )
 
 
@@ -352,29 +342,3 @@ def _get_traces_gapped(self, start_frame, end_frame, channel_indices):
 
 
 resample = define_function_handling_dict_from_class(source_class=ResampleRecording, name="resample")
-
-
-# Some helpers to do checks
-def check_nyquist(recording, resample_rate):
-    # Check that the original and requested sampling rates will not induce aliasing
-    # Basic test, compare the sampling frequency with the resample rate
-    sampling_frequency_check = recording.get_sampling_frequency() / 2 > resample_rate
-    # Check that the signal, if it has been filtered, is still not violating
-    if recording.is_filtered():
-        # Check if we have access to the highcut frequency
-        freq_max = list(recursive_key_finder(recording, "freq_max"))
-        if freq_max:
-            # Given that there might be more than one filter applied, keep the lowest
-            freq_max = min(freq_max)
-            lowpass_cutoff_check = freq_max / 2 > resample_rate
-        else:
-            # If has been filterd but unknown high cutoff, give warning and asume the best
-            warnings.warn("The recording is filtered, but we can't ensure that it complies with the Nyquist limit.")
-            lowpass_cutoff_check = True
-    else:
-        # If it hasn't been filtered, we only depend on the previous test
-        warnings.warn(
-            "The recording is not filtered, so cutoff frequencies cannot be checked. " "Use resampling with caution"
-        )
-        lowpass_cutoff_check = True
-    return all([sampling_frequency_check, lowpass_cutoff_check])

From 9fd05d7c7567fa2429d89c2a600b54ffda4719e0 Mon Sep 17 00:00:00 2001
From: Graham Findlay <graham.b.findlay@gmail.com>
Date: Thu, 10 Sep 2026 14:20:16 -0500
Subject: [PATCH 10/11] Use polyphase filtering consistently for resampling and
 decimation

Share scipy.signal.resample_poly between resample() and antialiased
decimate(), replacing the IIR / FFT approaches. Derive chunk
margins from FIR support.

Given a target sample rate, select the closest rational ratio under
max_denominator, report the achieved sampling rate, and warn when it
differs from the request.
---
 .../preprocessing/_decimation_tools.py        | 145 -------------
 .../preprocessing/_resampling_tools.py        |  58 +++++
 src/spikeinterface/preprocessing/decimate.py  | 112 ++++------
 src/spikeinterface/preprocessing/resample.py  | 199 +++++++-----------
 .../preprocessing/tests/test_decimate.py      |  83 ++------
 .../preprocessing/tests/test_resample.py      | 177 ++++++++--------
 6 files changed, 275 insertions(+), 499 deletions(-)
 delete mode 100644 src/spikeinterface/preprocessing/_decimation_tools.py
 create mode 100644 src/spikeinterface/preprocessing/_resampling_tools.py

diff --git a/src/spikeinterface/preprocessing/_decimation_tools.py b/src/spikeinterface/preprocessing/_decimation_tools.py
deleted file mode 100644
index e6b4aed6e6..0000000000
--- a/src/spikeinterface/preprocessing/_decimation_tools.py
+++ /dev/null
@@ -1,145 +0,0 @@
-"""
-Helpers for splitting a (potentially large) integer decimation factor into several balanced
-sub-factors, so that anti-aliased decimation can be applied as multiple stable scipy.signal.decimate
-passes. Shared by DecimateRecording and ResampleRecording.
-"""
-
-import math
-import warnings
-
-# scipy.signal.decimate uses an order-8 Chebyshev type I IIR filter by default, and its
-# documentation recommends decimating in several steps rather than a single step
-# for downsampling factors larger than this value.
-_MAX_SINGLE_PASS_DECIMATION = 13
-
-
-def get_resampling_margin(sampling_frequency, margin_ms, decimation_factors=None):
-    """Return the margin in input samples, using an automatic estimate when margin_ms is None.
-
-    The estimate is an extension of the pole-decay heuristic illustrated in SciPy's filtfilt
-    documentation (theirs is for a single stage, extended here to the cascade).
-
-    Automatic margins are at least 100 ms and aligned to the total decimation factor.
-    FFT resampling margins remain 100 ms because this estimate only applies to IIR decimation.
-    An explicit margin_ms overrides the estimate.
-    """
-    if margin_ms is not None:
-        if not math.isfinite(margin_ms) or margin_ms < 0:
-            raise ValueError("margin_ms must be finite and nonnegative, or None")
-        return int(margin_ms * sampling_frequency / 1000)
-
-    margin = math.ceil(0.1 * sampling_frequency)
-    if decimation_factors is None:
-        return margin
-
-    from scipy.signal import cheby1
-
-    # Estimation method: for each default Chebyshev IIR stage, 
-    # estimate settling as ceil(log(1e-6) / log(r)), where r is the largest pole magnitude. 
-    # Convert each stage's estimate to input samples and sum them.
-    cascade_margin = 0
-    input_stride = 1
-    for factor in decimation_factors:
-        _, poles, _ = cheby1(8, 0.05, 0.8 / factor, output="zpk")
-        radius = max(abs(poles))
-        if not 0 < radius < 1:
-            raise ValueError("Cannot estimate a stable decimation margin. Specify margin_ms explicitly.")
-        cascade_margin += input_stride * math.ceil(math.log(1e-6) / math.log(radius))
-        input_stride *= factor
-
-    margin = max(margin, cascade_margin)
-    return ((margin + input_stride - 1) // input_stride) * input_stride
-
-
-def _prime_factors(n):
-    """
-    Return the prime factors of a positive integer `n` (ascending, with multiplicity).
-
-    Examples
-    --------
-    >>> _prime_factors(60)
-    [2, 2, 3, 5]
-    >>> _prime_factors(17)
-    [17]
-    """
-    factors = []
-    divisor = 2
-    while divisor * divisor <= n:
-        while n % divisor == 0:
-            factors.append(divisor)
-            n //= divisor
-        divisor += 1
-    if n > 1:
-        factors.append(n)
-    return factors
-
-
-def _greedy_pack(primes_desc, num_bins, max_factor=_MAX_SINGLE_PASS_DECIMATION):
-    """
-    Greedily pack `primes_desc` (largest first) into `num_bins` bins, keeping each bin's
-    product <= `max_factor` and the bins as balanced as possible.
-
-    Returns the list of bin products, or None if some prime cannot be placed (i.e. `num_bins`
-    is too small to keep every bin <= the single-pass limit).
-
-    Examples
-    --------
-    Pack the prime factors of 48 into two balanced bins (6 and 8):
-
-    >>> _greedy_pack([3, 2, 2, 2, 2], 2)
-    [6, 8]
-
-    Two bins cannot hold 2 ** 7 = 128 without a bin exceeding the single-pass limit of 13:
-
-    >>> _greedy_pack([2, 2, 2, 2, 2, 2, 2], 2) is None
-    True
-    """
-    bins = [1] * num_bins
-    for prime in primes_desc:
-        fitting = [i for i in range(num_bins) if bins[i] * prime <= max_factor]
-        if not fitting:
-            return None
-        # Place into the smallest fitting bin (ties broken by index, for determinism).
-        target = min(fitting, key=lambda i: (bins[i], i))
-        bins[target] *= prime
-    return bins
-
-
-def get_balanced_decimation_factors(decimation_factor, max_factor=_MAX_SINGLE_PASS_DECIMATION):
-    """
-    Split `decimation_factor` into balanced sub-factors no greater than `max_factor`.
-
-    SciPy recommends multiple IIR decimation passes for factors above 13.
-    Balancing the factors (e.g. 48 -> [8, 6] rather than [12, 4]) further aids stability.
-    The product of the returned factors always equals `decimation_factor`.
-
-    If `decimation_factor` has a prime factor greater than 13 (e.g. a large prime such as 17),
-    no valid split exists, a warning is issued, and `[decimation_factor]` is returned;
-    it is the caller's responsibility to handle this (e.g., warn that a single,
-    potentially unstable, pass will be used).
-    """
-    if not isinstance(max_factor, int) or max_factor < 2:
-        raise ValueError("max_factor must be an integer greater than one")
-    if decimation_factor <= max_factor:
-        return [decimation_factor]
-
-    primes = _prime_factors(decimation_factor)
-    if max(primes) > max_factor:
-        warnings.warn(
-            f"`decimation_factor`={decimation_factor} cannot be split into anti-aliasing passes of <= {max_factor} "
-            f"(it has a prime factor > {max_factor}). A single `scipy.signal.decimate` pass will be used, "
-            f"which may be unstable. Consider a `decimation_factor` without large prime factors.",
-            stacklevel=2,
-        )
-        return [decimation_factor]
-
-    primes_desc = sorted(primes, reverse=True)
-    # Minimum number of passes so that, ideally, each pass decimates by <= max_factor.
-    num_passes = max(1, math.ceil(math.log(decimation_factor) / math.log(max_factor)))
-    while num_passes <= len(primes_desc):
-        bins = _greedy_pack(primes_desc, num_passes, max_factor)
-        if bins is not None:
-            return sorted(bins, reverse=True)
-        num_passes += 1
-    # Fallback: one prime per pass (always valid since every prime is <= max_factor).
-    return primes_desc
diff --git a/src/spikeinterface/preprocessing/_resampling_tools.py b/src/spikeinterface/preprocessing/_resampling_tools.py
new file mode 100644
index 0000000000..72d516a7f8
--- /dev/null
+++ b/src/spikeinterface/preprocessing/_resampling_tools.py
@@ -0,0 +1,58 @@
+"""Rational rate selection and FIR design for polyphase resampling."""
+
+import math
+import warnings
+from fractions import Fraction
+
+import numpy as np
+
+
+def get_resampling_factors(parent_rate, resample_rate, max_denominator):
+    """Select the closest ratio under the denominator limit and warn about rate differences."""
+    if not math.isfinite(parent_rate) or parent_rate <= 0:
+        raise ValueError("The parent sampling frequency must be finite and positive")
+    if not math.isfinite(resample_rate) or resample_rate <= 0:
+        raise ValueError("resample_rate must be finite and positive")
+    if not isinstance(max_denominator, int) or max_denominator < 1:
+        raise ValueError("max_denominator must be a positive integer")
+
+    ratio = (Fraction(float(resample_rate)) / Fraction(float(parent_rate))).limit_denominator(max_denominator)
+    up, down = ratio.numerator, ratio.denominator
+    if up == 0:
+        raise ValueError("The requested rate is too low for max_denominator; increase max_denominator")
+    achieved_rate = float(Fraction(float(parent_rate)) * ratio)
+    if not math.isclose(achieved_rate, resample_rate, rel_tol=1e-12, abs_tol=0):
+        error_ppm = (achieved_rate / resample_rate - 1) * 1e6
+        warnings.warn(
+            f"Requested resample_rate={resample_rate:.16g} Hz; the polyphase ratio {up}/{down} "
+            f"achieves {achieved_rate:.16g} Hz ({error_ppm:+.6g} ppm). The output sampling frequency "
+            "uses the achieved rate. Increase max_denominator for a closer approximation.",
+            stacklevel=2,
+        )
+    return up, down, achieved_rate
+
+
+def get_polyphase_filter(sampling_frequency, up, down, margin_ms):
+    """Design SciPy's default FIR and cover its support on an aligned input grid."""
+    from scipy.signal import firwin
+
+    if margin_ms is not None and (not math.isfinite(margin_ms) or margin_ms < 0):
+        raise ValueError("margin_ms must be finite and nonnegative, or None")
+
+    if up == down == 1:
+        return np.ones(1), 0
+
+    # Important! The multiplier 10 and Kaiser parameter 5.0 come directly from SciPy’s
+    # default resample_poly design. Don't change them!
+    half_length = 10 * max(up, down)
+    coefficients = firwin(
+        2 * half_length + 1,  # odd length gives a symmetric filter with a central sample
+        1.0 / max(up, down),
+        window=("kaiser", 5.0),
+    )
+
+    margin = (half_length + up - 1) // up  # Convert filter support to input samples
+    if margin_ms is not None:
+        margin = max(margin, math.ceil(margin_ms * sampling_frequency / 1000))
+    margin = ((margin + down - 1) // down) * down
+    return coefficients, margin
diff --git a/src/spikeinterface/preprocessing/decimate.py b/src/spikeinterface/preprocessing/decimate.py
index 32c0110272..358dbce460 100644
--- a/src/spikeinterface/preprocessing/decimate.py
+++ b/src/spikeinterface/preprocessing/decimate.py
@@ -7,7 +7,7 @@
 
 from .basepreprocessor import BasePreprocessor
 from .filter import fix_dtype
-from ._decimation_tools import get_balanced_decimation_factors, get_resampling_margin
+from ._resampling_tools import get_polyphase_filter
 from spikeinterface.core import BaseRecordingSegment, get_chunk_with_margin
 
 
@@ -18,7 +18,7 @@ class DecimateRecording(BasePreprocessor):
     By default this uses simple array slicing
     (``<parent_traces>[<decimation_offset>::<decimation_factor>]``), which is fast but applies no
     anti-aliasing filter and so might introduce aliasing, or skip across signal of interest. Set
-    `antialias=True` to low-pass filter before downsampling using ``scipy.signal.decimate`` (the
+    `antialias=True` to low-pass filter before downsampling using ``scipy.signal.resample_poly`` (the
     same anti-aliased decimation used by ``spikeinterface.preprocessing.ResampleRecording``).
 
     Parameters
@@ -36,19 +36,15 @@ class DecimateRecording(BasePreprocessor):
         The same decimation offset is applied to all segments from the parent recording.
     antialias : bool | None, default: None
         If True, apply an anti-aliasing low-pass filter before downsampling, using
-        ``scipy.signal.decimate``. When `decimation_factor` exceeds 13, the decimation is
-        automatically performed in several balanced sub-13 passes (e.g. a factor of 48 is applied
-        as 8 then 6), as scipy recommends, to keep the IIR anti-aliasing filter stable. If False,
+        ``scipy.signal.resample_poly`` with a Kaiser-windowed FIR filter. If False,
         traces are downsampled by plain array slicing with no filtering, and `margin_ms`
         is ignored. If omitted or None, currently behaves as False and emits a FutureWarning:
         a future release will enable antialiasing by default. Pass True or False explicitly
         to select the behavior and silence the transition warning.
     margin_ms : float | None, default: None
-        Margin in ms used on each side of every chunk to limit edge effects of the anti-aliasing
-        filter. Only used when `antialias=True`. The margin is internally rounded up to a whole
-        number of output samples so the filtered, downsampled traces stay aligned across chunks.
-        If None, a suitable margin estimate based on the filter properties is used, with
-        a minimum of 100 ms. A nonnegative value overrides the estimate.
+        Additional context in ms on each side of a chunk. Only used when `antialias=True`.
+        If None, use the FIR filter's finite support. An explicit nonnegative value requests
+        at least that much context; filter support and sample-grid alignment are always retained.
     dtype : dtype or None, default: None
         The dtype of the returned traces. If None, the dtype of the parent recording is used.
 
@@ -102,10 +98,10 @@ def __init__(
             )
             antialias = False
 
-        decimation_factors = get_balanced_decimation_factors(decimation_factor) if antialias else None
-
-        # Margin (in parent samples) to limit anti-aliasing filter edge effects.
-        margin = get_resampling_margin(self._orig_samp_freq, margin_ms, decimation_factors) if antialias else 0
+        if antialias:
+            filter_coefficients, margin = get_polyphase_filter(self._orig_samp_freq, 1, decimation_factor, margin_ms)
+        else:
+            filter_coefficients, margin = None, 0
 
         BasePreprocessor.__init__(self, recording, sampling_frequency=decimated_sampling_frequency, dtype=dtype)
 
@@ -120,7 +116,7 @@ def __init__(
                     self._dtype,
                     antialias,
                     margin,
-                    decimation_factors,
+                    filter_coefficients,
                 )
             )
 
@@ -145,7 +141,7 @@ def __init__(
         dtype,
         antialias=False,
         margin=0,
-        decimation_factors=None,
+        filter_coefficients=None,
     ):
         if parent_recording_segment._time_vector is not None:
             time_vector = parent_recording_segment._time_vector[decimation_offset::decimation_factor]
@@ -153,10 +149,9 @@ def __init__(
             t_start = None
         else:
             time_vector = None
-            if parent_recording_segment._t_start is None:
-                t_start = None
-            else:
-                t_start = parent_recording_segment._t_start + (decimation_offset / parent_rate)
+            t_start = parent_recording_segment._t_start
+            if decimation_offset:
+                t_start = (0.0 if t_start is None else t_start) + decimation_offset / parent_rate
 
         # Do not use BasePreprocessorSegment bcause we have to reset the sampling rate!
         BaseRecordingSegment.__init__(
@@ -168,12 +163,12 @@ def __init__(
         self._dtype = dtype
         self._antialias = antialias
         self._margin = margin
-        self._decimation_factors = decimation_factors if decimation_factors is not None else [decimation_factor]
+        self._filter_coefficients = filter_coefficients
 
     def get_num_samples(self):
         parent_n_samp = self._parent_segment.get_num_samples()
         assert self._decimation_offset < parent_n_samp  # Sanity check (already enforced). Formula changes otherwise
-        return int(np.ceil((parent_n_samp - self._decimation_offset) / self._decimation_factor))
+        return (parent_n_samp - self._decimation_offset + self._decimation_factor - 1) // self._decimation_factor
 
     def get_traces(self, start_frame, end_frame, channel_indices):
         if not self._antialias:
@@ -188,69 +183,41 @@ def get_traces(self, start_frame, end_frame, channel_indices):
                 :: self._decimation_factor
             ].astype(self._dtype)
 
-        # Anti-aliased decimation as a cascade of balanced scipy.signal.decimate passes.
-        return get_antialiased_decimated_traces(
+        return get_polyphase_resampled_traces(
             self._parent_segment,
             start_frame,
             end_frame,
             channel_indices,
+            1,
             self._decimation_factor,
-            self._decimation_factors,
             self._margin,
             self._dtype,
+            self._filter_coefficients,
             decimation_offset=self._decimation_offset,
         )
 
 
-def get_antialiased_decimated_traces(
+def get_polyphase_resampled_traces(
     parent_segment,
     start_frame,
     end_frame,
     channel_indices,
-    decimation_factor,
-    decimation_factors,
+    up,
+    down,
     margin,
     dtype,
+    filter_coefficients,
     decimation_offset=0,
 ):
-    """
-    Fetch a margined chunk from `parent_segment` and decimate it by `decimation_factor`, applied
-    as a cascade of the balanced `decimation_factors` passes of ``scipy.signal.decimate``.
+    """Resample a chunk, with reflected boundary padding."""
+    from scipy.signal import resample_poly
 
-    The margin is rounded up to a multiple of the total `decimation_factor` so that
-    ``left_margin // decimation_factor`` is exact; combined with scipy's default
-    ``zero_phase=True`` (output sample i maps to filtered input sample i * factor), this keeps the
-    downsampled traces aligned across chunks (a chunked read matches a full read). Exactly
-    ``end_frame - start_frame`` decimated samples are returned.
+    if end_frame <= start_frame:
+        return parent_segment.get_traces(0, 0, channel_indices).astype(dtype)
 
-    Parameters
-    ----------
-    parent_segment : BaseRecordingSegment
-        The parent segment to read (full-rate) traces from.
-    start_frame, end_frame : int
-        Output (decimated) frame range to return.
-    channel_indices : slice | list | np.ndarray | None
-        Channels to read, forwarded to the parent segment.
-    decimation_factor : int
-        The total decimation factor (the product of `decimation_factors`).
-    decimation_factors : list[int]
-        The per-pass sub-factors (each <= 13), e.g. from `get_balanced_decimation_factors`.
-    margin : int
-        Margin in parent samples used to limit anti-aliasing filter edge effects. Rounded up
-        internally to a multiple of `decimation_factor`.
-    dtype : np.dtype | str
-        Output dtype. Integer output is rounded and clipped to its range.
-    decimation_offset : int, default: 0
-        Index of the first parent frame, applied to the first output sample only.
-    """
-    from scipy import signal
-
-    q = decimation_factor
-    parent_start_frame = decimation_offset + start_frame * q
-    parent_end_frame = parent_start_frame + (end_frame - start_frame) * q
-    # Round the margin up to a multiple of q so that left_margin // q is exact.
-    margin = ((margin + q - 1) // q) * q
-    parent_traces, left_margin, right_margin = get_chunk_with_margin(
+    parent_start_frame = decimation_offset + (start_frame // up) * down
+    parent_end_frame = decimation_offset + ((end_frame + up - 1) // up) * down
+    parent_traces, left_margin, _ = get_chunk_with_margin(
         parent_segment,
         parent_start_frame,
         parent_end_frame,
@@ -259,13 +226,16 @@ def get_antialiased_decimated_traces(
         add_reflect_padding=True,
     )
     working_dtype = np.result_type(parent_traces.dtype, dtype, np.float32)
-    decimated_traces = parent_traces.astype(working_dtype, copy=False)
-    for sub_q in decimation_factors:
-        decimated_traces = signal.decimate(decimated_traces, q=sub_q, axis=0)
-    start_drop = left_margin // q
-    n_out = end_frame - start_frame
-    decimated_traces = decimated_traces[start_drop : start_drop + n_out]
-    return _cast_resampled_traces(decimated_traces, dtype)
+    traces = resample_poly(
+        parent_traces.astype(working_dtype, copy=False),
+        up,
+        down,
+        axis=0,
+        window=filter_coefficients.astype(working_dtype, copy=False),
+    )
+    start_drop = start_frame % up + left_margin * up // down
+    traces = traces[start_drop : start_drop + end_frame - start_frame]
+    return _cast_resampled_traces(traces, dtype)
 
 
 def _cast_resampled_traces(traces, dtype):
diff --git a/src/spikeinterface/preprocessing/resample.py b/src/spikeinterface/preprocessing/resample.py
index f2294f123d..c4c9b11abd 100644
--- a/src/spikeinterface/preprocessing/resample.py
+++ b/src/spikeinterface/preprocessing/resample.py
@@ -5,9 +5,9 @@
 
 from .basepreprocessor import BasePreprocessor
 from .filter import fix_dtype
-from ._decimation_tools import get_balanced_decimation_factors, get_resampling_margin
-from .decimate import get_antialiased_decimated_traces, _cast_resampled_traces
-from spikeinterface.core import get_chunk_with_margin, BaseRecordingSegment
+from ._resampling_tools import get_resampling_factors, get_polyphase_filter
+from .decimate import get_polyphase_resampled_traces
+from spikeinterface.core import BaseRecordingSegment
 from spikeinterface.core.frameslicerecording import FrameSliceRecordingSegment
 
 
@@ -15,19 +15,18 @@ class ResampleRecording(BasePreprocessor):
     """
     Resample the recording extractor traces.
 
-    If the parent sampling rate is an exact integer multiple of `resample_rate`, the
-    ``signal.decimate`` method from scipy is used (anti-aliased decimation). In other cases
-    ``signal.resample`` is used, in which case the resulting signal can have issues on the edges,
-    mainly on the rightmost. See Notes for a caveat on how the integer multiple is detected.
+    Uses ``scipy.signal.resample_poly`` with a Kaiser-windowed FIR filter for both
+    downsampling and upsampling. Detected gaps are handled section by section.
 
     Parameters
     ----------
     recording : Recording
         The recording extractor to be re-referenced
     resample_rate : int | float
-        The resampling frequency. Integer ratios (parent_rate / resample_rate) use
-        ``scipy.signal.decimate``; non-integer ratios use ``scipy.signal.resample`` (FFT-based),
-        which can have edge effects, mainly on the rightmost samples.
+        The requested sampling frequency. The closest output/input rate ratio with 
+        denominator at most `max_denominator` is selected. The output reports the 
+        achieved rate, while the requested rate is retained in serialization kwargs.
+        A relative difference exceeding 1e-12 emits a warning. 
     gap_tolerance_ms : float | None, default: None
         Maximum acceptable gap size in milliseconds for automatic segmentation.
 
@@ -51,15 +50,19 @@ class ResampleRecording(BasePreprocessor):
         - 1.0: Tolerate gaps up to 1 ms, split on larger gaps
         - 100.0: Only major pauses (>100 ms) create sections
     margin_ms : float | None, default: None
-        Margin in ms for computations, used to decrease edge effects. If None, integer-factor
-        decimation estimates the cascade's settling margin from its filter poles, with a
-        minimum of 100 ms. FFT resampling uses 100 ms. A nonnegative value overrides the
-        estimate. The estimate is a heuristic, not a bound on output error.
+        Additional context in ms on each side of a chunk. If None, use the FIR filter's
+        finite support. An explicit nonnegative value requests at least that much context;
+        filter support is always retained within each section.
     dtype : dtype or None, default: None
         The dtype of the returned traces. If None, the dtype of the parent recording is used.
         Integer output is rounded and clipped to the dtype range. Nonfinite resampled
         output raises a ValueError before conversion.
 
+    max_denominator : int, default: 10000
+        Maximum denominator of the rational output/input rate ratio. Increasing this can
+        improve rate accuracy, but can also increase filter length and the aligned input
+        span needed for a chunk. 
+
     Returns
     -------
     resample_recording : ResampleRecording
@@ -67,16 +70,15 @@ class ResampleRecording(BasePreprocessor):
 
     Notes
     -----
-    The (anti-aliased) decimation path is selected by an exact check,
-    ``parent_rate % resample_rate == 0``. This only detects an integer downsampling factor when
-    both rates make that modulo exactly zero. If either the parent rate or `resample_rate` is a
-    non-integer float (i.e. ``float(int(x)) != float(x)``), a conceptually integer ratio can go
-    undetected and silently fall back to the FFT-based ``scipy.signal.resample`` path. For example,
-    decimating a 625 Hz recording by a factor of 6 means a target of 104.1666... Hz, and
-    ``625 % 104.1666... != 0``, so the integer-decimation path is not taken (whereas a factor of 5,
-    i.e. a 125 Hz target, is detected since ``625 % 125 == 0``). To force anti-aliased integer
-    decimation by a known factor regardless of the rates, use
-    ``spikeinterface.preprocessing.DecimateRecording`` with ``antialias=True``.
+    Each section returns ``ceil(num_input_samples * up / down)`` samples, matching SciPy
+    and ``decimate()``. Output timestamps use the same rational grid as the traces. 
+    Explicit parent timestamps are sampled or interpolated within each section. 
+    Output positions beyond the last input sample extrapolate its timestamp by less 
+    than one nominal input period.
+
+    For example, resampling from 30000.01 Hz to a requested 2500 Hz with the default 
+    denominator limit selects up=1 and down=12. The output reports 2500.000833333333 Hz, 
+    and a warning reports the difference of approximately +0.333333 ppm.
 
     """
 
@@ -87,33 +89,29 @@ def __init__(
         gap_tolerance_ms=None,
         margin_ms=None,
         dtype=None,
+        max_denominator=10000,
     ):
         self._orig_samp_freq = recording.get_sampling_frequency()
-        self._resample_rate = resample_rate
-        self._sampling_frequency = resample_rate
-        # Exact integer-factor downsampling uses one or more antialiased decimation passes.
-        if self._orig_samp_freq % resample_rate == 0:
-            decimation_factor = int(self._orig_samp_freq / resample_rate)
-            decimation_factors = get_balanced_decimation_factors(decimation_factor)
-        else:
-            decimation_factors = None
+        up, down, achieved_rate = get_resampling_factors(self._orig_samp_freq, resample_rate, max_denominator)
+        self._resample_rate = achieved_rate
         # fix_dtype not always returns the str, make sure it does
         dtype = fix_dtype(recording, dtype).str
 
-        # Get a margin to avoid issues later
-        margin = get_resampling_margin(self._orig_samp_freq, margin_ms, decimation_factors)
+        filter_coefficients, margin = get_polyphase_filter(self._orig_samp_freq, up, down, margin_ms)
 
-        BasePreprocessor.__init__(self, recording, sampling_frequency=resample_rate, dtype=dtype)
+        BasePreprocessor.__init__(self, recording, sampling_frequency=achieved_rate, dtype=dtype)
         for parent_segment in recording.segments:
             self.add_recording_segment(
                 ResampleRecordingSegment(
                     parent_segment,
-                    resample_rate,
+                    achieved_rate,
                     recording.get_sampling_frequency(),
                     margin,
                     dtype,
                     gap_tolerance_ms,
-                    decimation_factors,
+                    up,
+                    down,
+                    filter_coefficients,
                 )
             )
 
@@ -123,6 +121,7 @@ def __init__(
             gap_tolerance_ms=gap_tolerance_ms,
             margin_ms=margin_ms,
             dtype=dtype,
+            max_denominator=max_denominator,
         )
 
 
@@ -134,8 +133,10 @@ def __init__(
         parent_rate,
         margin,
         dtype,
-        gap_tolerance_ms=None,
-        decimation_factors=None,
+        gap_tolerance_ms,
+        up,
+        down,
+        filter_coefficients,
     ):
         self._resample_rate = resample_rate
         self._parent_segment = parent_recording_segment
@@ -143,9 +144,9 @@ def __init__(
         self._margin = margin
         self._dtype = dtype
         self._has_gaps = False
-        # Per-pass integer decimation factors when the ratio is an exact integer, else None
-        # (non-integer ratio -> FFT-based scipy.signal.resample).
-        self._decimation_factors = decimation_factors
+        self._up = up
+        self._down = down
+        self._filter_coefficients = filter_coefficients
 
         # Compute time_vector or t_start, following the pattern from DecimateRecordingSegment.
         # Do not use BasePreprocessorSegment because we have to reset the sampling rate!
@@ -196,7 +197,7 @@ def __init__(
             K = len(sec_boundaries_parent)
             sec_n_out = np.array(
                 [
-                    int((sec_boundaries_parent[k, 1] - sec_boundaries_parent[k, 0]) / parent_rate * resample_rate)
+                    (int(sec_boundaries_parent[k, 1] - sec_boundaries_parent[k, 0]) * up + down - 1) // down
                     for k in range(K)
                 ],
                 dtype=np.int64,
@@ -209,49 +210,11 @@ def __init__(
             self._sec_boundaries_output = sec_boundaries_output
             self._sec_n_out = sec_n_out
 
-            # Compute time_vector
-            n_out = int(len(parent_tv) / parent_rate * resample_rate)
-
-            if parent_rate % resample_rate == 0:
-                q_int = int(parent_rate / resample_rate)
-                if not self._has_gaps:
-                    time_vector = parent_tv[::q_int][:n_out]
-                else:
-                    # Section-wise slicing to keep time_vector consistent
-                    # with _sec_boundaries_output
-                    tv_pieces = []
-                    for k in range(K):
-                        p_start, p_end = sec_boundaries_parent[k]
-                        n_out_k = sec_n_out[k]
-                        if n_out_k == 0:
-                            continue
-                        tv_pieces.append(parent_tv[p_start:p_end:q_int][:n_out_k])
-                    time_vector = np.concatenate(tv_pieces)
-            elif not self._has_gaps:
-                # Non-integer ratio, no gaps: existing fast path
-                warnings.warn(
-                    "Resampling with a non-integer ratio requires interpolating the time_vector. "
-                    "An integer ratio (parent_rate / resample_rate) is more performant."
-                )
-                parent_indices = np.linspace(0, len(parent_tv) - 1, n_out)
-                time_vector = np.interp(parent_indices, np.arange(len(parent_tv)), parent_tv)
-            else:
-                # Non-integer ratio with gaps: per-section interpolation
-                warnings.warn(
-                    "Resampling with a non-integer ratio requires interpolating the time_vector. "
-                    "An integer ratio (parent_rate / resample_rate) is more performant."
-                )
-                tv_pieces = []
-                for k in range(K):
-                    p_start, p_end = sec_boundaries_parent[k]
-                    n_out_k = sec_n_out[k]
-                    if n_out_k == 0:
-                        continue
-                    sec_parent_tv = parent_tv[p_start:p_end]
-                    sec_len = p_end - p_start
-                    sec_indices = np.linspace(0, sec_len - 1, n_out_k)
-                    tv_pieces.append(np.interp(sec_indices, np.arange(sec_len), sec_parent_tv))
-                time_vector = np.concatenate(tv_pieces)
+            tv_pieces = [
+                _resample_time_vector(parent_tv[p_start:p_end], up, down, parent_rate)
+                for p_start, p_end in sec_boundaries_parent
+            ]
+            time_vector = np.concatenate(tv_pieces) if self._has_gaps else tv_pieces[0]
 
             BaseRecordingSegment.__init__(self, sampling_frequency=None, t_start=None, time_vector=time_vector)
         else:
@@ -262,57 +225,29 @@ def __init__(
     def get_num_samples(self):
         if self._time_vector is not None:
             return len(self._time_vector)
-        return int(self._parent_segment.get_num_samples() / self._parent_rate * self._resample_rate)
+        n = self._parent_segment.get_num_samples()
+        return (n * self._up + self._down - 1) // self._down
 
     def get_traces(self, start_frame, end_frame, channel_indices):
+        if end_frame <= start_frame:
+            return self._parent_segment.get_traces(0, 0, channel_indices).astype(self._dtype)
         if self._has_gaps:
             return self._get_traces_gapped(start_frame, end_frame, channel_indices)
 
         return self._get_resampled_traces(self._parent_segment, start_frame, end_frame, channel_indices)
 
     def _get_resampled_traces(self, parent_segment, start_frame, end_frame, channel_indices):
-        if self._decimation_factors is not None:
-            # Integer-factor downsampling, no gaps: one or more antialiased decimation passes.
-            decimation_factor = int(self._parent_rate / self._resample_rate)
-            return get_antialiased_decimated_traces(
-                parent_segment,
-                start_frame,
-                end_frame,
-                channel_indices,
-                decimation_factor,
-                self._decimation_factors,
-                self._margin,
-                self._dtype,
-            )
-
-        # Non-integer-ratio downsampling or upsampling, no gaps: FFT with proportional margins.
-        from scipy import signal
-
-        parent_start_frame, parent_end_frame = [
-            int((frame / self._resample_rate) * self._parent_rate) for frame in [start_frame, end_frame]
-        ]
-        parent_traces, left_margin, right_margin = get_chunk_with_margin(
+        return get_polyphase_resampled_traces(
             parent_segment,
-            parent_start_frame,
-            parent_end_frame,
+            start_frame,
+            end_frame,
             channel_indices,
+            self._up,
+            self._down,
             self._margin,
-            add_reflect_padding=True,
+            self._dtype,
+            self._filter_coefficients,
         )
-        working_dtype = np.result_type(parent_traces.dtype, self._dtype, np.float32)
-        parent_traces = parent_traces.astype(working_dtype, copy=False)
-        # get left and right margins for the resampled case
-        left_margin_rs, right_margin_rs = [
-            int((margin / self._parent_rate) * self._resample_rate) for margin in [left_margin, right_margin]
-        ]
-
-        # get the size for the resampled traces
-        num = int((end_frame + right_margin_rs) - (start_frame - left_margin_rs))
-        resampled_traces = signal.resample(parent_traces, num, axis=0)
-
-        # now take care of the edges
-        resampled_traces = resampled_traces[left_margin_rs : num - right_margin_rs]
-        return _cast_resampled_traces(resampled_traces, self._dtype)
 
     def _get_traces_gapped(self, start_frame, end_frame, channel_indices):
         """Resample each section with margins bounded by its own samples."""
@@ -341,4 +276,18 @@ def _get_traces_gapped(self, start_frame, end_frame, channel_indices):
         return result
 
 
+def _resample_time_vector(parent_times, up, down, parent_rate):
+    """Map the rational sample grid onto timestamps without crossing section boundaries."""
+    if up == 1:
+        return parent_times[::down]
+    n_out = (len(parent_times) * up + down - 1) // down
+    positions = np.arange(n_out, dtype=np.int64) * down
+    left = positions // up
+    weight = (positions % up) / up
+    right = np.minimum(left + 1, len(parent_times) - 1)
+    intervals = parent_times[right] - parent_times[left]
+    intervals[left == len(parent_times) - 1] = 1.0 / parent_rate
+    return parent_times[left] + weight * intervals
+
+
 resample = define_function_handling_dict_from_class(source_class=ResampleRecording, name="resample")
diff --git a/src/spikeinterface/preprocessing/tests/test_decimate.py b/src/spikeinterface/preprocessing/tests/test_decimate.py
index c0f535d61f..b62ec43bdf 100644
--- a/src/spikeinterface/preprocessing/tests/test_decimate.py
+++ b/src/spikeinterface/preprocessing/tests/test_decimate.py
@@ -5,7 +5,7 @@
 
 from spikeinterface import NumpyRecording
 from spikeinterface.core import generate_recording, load
-from spikeinterface.preprocessing.decimate import DecimateRecording, decimate, get_balanced_decimation_factors
+from spikeinterface.preprocessing.decimate import DecimateRecording, decimate
 from spikeinterface.preprocessing.resample import ResampleRecording
 from spikeinterface.preprocessing.tests.test_resample import create_sinusoidal_traces
 import numpy as np
@@ -89,46 +89,29 @@ def test_decimate_with_times(antialias):
         )
 
 
-@pytest.mark.parametrize(
-    "decimation_factor, max_factor, expected",
-    [
-        (1, 13, [1]),
-        (7, 13, [7]),
-        (13, 13, [13]),
-        (48, 13, [8, 6]),
-        (50, 13, [10, 5]),
-        (60, 13, [10, 6]),
-        (100, 13, [10, 10]),
-        (17, 13, [17]),  # prime > 13: cannot be split
-        (23, 13, [23]),
-        (48, 8, [8, 6]),
-        (48, 4, [4, 4, 3]),
-        (10, 4, [10]),  # prime > the custom limit: cannot be split
-    ],
-)
-def test_balanced_decimation_factors(decimation_factor, max_factor, expected):
-    if expected == [decimation_factor] and decimation_factor > max_factor:
-        with pytest.warns(UserWarning, match=f"prime factor > {max_factor}"):
-            factors = get_balanced_decimation_factors(decimation_factor, max_factor=max_factor)
-    else:
-        factors = get_balanced_decimation_factors(decimation_factor, max_factor=max_factor)
-    assert factors == expected
-    # The product of the sub-factors always reconstructs the requested factor.
-    assert int(np.prod(factors)) == decimation_factor
-    # Every pass respects the limit unless no valid split exists.
-    if len(factors) > 1:
-        assert all(f <= max_factor for f in factors)
+@pytest.mark.parametrize("factor", [1, 7, 17, 48, 300])
+def test_decimate_polyphase(factor):
+    from scipy.signal import resample_poly
+
+    traces = np.random.default_rng(4621).standard_normal((1001, 2))
+    rec = NumpyRecording(traces, 30000)
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always")
+        decimated = decimate(rec, factor, antialias=True)
+        resampled = ResampleRecording(rec, 30000 / factor)
+    assert not caught
+    expected = resample_poly(traces, 1, factor, axis=0, padtype="reflect")
+    np.testing.assert_allclose(decimated.get_traces(), expected, rtol=1e-12, atol=1e-12)
+    np.testing.assert_array_equal(resampled.get_traces(), decimated.get_traces())
+    np.testing.assert_array_equal(resampled.get_times(), decimated.get_times())
 
 
 @pytest.mark.parametrize("decimation_factor", [6, 10, 48, 300])
 def test_decimate_antialias_by_chunks(decimation_factor):
-    # Mirror test_resample_by_chunks: chunked reads must match a full read once the
-    # anti-aliasing margins are accounted for. Factor 48 exercises the internal multi-pass.
     sampling_frequency = int(3e4)
     duration = 30
     traces, _ = create_sinusoidal_traces(sampling_frequency, duration, freqs_n=10, max_freq=1000, dtype=np.float32)
     parent_rec = NumpyRecording(traces, sampling_frequency)
-    rms = np.sqrt(np.mean(parent_rec.get_traces() ** 2))
     decimated_rate = sampling_frequency / decimation_factor
 
     for margin_ms in [None, 100, 1000]:
@@ -139,13 +122,7 @@ def test_decimate_antialias_by_chunks(decimation_factor):
         traces2 = rec2.get_traces()
         traces3 = rec3.get_traces()
 
-        # Drop the first and last chunk before comparing (as in test_resample_by_chunks).
-        sl = slice(chunk_size, -chunk_size)
-        error_mean = np.sqrt(np.mean((traces2[sl] - traces3[sl]) ** 2))
-        error_max = np.sqrt(np.max((traces2[sl] - traces3[sl]) ** 2))
-
-        assert error_mean / rms < 0.01
-        assert error_max / rms < 0.05
+        np.testing.assert_array_equal(traces2, traces3)
 
 
 @pytest.mark.parametrize("decimation_factor", [6, 10])
@@ -163,6 +140,8 @@ def test_decimate_antialias_with_offset(decimation_factor, decimation_offset):
         parent_rec, decimation_factor, decimation_offset=decimation_offset, antialias=False, dtype="float32"
     )
 
+    np.testing.assert_allclose(dec_aa.get_times(), parent_rec.get_times()[decimation_offset::decimation_factor])
+
     # The anti-aliasing path returns the same number of samples as plain slicing.
     parent_n = parent_rec.get_num_samples()
     expected_n = int(np.ceil((parent_n - decimation_offset) / decimation_factor))
@@ -174,7 +153,7 @@ def test_decimate_antialias_with_offset(decimation_factor, decimation_offset):
     assert corr > 0.95
 
 
-def test_decimate_antialias_multipass():
+def test_decimate_polyphase_serialization():
     sampling_frequency = 30000
     decimation_factor = 48
     traces, _ = create_sinusoidal_traces(sampling_frequency, duration=10, freqs_n=8, max_freq=200, dtype=np.float32)
@@ -182,14 +161,9 @@ def test_decimate_antialias_multipass():
 
     dec = decimate(parent_rec, decimation_factor, antialias=True)
 
-    # Multi-pass happens internally: a single DecimateRecording carries the full factor.
     assert isinstance(dec, DecimateRecording)
     assert dec._kwargs["decimation_factor"] == decimation_factor
 
-    segment = dec.segments[0]
-    assert int(np.prod(segment._decimation_factors)) == decimation_factor
-    assert all(f <= 13 for f in segment._decimation_factors)
-
     parent_n = parent_rec.get_num_samples()
     assert dec.get_num_samples() == int(np.ceil(parent_n / decimation_factor))
 
@@ -198,22 +172,5 @@ def test_decimate_antialias_multipass():
     np.testing.assert_allclose(dec_loaded.get_traces(), dec.get_traces())
 
 
-def test_decimate_antialias_large_prime_warns():
-    rec = generate_recording(durations=[2.0], num_channels=2, sampling_frequency=34000)
-    for preprocess in [
-        lambda: DecimateRecording(rec, 17, antialias=True),
-        lambda: ResampleRecording(rec, rec.get_sampling_frequency() / 17),
-    ]:
-        with pytest.warns(UserWarning, match="prime factor > 13") as caught:
-            dec = preprocess()
-        assert len(caught) == 1
-        assert dec.segments[0]._decimation_factors == [17]
-
-    with warnings.catch_warnings(record=True) as caught:
-        warnings.simplefilter("always")
-        DecimateRecording(rec, 17, antialias=False)
-    assert not caught
-
-
 if __name__ == "__main__":
     test_decimate()
diff --git a/src/spikeinterface/preprocessing/tests/test_resample.py b/src/spikeinterface/preprocessing/tests/test_resample.py
index 1d9a59c88a..e116b90bd7 100644
--- a/src/spikeinterface/preprocessing/tests/test_resample.py
+++ b/src/spikeinterface/preprocessing/tests/test_resample.py
@@ -100,7 +100,7 @@ def test_resample_freq_domain():
     parent_rec = NumpyRecording(traces, sampling_frequency)
     # Different resampling frequencies, always below Niquist
     resamp_fss = (np.linspace(0.1, 0.45, 10) * sampling_frequency).astype(int)
-    resamp_recs = [resample(parent_rec, resamp_fs) for resamp_fs in resamp_fss]
+    resamp_recs = [resample(parent_rec, resamp_fs, max_denominator=30000) for resamp_fs in resamp_fss]
     # First set of tests, we are updating frames and time duration correctly
 
     # that they all have the correct number of frames:
@@ -155,74 +155,67 @@ def test_resample_freq_domain():
 
 
 def test_resample_by_chunks():
-    # Now tests the margin effects and the chunk_sizes for sanity
-    # The same as in the phase_shift tests.
-    sampling_frequency = int(3e4)
-    duration = 30
-    freqs_n = 10
-    # ~ dtype = np.int16
-    dtype = np.float32
-    max_freq = 1000
-    traces, [freqs_vals, amps_vals, phase_shifts] = create_sinusoidal_traces(
-        sampling_frequency, duration, freqs_n, max_freq, dtype
-    )
-    parent_rec = NumpyRecording(traces, sampling_frequency)
-    rms = np.sqrt(np.mean(parent_rec.get_traces() ** 2))
-    # The chunk_size must be always at least some 1 second of the resample, else it breaks.
-    # Does this makes sense?
-    # Also, sometimes decimate might give warnings about filter designs
-    resample_rates = [1000, 2000]  # [500, 1000, 2500]
-    margins_ms = [100, 1000]  # [100, 200, 1000]
-    chunk_durations = [0.5, 1]  # [1, 2, 3]
-
-    for resample_rate in resample_rates:
-        for margin_ms in margins_ms:
-            for chunk_size in [int(resample_rate * chunk_multi) for chunk_multi in chunk_durations]:
-                # print(f'resmple_rate = {resample_rate}; margin_ms = {margin_ms}; chunk_size={chunk_size}')
-                rec2 = resample(parent_rec, resample_rate, margin_ms=margin_ms)
-                # save by chunk rec3 is the cached version
-                rec3 = rec2.save(format="memory", chunk_size=chunk_size, n_jobs=1, progress_bar=False)
-
-                traces2 = rec2.get_traces()
-                traces3 = rec3.get_traces()
-
-                # error between full and chunked
-                # for error first and last chunk is removed
-                sl = slice(chunk_size, -chunk_size)
-                error_mean = np.sqrt(np.mean((traces2[sl] - traces3[sl]) ** 2))
-                error_max = np.sqrt(np.max((traces2[sl] - traces3[sl]) ** 2))
-
-                # this will never be possible:
-                # assert np.allclose(traces2, traces3)
-                # so we check that the diff between chunk processing and not chunked is small
-                # print()
-                # print(dtype, margin_ms, chunk_size)
-                # print(error_mean, rms, error_mean / rms)
-                # print(error_max, rms, error_max / rms)
-                # The original thrshold are too restrictive, but in all cases
-                # The signals look quite similar, with error that are small enough
-                # But, when using signal.resample, the last edge becomes too noisy
-
-                assert error_mean / rms < 0.01
-                assert error_max / rms < 0.05
-
-                if DEBUG:
-                    fig, axs = plt.subplots(nrows=2, sharex=True)
-                    fig.suptitle(
-                        f"Resample rate {resample_rate}\nMargin {margin_ms}\nChunk size {chunk_size}\n error mean(%) {error_mean / rms}  error max(%){error_max / rms} "
-                    )
-                    ax = axs[0]
-                    ax.plot(traces2, color="g", label="no chunk")
-                    ax.plot(traces3, color="r", label=f"chunked")
-                    for i in range(traces2.shape[0] // chunk_size):
-                        ax.axvline(chunk_size * i, color="k", alpha=0.4)
-                    ax.legend()
-                    ax = axs[1]
-                    ax.plot(traces3 - traces2)
-                    for i in range(traces2.shape[0] // chunk_size):
-                        ax.axvline(chunk_size * i, color="k", alpha=0.4)
-
-                    plt.show()
+    traces = np.random.default_rng(4621).standard_normal((60001, 2)).astype("float32")
+    parent_rec = NumpyRecording(traces, 30000)
+    for rate in [1000, 700, 45000, 333.3]:
+        for margin_ms in [None, 0, 100]:
+            processed = resample(parent_rec, rate, margin_ms=margin_ms)
+            saved = processed.save(format="memory", chunk_size=137, n_jobs=1, progress_bar=False)
+            np.testing.assert_array_equal(saved.get_traces(), processed.get_traces())
+
+
+def test_resample_rational_grid():
+    # "the rational sample grid" = "the positions of output samples expressed in input-sample 
+    # coordinates, using the ratio up/down"
+    from contextlib import nullcontext
+    from scipy.signal import resample_poly
+    from spikeinterface.core import load
+
+    traces = np.random.default_rng(4621).standard_normal((1001, 2))
+    for parent_rate, requested_rate, up, down, max_denominator in [
+        (1000, 333.3, 3333, 10000, 10000),
+        (1000, 333.3, 1, 3, 10),
+        (1000, 1500, 3, 2, 10000),
+        (30000.03, 1000, 1, 30, 10000),
+        (625, 625 / 6, 1, 6, 10000),
+        (1000, 1000, 1, 1, 10000),
+    ]:
+        parent = NumpyRecording(traces, parent_rate, t_starts=[10.0])
+        expected = resample_poly(traces, up, down, axis=0, padtype="reflect")
+        achieved_rate = parent_rate * up / down
+        approximate = not np.isclose(achieved_rate, requested_rate, rtol=1e-12, atol=0)
+        for explicit_times in [False, True]:
+            if explicit_times:
+                parent.set_times(parent.get_times(), with_warning=False)
+            with pytest.warns(UserWarning, match="achieves") if approximate else nullcontext():
+                processed = resample(parent, requested_rate, max_denominator=max_denominator)
+            assert processed.get_sampling_frequency() == pytest.approx(achieved_rate, rel=1e-15)
+            np.testing.assert_allclose(processed.get_traces(), expected, rtol=0, atol=1e-12)
+            np.testing.assert_allclose(
+                processed.get_times(), 10.0 + np.arange(len(expected)) / achieved_rate, rtol=0, atol=1e-12
+            )
+            assert processed._kwargs["resample_rate"] == requested_rate
+            assert processed._kwargs["max_denominator"] == max_denominator
+            with pytest.warns(UserWarning, match="achieves") if approximate else nullcontext():
+                restored = load(processed.to_dict())
+            assert restored.get_sampling_frequency() == processed.get_sampling_frequency()
+            np.testing.assert_array_equal(restored.get_traces(), processed.get_traces())
+
+
+def test_resample_short_sections():
+    for rate in [100, 1500]:
+        parent = NumpyRecording(np.array([[1.0], [9.0]]), 1000)
+        parent.set_times(np.array([10.0, 20.0]), with_warning=False)
+        processed = resample(parent, rate, gap_tolerance_ms=0)
+        sections = [resample(parent.frame_slice(i, i + 1), rate) for i in range(2)]
+        expected = np.concatenate([section.get_traces() for section in sections])
+        np.testing.assert_array_equal(processed.get_traces(), expected)
+        np.testing.assert_array_equal(processed.get_times(), np.concatenate([s.get_times() for s in sections]))
+        for i in range(len(expected)):
+            np.testing.assert_array_equal(processed.get_traces(start_frame=i, end_frame=i + 1), expected[i : i + 1])
+        assert processed.get_traces(start_frame=1, end_frame=1).shape == (0, 1)
+        empty = resample(NumpyRecording(np.empty((0, 1)), 1000), rate)
+        assert empty.get_traces().shape == (0, 1)
 
 
 def test_resample_preserves_t_start():
@@ -301,17 +294,12 @@ def test_resample_preserves_time_vector_non_integer_ratio():
     time_vector = np.arange(n_samples, dtype="float64") / sampling_frequency + 10.0
     parent_rec.set_times(time_vector)
 
-    import warnings as _warnings
-
-    with _warnings.catch_warnings(record=True) as w:
-        _warnings.simplefilter("always")
-        resampled = resample(parent_rec, resample_rate)
-        assert any("non-integer ratio" in str(warning.message).lower() for warning in w)
+    resampled = resample(parent_rec, resample_rate)
 
     assert resampled.has_time_vector()
     resampled_times = resampled.get_times()
     assert len(resampled_times) == resampled.get_num_samples()
-    assert np.isclose(resampled_times[0], 10.0, atol=1.0 / sampling_frequency)
+    np.testing.assert_allclose(resampled_times, 10.0 + np.arange(len(resampled_times)) / resample_rate, atol=1e-12)
 
 
 def test_resample_errors_on_gaps_by_default():
@@ -360,24 +348,24 @@ def test_resample_preserves_gaps_non_integer_ratio():
     [
         700,  # non-integer ratio (30000 / 700 ~= 42.857)
         500,  # integer ratio (30000 / 500 = 60)
-        625,  # factor 48 does not divide the default 3000-sample margin
+        625,  # integer-factor downsampling
+        45000,  # upsampling
     ],
-    ids=["non_integer_ratio", "integer_ratio", "integer_ratio_unaligned_margin"],
+    ids=["non_integer_ratio", "integer_ratio", "factor_48", "upsampling"],
 )
 def test_resample_traces_across_gap(resample_rate):
     """Section-wise resampling should match individually resampled sections.
 
     Build a gapped recording, resample it with gap_tolerance_ms, and verify
     that each section's output matches what you'd get by resampling that
-    section alone (without the gap).  This confirms that _get_traces_gapped
-    does not apply FFT processing (or decimate filtering) across gap boundaries.
+    section alone (without the gap). Filtering must not cross gap boundaries.
     """
     sampling_frequency = 30000
     sec_duration = 2.0
     gap_s = 5.0
 
-    n1 = int(sec_duration * sampling_frequency)
-    n2 = int(sec_duration * sampling_frequency)
+    n1 = int(sec_duration * sampling_frequency) + 7
+    n2 = int(sec_duration * sampling_frequency) + 13
 
     # Build random traces (more realistic than a sinusoid)
     rng = np.random.default_rng(42)
@@ -437,19 +425,18 @@ def test_resample_traces_across_gap(resample_rate):
     assert gapped_s2.shape == ref_traces2.shape, f"Section 2 shape mismatch: {gapped_s2.shape} vs {ref_traces2.shape}"
     np.testing.assert_allclose(gapped_s2, ref_traces2, rtol=1e-5, atol=1e-5)
 
-    if sampling_frequency % resample_rate == 0:
-        # Exercise reads within a section and across the gap, including channel selection.
-        for start, end in [(123, 456), (n_out_1 - 123, n_out_1 + 456)]:
-            expected_pieces = []
-            for offset, section in [(0, resampled1), (n_out_1, resampled2)]:
-                local_start = max(0, start - offset)
-                local_end = min(section.get_num_samples(), end - offset)
-                if local_start < local_end:
-                    expected_pieces.append(
-                        section.get_traces(start_frame=local_start, end_frame=local_end, channel_ids=[1])
-                    )
-            actual = resampled.get_traces(start_frame=start, end_frame=end, channel_ids=[1])
-            np.testing.assert_allclose(actual, np.concatenate(expected_pieces), rtol=1e-5, atol=1e-5)
+    # Exercise reads within a section and across the gap, including channel selection.
+    for start, end in [(123, 456), (n_out_1 - 123, n_out_1 + 456)]:
+        expected_pieces = []
+        for offset, section in [(0, resampled1), (n_out_1, resampled2)]:
+            local_start = max(0, start - offset)
+            local_end = min(section.get_num_samples(), end - offset)
+            if local_start < local_end:
+                expected_pieces.append(
+                    section.get_traces(start_frame=local_start, end_frame=local_end, channel_ids=[1])
+                )
+        actual = resampled.get_traces(start_frame=start, end_frame=end, channel_ids=[1])
+        np.testing.assert_allclose(actual, np.concatenate(expected_pieces), rtol=1e-5, atol=1e-5)
 
 
 @pytest.mark.parametrize("resample_rate", [700, 625])

From 1b37035fb32d99d331d050627afeed37509e851d Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 10 Sep 2026 19:26:31 +0000
Subject: [PATCH 11/11] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 src/spikeinterface/preprocessing/resample.py   | 18 +++++++++---------
 .../preprocessing/tests/test_resample.py       |  2 +-
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/spikeinterface/preprocessing/resample.py b/src/spikeinterface/preprocessing/resample.py
index c4c9b11abd..3123a4d69f 100644
--- a/src/spikeinterface/preprocessing/resample.py
+++ b/src/spikeinterface/preprocessing/resample.py
@@ -23,10 +23,10 @@ class ResampleRecording(BasePreprocessor):
     recording : Recording
         The recording extractor to be re-referenced
     resample_rate : int | float
-        The requested sampling frequency. The closest output/input rate ratio with 
-        denominator at most `max_denominator` is selected. The output reports the 
+        The requested sampling frequency. The closest output/input rate ratio with
+        denominator at most `max_denominator` is selected. The output reports the
         achieved rate, while the requested rate is retained in serialization kwargs.
-        A relative difference exceeding 1e-12 emits a warning. 
+        A relative difference exceeding 1e-12 emits a warning.
     gap_tolerance_ms : float | None, default: None
         Maximum acceptable gap size in milliseconds for automatic segmentation.
 
@@ -61,7 +61,7 @@ class ResampleRecording(BasePreprocessor):
     max_denominator : int, default: 10000
         Maximum denominator of the rational output/input rate ratio. Increasing this can
         improve rate accuracy, but can also increase filter length and the aligned input
-        span needed for a chunk. 
+        span needed for a chunk.
 
     Returns
     -------
@@ -71,13 +71,13 @@ class ResampleRecording(BasePreprocessor):
     Notes
     -----
     Each section returns ``ceil(num_input_samples * up / down)`` samples, matching SciPy
-    and ``decimate()``. Output timestamps use the same rational grid as the traces. 
-    Explicit parent timestamps are sampled or interpolated within each section. 
-    Output positions beyond the last input sample extrapolate its timestamp by less 
+    and ``decimate()``. Output timestamps use the same rational grid as the traces.
+    Explicit parent timestamps are sampled or interpolated within each section.
+    Output positions beyond the last input sample extrapolate its timestamp by less
     than one nominal input period.
 
-    For example, resampling from 30000.01 Hz to a requested 2500 Hz with the default 
-    denominator limit selects up=1 and down=12. The output reports 2500.000833333333 Hz, 
+    For example, resampling from 30000.01 Hz to a requested 2500 Hz with the default
+    denominator limit selects up=1 and down=12. The output reports 2500.000833333333 Hz,
     and a warning reports the difference of approximately +0.333333 ppm.
 
     """
diff --git a/src/spikeinterface/preprocessing/tests/test_resample.py b/src/spikeinterface/preprocessing/tests/test_resample.py
index e116b90bd7..b9450b92e3 100644
--- a/src/spikeinterface/preprocessing/tests/test_resample.py
+++ b/src/spikeinterface/preprocessing/tests/test_resample.py
@@ -165,7 +165,7 @@ def test_resample_by_chunks():
 
 
 def test_resample_rational_grid():
-    # "the rational sample grid" = "the positions of output samples expressed in input-sample 
+    # "the rational sample grid" = "the positions of output samples expressed in input-sample
     # coordinates, using the ratio up/down"
     from contextlib import nullcontext
     from scipy.signal import resample_poly