diff --git a/doc/api.rst b/doc/api.rst index a51e26ac91..2e26aba251 100755 --- a/doc/api.rst +++ b/doc/api.rst @@ -217,6 +217,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/doc/modules/preprocessing.rst b/doc/modules/preprocessing.rst index e587ac45a5..97958cd4d0 100644 --- a/doc/modules/preprocessing.rst +++ b/doc/modules/preprocessing.rst @@ -53,6 +53,7 @@ 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. + The Preprocessing Pipeline -------------------------- 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 716cc46f70..358dbce460 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,19 @@ from .basepreprocessor import BasePreprocessor from .filter import fix_dtype -from spikeinterface.core import BaseRecordingSegment +from ._resampling_tools import get_polyphase_filter +from spikeinterface.core import BaseRecordingSegment, get_chunk_with_margin 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 + (``[::]``), 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.resample_poly`` (the + same anti-aliased decimation used by ``spikeinterface.preprocessing.ResampleRecording``). Parameters ---------- @@ -29,12 +34,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 | None, default: None + If True, apply an anti-aliasing low-pass filter before downsampling, using + ``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 + 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. 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: ``` = [::]``` """ @@ -44,14 +62,17 @@ def __init__( recording, decimation_factor, decimation_offset=0, + antialias=None, + margin_ms=None, + dtype=None, ): # Original sampling frequency self._orig_samp_freq = recording.get_sampling_frequency() 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())] ) @@ -63,7 +84,26 @@ 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 + + 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 + + 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) for parent_segment in recording.segments: self.add_recording_segment( @@ -74,6 +114,9 @@ def __init__( decimation_factor, decimation_offset, self._dtype, + antialias, + margin, + filter_coefficients, ) ) @@ -81,6 +124,9 @@ def __init__( recording=recording, decimation_factor=decimation_factor, decimation_offset=decimation_offset, + antialias=antialias, + margin_ms=margin_ms, + dtype=dtype, ) @@ -93,6 +139,9 @@ def __init__( decimation_factor, decimation_offset, dtype, + antialias=False, + margin=0, + filter_coefficients=None, ): if parent_recording_segment._time_vector is not None: time_vector = parent_recording_segment._time_vector[decimation_offset::decimation_factor] @@ -100,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__( @@ -113,25 +161,106 @@ def __init__( self._decimation_factor = decimation_factor self._decimation_offset = decimation_offset self._dtype = dtype + self._antialias = antialias + self._margin = margin + 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): - # 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) + + return get_polyphase_resampled_traces( + self._parent_segment, + start_frame, + end_frame, channel_indices, - )[ - :: self._decimation_factor - ].astype(self._dtype) + 1, + self._decimation_factor, + self._margin, + self._dtype, + self._filter_coefficients, + decimation_offset=self._decimation_offset, + ) + + +def get_polyphase_resampled_traces( + parent_segment, + start_frame, + end_frame, + channel_indices, + up, + down, + margin, + dtype, + filter_coefficients, + decimation_offset=0, +): + """Resample a chunk, with reflected boundary padding.""" + from scipy.signal import resample_poly + + if end_frame <= start_frame: + return parent_segment.get_traces(0, 0, channel_indices).astype(dtype) + + 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, + channel_indices, + margin, + add_reflect_padding=True, + ) + working_dtype = np.result_type(parent_traces.dtype, dtype, np.float32) + 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): + """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 a801c45eff..3123a4d69f 100644 --- a/src/spikeinterface/preprocessing/resample.py +++ b/src/spikeinterface/preprocessing/resample.py @@ -1,31 +1,32 @@ 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 -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 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. + 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 - The resampling frequency + 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 + 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. @@ -48,18 +49,37 @@ 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, will be used to decrease edge effects. + margin_ms : float | None, default: None + 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. - skip_checks : bool, default: False - If True, checks on sampling frequencies and cutoff filter frequencies are skipped + 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 The resampled recording extractor object. + 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 + 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. + """ def __init__( @@ -67,36 +87,31 @@ def __init__( recording, resample_rate, gap_tolerance_ms=None, - margin_ms=100.0, + margin_ms=None, dtype=None, - skip_checks=False, + max_denominator=10000, ): - # 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 + 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 - # 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 = int(margin_ms * recording.get_sampling_frequency() / 1000) + 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, + up, + down, + filter_coefficients, ) ) @@ -106,7 +121,7 @@ def __init__( gap_tolerance_ms=gap_tolerance_ms, margin_ms=margin_ms, dtype=dtype, - skip_checks=skip_checks, + max_denominator=max_denominator, ) @@ -118,7 +133,10 @@ def __init__( parent_rate, margin, dtype, - gap_tolerance_ms=None, + gap_tolerance_ms, + up, + down, + filter_coefficients, ): self._resample_rate = resample_rate self._parent_segment = parent_recording_segment @@ -126,6 +144,9 @@ def __init__( self._margin = margin self._dtype = dtype self._has_gaps = False + 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! @@ -176,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, @@ -189,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: @@ -242,177 +225,69 @@ 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) - # Original code path: no gaps (or no time_vector) - # get parent traces with margin - 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( - self._parent_segment, - parent_start_frame, - parent_end_frame, + 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): + return get_polyphase_resampled_traces( + parent_segment, + start_frame, + end_frame, channel_indices, + self._up, + self._down, self._margin, - add_reflect_padding=True, - dtype=np.float32, + self._dtype, + self._filter_coefficients, ) - # 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 in case of resample: - 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) - - # now take care of the edges - resampled_traces = resampled_traces[left_margin_rs : num - right_margin_rs] - return resampled_traces.astype(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._parent_rate % self._resample_rate) == 0 - - 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 - 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) - - 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 = self._margin - left_margin - pad_right = self._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 - - # 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) - - # 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 - if is_integer_ratio: - q = int(self._parent_rate / self._resample_rate) - resampled = decimate(parent_traces, q=q, axis=0) - if np.any(np.isnan(resampled)): - resampled = resample(parent_traces, num, axis=0) - else: - resampled = resample(parent_traces, num, axis=0) + return result - # 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] +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") - - -# 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]) diff --git a/src/spikeinterface/preprocessing/tests/test_decimate.py b/src/spikeinterface/preprocessing/tests/test_decimate.py index 93f70ea9bd..b62ec43bdf 100644 --- a/src/spikeinterface/preprocessing/tests/test_decimate.py +++ b/src/spikeinterface/preprocessing/tests/test_decimate.py @@ -1,12 +1,24 @@ +import warnings + import pytest 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 +from spikeinterface.preprocessing.resample import ResampleRecording +from spikeinterface.preprocessing.tests.test_resample import create_sinusoidal_traces 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]) @@ -45,7 +57,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 +68,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 +81,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 +89,88 @@ def test_decimate_with_times(): ) +@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): + 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) + decimated_rate = sampling_frequency / decimation_factor + + 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) + + traces2 = rec2.get_traces() + traces3 = rec3.get_traces() + + np.testing.assert_array_equal(traces2, traces3) + + +@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" + ) + + 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)) + 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_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) + parent_rec = NumpyRecording(traces, sampling_frequency) + + dec = decimate(parent_rec, decimation_factor, antialias=True) + + assert isinstance(dec, DecimateRecording) + assert dec._kwargs["decimation_factor"] == decimation_factor + + 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()) + + if __name__ == "__main__": test_decimate() diff --git a/src/spikeinterface/preprocessing/tests/test_resample.py b/src/spikeinterface/preprocessing/tests/test_resample.py index 089905aaee..b9450b92e3 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,23 +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, # integer-factor downsampling + 45000, # upsampling ], - ids=["non_integer_ratio", "integer_ratio"], + 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) @@ -436,11 +425,24 @@ 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) - -def test_resample_gapped_chunked_consistency(): + # 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]) +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