Add antialiased decimation and allow non-integer resampling rates

[grahamfindlay]

Implements #4605 : Adds antialised decimation (multipass with balanced factors, when possible) to both DecimateRecording and ResampleRecording. Opt-in for DecimateRecording with the antialias parameter (preserves existing default behavior). For ResampleRecording, antialiasing is always applied. For a future release where existing default behavior can be changed, I would recommend making antalias=True the default.

For simplicity, we could just remove the fast decimation path from ResampleRecording entirely and instead, if an integer resampling ratio is detected, just suggest to the user that they use DecimateRecording instead. But, that would change existing behavior, and I don't think this is too terribly complicated as-is.

[grahamfindlay]

@alejoe91 Fixed (surface-level) conflicts introduced by merging #4499 into main. All good to go.

[JoeZiminski]

Hi @grahamfindlay this looks great! Some small suggestions on the code. I didn't have a chance to play around locally and look at the tests but will do that early next week! In the meantime just some general questions:

1. Do you think it is worth putting a warning on antialias=False with possibly a timeline to switching the default to True? I would be tempted to change it for 105 but it's such a big change it might trip people up! But I guess it is important for people to use the AA filter especially if going from AP -> LFP band which I guess is the most common use of the decimate function?

2. What do you think of resample_poly? Would solve a few issues (more predictable artefact than FFT method, natively handles large decimation factors, can funnel all resample calls to it rather than splitting decimate by integer/non-integer)) but maybe tricky to always find nice up/down arguments from the resample ratio depending on the sampling rates.

3. Currently it warns on NaN if found after decimate, but should we error? I think (?) that NaN in traces are not supported in downstream methods e.g. np.nan... methods are not typically used.

4. The methods to handle time-gaps in the recording is really elegant. In the wider-picture of SI though, I wonder if it is necessary to support this mode, and if the module can become more readable by stripping it out. The reason being that AFAIK SI makes no guarantees on supporting gaps in the time_vector for preprocessing, and most steps do not consider it nor work properly for it (e.g. filtering steps, or even decimate.py). My understanding is that in SI time gaps should be handled by the user by using segments or separate recordings. I think in the resample.py case we could just check for time gaps and error accordingly. That being said maybe there are important use cases for this I am overlooking.

Cheers!

[JoeZiminski]

Could this check be centralised in _decimation_tools.py and also called from decimate.py? (and/or antialias_factors there called decimation_factors for consistencey)

[grahamfindlay]

yes, I would prefer that -- will do

[JoeZiminski]

Suggested change

	# Non-integer ratio, no gaps: FFT-based resampling with proportional margins.
	# Non-integer ratio downsampling or upsampling with no gaps: FFT-based resampling with proportional margins.

[JoeZiminski]

Suggested change

	# Integer ratio, no gaps: anti-aliased multi-pass decimation.
	# Integer ratio downsampling, no gaps: anti-aliased multi-pass decimation.

[JoeZiminski]

Here if Nan detected we warn, but if Nan detected in the section implementation we resample. Can these align?

[JoeZiminski]

Is it worth doing something like:

    np.round(decimated_traces, out=decimated_traces)
    info = np.iinfo(dtype)
    np.clip(decimated_traces, info.min, info.max, out=decimated_traces)

just in case overflow, although I guess the risk is very low with low-pass filter unless we encounter ringing artefact

[JoeZiminski]

This function is build around the 13 recommendation of Scipy. Is it easy to / worth generalizing, and passing _MAX_SINGLE_PASS_DECIMATION as an argument?

[JoeZiminski]

This margin seems a good catch-all. I guess the main use case of this will be going from AP band to LFP band? I think this margin will only be a problem based on going to very low sampling rates:

(this plot was generated by ChatGPT using the below sensible code to measure the impulse response from the corresponding 8th order Bessel). I guess we could add a warning for low reample_rate and/or perform a check similar to the one below directly from the filter (probably overkill)?

Details

import numpy as np
from scipy import signal


def estimate_iir_margin(fs, fout, eps=1e-6):
    q = fs / fout
    assert np.isclose(q, round(q))
    q = int(round(q))

    # Same default IIR design used by scipy.signal.decimate
    sos = signal.cheby1(
        8,
        0.05,
        0.8 / q,
        output="sos",
    )

    # Long impulse, far away from either edge
    n = 200_000
    x = np.zeros(n)
    centre = n // 2
    x[centre] = 1.0

    # Effective zero-phase filter
    h = signal.sosfiltfilt(sos, x)

    threshold = eps * np.max(np.abs(h))

    indices = np.flatnonzero(np.abs(h) > threshold)

    radius_samples = max(
        centre - indices[0],
        indices[-1] - centre,
    )

    margin_ms = radius_samples / fs * 1000

    return q, radius_samples, margin_ms


fs = 30_000

for fout in [2500, 1000, 500, 250, 100]:
    q, samples, margin_ms = estimate_iir_margin(fs, fout)

    print(
        f"{fs} -> {fout:4} Hz "
        f"q={q:3}, "
        f"margin={samples:5} samples = {margin_ms:.1f} ms"
    )

[grahamfindlay]

I think I just copied this margin code from what was already used elsewhere in the codebase, for consistency.
Yes, it is definitely worth checking. I think the code there that produced the table uses Chebyshev type I, not Bessel, and measures a single stage rather than the cascade. But in any case, I would test the actual cascade and then choose an automatic margin (or a warning threshold). Probably also makes sense to retain some kind of override. Will do.

[JoeZiminski]

Suggested change

	# When the parent rate is an exact integer multiple of resample_rate, get_traces uses
	# When the parent rate is an exact integer multiple of resample_rate (i.e. downsampling with an integer factor), get_traces uses

[JoeZiminski]

Could this move to decimate.py? (if it can be imported from resample.py). Not a huge thing but I think the convention is for get-traces like functionality to live in next to the segments.

[grahamfindlay]

Yes, will do

[grahamfindlay]

Hi @grahamfindlay this looks great! Some small suggestions on the code. I didn't have a chance to play around locally and look at the tests but will do that early next week! In the meantime just some general questions:

1. Do you think it is worth putting a warning on antialias=False with possibly a timeline to switching the default to True? I would be tempted to change it for 105 but it's such a big change it might trip people up! But I guess it is important for people to use the AA filter especially if going from AP -> LFP band which I guess is the most common use of the decimate function?

Yeah, I think that a transition warning is a good idea. I might distinguish an omitted argument from an explicit antialias=False, though, in case someone has deliberately applied an appropriate low-pass filter upstream and doesn't deserve repeated warnings. I will add this.

2. What do you think of resample_poly? Would solve a few issues (more predictable artefact than FFT method, natively handles large decimation factors, can funnel all resample calls to it rather than splitting decimate by integer/non-integer)) but maybe tricky to always find nice up/down arguments from the resample ratio depending on the sampling rates.

I definitely would support the use of resample_poly. It's closer to what I actually use with my own data. I think the most disruptive thing about it would be that currently we always deliver the exact requested target resample rate¹, whereas with resample_poly we would only promise to deliver something very close. Personally, I'd rather have that, but maybe others feel differently. Happy to implement it if anyone gives me the green light.

3. Currently it warns on NaN if found after decimate, but should we error? I think (?) that NaN in traces are not supported in downstream methods e.g. np.nan... methods are not typically used.

Yes, I agree. Actually, I think it's even worse than NaN not being supported downstream: casting to an integer type will turn the NaN to zeros, which will let downstream operations proceed where they probably should not... FWIW, I would also remove the silent FFT retry from the section-wise path. I know it's nice when things "just work" for the user, but I don't love the idea that the resampling method for a chunk can change silently because of a rare numerical failure.

4. The methods to handle time-gaps in the recording is really elegant. In the wider-picture of SI though, I wonder if it is necessary to support this mode, and if the module can become more readable by stripping it out. The reason being that AFAIK SI makes no guarantees on supporting gaps in the time_vector for preprocessing, and most steps do not consider it nor work properly for it (e.g. filtering steps, or even decimate.py). My understanding is that in SI time gaps should be handled by the user by using segments or separate recordings. I think in the resample.py case we could just check for time gaps and error accordingly. That being said maybe there are important use cases for this I am overlooking.

My (and other people's) chronic, continuous, multi-day recordings commonly have small gaps. I am very much open to representing continuous stretches as segments or separate recordings, as long as it provides a supported route through preprocessing, sorting, postprocessing, analyzer objects, etc. Last I checked, adding limited qualified support for gaps within a segment was by far the easiest and least disruptive way to achieve this, even though SI does generally expect contiguous segments. It is probably worth mentioning that this expectation is rarely if ever checked or enforced, and there is explicit support for and discussion of missing frames and irregularly-sampled data scattered all over the place, e.g. in tutorials, tests, in the NWB recording extractor, the Blackrock reader, etc. I guess that is part of why the time vector exists, rather than always assuming a start time + fixed sample rate. So I would not say there is any contract in SI that says a segment must be gap-free and regularly-sampled, even if that is asserted in the docs once or twice. It is maybe most accurate to say that all operations are expected to work for gap-free, regularly-sampled segments, and support for imperfect segments depends on the operation.

It is also worth mentioning that motion estimation, lots of sorter-specific processing, and timing restoration/persistence (e.g. from analyzers), and exporters (e.g. Phy) are not currently able to correctly deal with time-gapped multisegment recordings. Using separate recordings is technically possible, but then there needs to be machinery to handle independently assigned unit IDs across recordings, and many other things.

All this is just to say that, yes, supporting a segment with gaps here is a little bit of an anti-pattern, and in my ideal world we would have a bit more clarity about what a segment actually is and isn't, but if we stop handling segments with gaps correctly in the few places where it is currently needed, that opens up a much wider discussion about how to add support back for these kinds of data.

Footnotes

1. Actually, we don't. Technically, the noninteger-ratio + time vector path can stretch the output timestamps a bit. For example, if someone gives 1000 samples at 1000 Hz with explicit regular timestamps, and asks for a 333.3 Hz output, they get back 333 points spanning 0 to ~0.999 s, which implies a rate slightly less than 333.3 Hz. But the output rate we report back to the user is 333.3 Hz. 🤫 . ↩

[grahamfindlay]

@JoeZiminski Just pushed a batch of commits that hopefully address all of the above. The first 9 commits (through 7d1f053) keep the IIR/FFT filtering/resampling. The 10th and final commit clobbers most of that and replaces it with the FIR/polyphase resampling. So the last commit could be rejected if someone is reluctant to change.