FIR phase-shift + extract apply_raised_cosine_taper from get_chunk_with_margin

[galenlynch]

phase_shift delays each channel by a fraction of a sample, to undo the ADC's staggered sampling. This PR adds method="fir", which applies the same delay with a 16-tap filter instead of an FFT of the whole chunk. On 1 s of 384 int16 channels it takes 27 ms on one thread, against 265 ms for the FFT method. The default stays method="fft".

method="fir" convolves each channel with 16 taps of a sinc

The FFT method's multiply by $e^{i 2\pi f d_c}$ is a convolution with $\mathrm{sinc}(k - d_c)$. That kernel decays quickly, so the FIR method keeps 16 taps of it, Kaiser-windowed (β = 8.6) and scaled to unit DC gain.

• The margin is 8 samples instead of 40 ms, with no taper.
• The kernels run in numba, single-threaded with nogil, and read int16 directly. Without numba, a numpy path gives the same output in 143 ms.
• dtype=None rounds back to the parent's dtype, as the FFT method does; dtype="float32" writes float32.
• n_taps sets the kernel length.

FIR and FFT differ by 34 ppm in the spike band

On 10 s of public Neuropixels 2.0 data (script), the two methods differ by 34 ppm of the signal rms in 300–6,000 Hz and 100 ppm in 300–10,000 Hz. Above 10 kHz the FIR rolls off, where the probe's anti-aliasing filter has already removed the signal. FIR output is identical whether computed in chunks or in one call.

apply_raised_cosine_taper is split out of get_chunk_with_margin

The taper moves into apply_raised_cosine_taper(data, margin, *, inplace=True), so the FIR method can fetch integer chunks without it. window_on_margin=True still works, with a DeprecationWarning, and the FFT method's output is bit-identical.

Tests

• FIR matches FFT within 1% rms on synthetic tones, and recovers a known delay.
• The numpy fallback matches numba for int16 and float32.
• The kernels are in np.convolve order.
• The FFT method is unchanged by the taper refactor.

[galenlynch]

Test failures are unrelated to this PR

[alejoe91]

@samuelgarcia @oliche can you take a look at this implementation?

[oliche]

LGTM

0.2% is not a great relative tol, but this is fine for spike sorting, and this implementation is very valuable for real-time applications !

But I wonder if for waveform extraction we shouldn't push a bit further the analysis. The mismatch with the true FFT implementation will be frequency dependent and worst around Nyquist (which I think is why on a synthetic it is higher than on real data). Here we would reassure neuroscientists by showing waveforms with both techniques applied to see if there is any distortion that could be picked up by downstream models !

Then we could confidently point the users to this method for both spike sorting and waveform extraction.

[oliche]

I don't think we are trying to optimize for sidelobes nor frequency response here, we just want apodization to only modify the edges.

Ideally those are also complementary (ie we could reconstruct the signal by summing batch n-1 with batch n at the edges).

For this the hamming window is not a good choice as it has a non-zero value at its extrema.

I'd suggest a much simpler cosine window, as done in the original IBL implementation !

[galenlynch]

Thanks for taking a look @oliche! I don't understand your last comment, are you referring to the apodization in the FFT phase shift implementation (which I think is what IBL did)?

If it was in response to the window commit, that was a dead-end investigation on my part on the FIR sinc function window — I was wondering why the 0.2% error was so high, and it turned out that it was because I was not normalizing the windowed sinc functions for different phase shifts.

Whether or not 0.2% tolerance matters for biology is an interesting question and caused me to do some investigation and catch the normalization bug and reduce it to 0.0145% while reducing the number of FIR taps to 16 FIR taps. I think this is sufficient for Neuropixels 1.0 and 2.0. The reason is, briefly, because of the limited bandwidth of Neuropixels (0.3 Hz - 10 kHz for NPX 2.0), their input-referred noise, and their limited ADC precision.

To recap the overall approach of the PR: the existing FFT approach to interpolation is equivalent to circular convolution of the signal with a periodic sinc function. This PR speeds that up by convolving with truncated sinc functions, but we first need to window them to prevent Gibbs ripples in the frequency domain and therefore distortion of the interpolated signal. The window function here does matter since Neuropixels have a pass band of 0.3 Hz - 10 kHz: we can use very small sinc functions if our window function does not attenuate the pass band. Hann and cosine windows have 0 at the endpoints and effectively waste two FIR taps as a result.

---

To summarize my investigations: the FIR phase shift's in-band fidelity is ≥12-bit-equivalent (better than the NPX 2.0 ADC itself), its rolloff sits exactly where the analog anti-aliasing filter already rolled off, and its impact on real spike templates is ~200× below the Neuropixel's input-referred noise floor. This makes it safe for waveform extraction and arguably a good choice for everyone, not just real-time pipelines.

The acquisition hardware has a limited passband of useful signal

The NPX 2.0 has a fs=30 kHz and pass band up to 10 kHz (after which the anti-aliasing filter attenuates electrode signal), so we need a window that has a flat pass band to ~0.33·fs and rolls off above. The FIR's rolloff above 10 kHz only attenuates content the analog filter already attenuated. The K=16 Kaiser-windowed sinc achieves just that with the smallest number of taps. We can compare the passband error to discretization noise for various bit depths, where NPX 2.0 has 12 bits, and see that with this Kaiser window K=12 → 7.5 bits, K=16 → 12.8 bits, K=32 only adds ~4 more (past the ADC's own quantization floor).

Filter fidelity in-band (analytical)

Computed from the kernel coefficients via scipy.signal.freqz, averaged over the 16 NPX inter-sample shifts d = k/16:

Band	average ε over shifts	worst-shift (d = ½)	bit-equivalent
0.3 – 6 kHz (spike sub-band)	42 ppm (0.0042 %)	60 ppm	~14.5 bits
0.3 – 10 kHz (NPX 2 passband)	145 ppm (0.0145 %)	205 ppm	~12.8 bits

Worst-case channel in the passband is 205 ppm = 12.3 bits, which is above the 244 ppm = 12-bit threshold of the NPX 2.0 ADC itself.

What the FIR actually does to the recording

The FIR's perturbation is multiplicative at each frequency: whatever sits at frequency f in the recording — spike content, LFP residual, thermal noise, amplifier flicker — gets multiplied by 1 ± ε(f), where ε(f) is exactly the analytical error spectrum |H_FIR(f) − e^{−iωd}| from above. So the FIR's data impact is independent of the recording's spectrum; only the filter's own error matters, and that's known a priori from the kernel.

Per-frequency bounds (NPX-averaged, K=16):

Subband	peak ε(f)	what's there in a real recording
0.3 – 3 kHz	< 80 ppm (0.008 %)	spike content + input-referred noise
3 – 6 kHz	< 100 ppm (0.010 %)	high-f spike tails + input-referred noise
6 – 9 kHz	< 100 ppm (0.010 %)	input-referred noise only
9 – 10 kHz	up to 1761 ppm (0.176 %)	input-referred noise only, near rolloff onset
> 10 kHz	large	already attenuated by the analog AAF

So whatever the spectrum of any individual recording, the FIR preserves each frequency in the spike-content band to within ~100 ppm (0.01 %), and to within ~200 ppm (0.02 %) RMS across the full passband (NPX-averaged 145 ppm = 0.0145 %, worst-shift 205 ppm = 0.0205 %). For reference: 12-bit ADC quantization is 244 ppm (0.0244 %) — the FIR perturbs each frequency by less than the ADC's own discretization noise, over the entire band where signal lives.

The end-to-end validation below — applied to a real recording's actual spectrum — is the empirical confirmation that this per-frequency picture propagates correctly through the rest of the preprocessing chain.

End-to-end on real data

Both pipelines (PhaseShift → Highpass(300 Hz) → CMR(median)) on 10 min of my NPX 2.0 data, with 219 units passing default QC + classifier ≠ noise, median 2 872 spikes per unit. Here are the spectra of the data processed both ways, with the residual:

• Per-spike RMS diff (no averaging): 1.6 µV, flat across all N — ~4× below the 6.8 µV input-referred noise floor (all-in: thermal + amplifier + 1/f + quantization), even on a single spike

• Mean-template RMS diff (after N-spike average): 35 nV — follows 1/√N exactly across 4 decades, ~200× below the input-referred noise floor

• Median template relative error: 0.24 % (99 %ile: 3.6 %, entirely low-N small-amplitude units; no high-SNR unit in the tail)

• Coherent/incoherent decomposition: ~99 % of the residual is incoherent across spikes — broadband noise that averages out as 1/√N, not spike-locked distortion

Mean templates from the two pipelines are visually indistinguishable on every unit I plotted; the 50× amplified residual trace is just noise.

Performance of 16 tap FIR filter

Moving from 32 (main PR text) to 16 taps (this message) doesn't change real-world performance on my machine since it's memory bandwidth limited (dual-channel DDR5) but high performance machines with more memory bandwidth would benefit.

Just for fun, here's a roofline plot of the different FIR taps vs the existing FFT phase shift algorithm:

And the throughput on my machine for the different algorithms themsleves (in reality get_traces adds significant overhead that obscures these performance wins). You can see that delivered performance is capped by memory bandwidth here, despite the arithmetic intensity dropping significantly.

[galenlynch]

@oliche Not sure if I answered your question, or if you were asking a different one?

[oliche]

Yes, at least in the sense that there is no visible distortion on the waveforms themselves and that there is a proper comparison test.

The report here is quite a bit to parse, one last thing I didn't understand is what happens at the seams. Would it make sense to ask the agent to look specifically at seams (ie chunk edges) ?

Other than that, it is quite an optimization !

[galenlynch]

@oliche Sorry if the report was dense. To be clear, even though I did use Claude to make the code changes I thought of the changes and drove the analysis. Perhaps I should ask: does this project have a stance on the use of coding agents? My PRs have received some push back for using them (i.e. "ask the agent about X"). I'm happy to not use them if that would make you and the other maintainers more comfortable.

I worry that we're talking past each other a bit, so I hope you're not offended if I just re-explain the objective of this PR, which will likely be stuff you already know (you're the original author of the existing FFT algorithm I believe). I just want to make sure we're on the same page. Here is a hopefully more clear summary of the code changes and the follow-on analyses that I did to address your concern about relative accuracy of the two phase-shift algorithms.

As you well know, we need to shift different recording channels by a fraction of a sample to account for the multiplexers of Neuropixels sampling inputs in sequence during a sample period. The existing implementation accomplishes this sub-sample shift of a channel's data by taking the Fourier transform of an entire ~30,000 sample recording chunk and multiplying it with a linear phase complex exponential in the frequency domain. That's correct, but it is computationally expensive (~5 * log2N + 15 flops per sample, where N is 30,000 by default in this package, so about 90 flops/sample).

By Whittaker-Shannon, that phase shift operation is exactly equivalent to doing time-domain filtering. This time-domain filtering algorithm is the approach of this PR. The reason to do time-domain filtering instead of frequency-domain multiplication is that it is more computationally efficient for the frequency response we care about. The complex linear phase ramp is exactly equivalent to convolving with an infinite sinc function. Because of finite sample sizes, we can only uses truncated sinc functions. The shorter the truncated sinc, the more we effectively low pass the data, with a sinc the length of the recording chunk being a brick-wall filter to the Nyquist frequency. If we use shorter truncations of the sinc function and apply a window function to avoid Gibbs ringing, we can use a shorter finite impulse response (FIR) filter at the cost of attenuating high-frequency response.

There are two important design constraints here: Neuropixels have a pass-band that goes up to 10 kHz (i.e. even though nyquist is 15 kHz for these probes, the 10-15 kHz band tells us nothing about biological signals), and Neuropixels 2.0 have a input-referred noise floor (from its spec sheet) of 6.8 µVrms in the AP band. This means that the actual signals we care about are below 10 kHz, and there is a real noise floor to bound our digital signal processing decisions. As the section of my previous post, The acquisition hardware has a limited passband of useful signal, shows, a 16 tap FIR sinc filter only attenuates frequencies beyond the pass-band, i.e. it does nothing to biological signals. A 16 tap FIR done in the time domain requires 32 flops/sample, so is around 1/3rd the flop cost of the FFT approach. Even though an efficient algorithm is ~3x faster using naive time-domain FIR filtering, the original benchmarking shows 94x speedups mostly due to RAM traffic etc.

The next section of my previous post, End-to-end on real data, shows that the raw spectrum of my actual data (blue) and the spectrum of my data after going through this 16 tap FIR filter are nearly identical below 10 kHz. The difference between the two spectra (red) is many orders of magnitude below both the stated noise floor of the Neuropixels 2.0 that acquired the data (gray line) and even what is physically plausible to care about - the electrical noise caused by the electrode itself (dashed gray line). This difference becomes larger above 10 kHz but is still frankly small enough to not matter. The average spike plots below show that the spikes extracted with the two phase shift algorithms are visually indistinguishable, and the green difference traces are a) ~four orders of magnitude smaller than the spikes and b) high-frequency ringing that probably doesn't matter for waveform classification. The coherence/incoherence analysis below that just shows that actual spikes have frequency content far below 10 kHz, so the FIR phase shift doesn't affect them.

To answer your question about the seams, I'm not 100% sure I understand what you mean. Are you asking about the edges of the recording chunks? As you know, the FFT approach takes ~40ms (1,200 sample) margins (zero padding at recording boundaries if there's no data), and applies a raised-cosine window to the data to avoid the circular nature of FFTs. The 16-tap FIR sinc requires 8 margin samples to produce valid outputs, and does not require any multiplication of the data. It uses the existing margin machinery of SpikeInterface to get the flanking samples.

---

A potential future benefit of this FIR phase shift approach is that this could be totally eliminated in a future PR. I have a stalled PR that attempts to parallelize some of the other preprocessing steps (#4564). I have yet another PR queued that relies on that one that introduces overlap-save Fourier-domain FIR filtering for high-pass/band-pass filtering as an alternative to IIR filtering. I implemented a similar algorithm for Julia's DSP library and the PR I have queued here shows huge speed improvements over this package's existing time-domain IIR filtering. The very neat thing about this is that, as I've laid out with this PR's characterization, FIR phase shifting necessarily low-passes your data. It turns out that you can combine any sort of frequency response filtering (high pass, band pass, etc) with the phase-shifting for free with frequency-domain FIR filtering. This would mean that instead of having separate phase shift and high pass steps, each of which needs to independently go through each sample of a chunk, you could merge the two kernels and phase shift and high pass in one pass of the data for literally no extra compute per sample.

However, I am not sure if there is much appetite for these numerically nearly-equivalent but vastly more performant changes in this project. If this PR is merged and if #4564 is merged, what would be the right way to propose this overlap-save FIR PR and maybe yet another merged phase shift and filtering PR?

[oliche]

Hi Galen,

No worries, it's more a matter of bandwidth (human bandwidth I mean).

For those implementations, it is good to have a test on RMSE (which you did, plus looking at the waveforms), and also a very specific check for no visible artefacts at the seams. As the filters are different it is always worth checking. It seems this is also done.

For this specific PR I'm all in, the improvement is worth a merge and switch to a default pipeline implementation.

For the PR about designing the filter to do both rephasing and filtering this is also a great idea. I would suggest to have a simpler "convolve" or "time_domain_filter" processing step that takes an arbitrary filter in one of scipy formats, and to perform the filter design separately of its application.

Happy to review the next one !
(one good practice is to post a 10 liner that implements the old and new methods on the same dataset for the reviewer use !)

[galenlynch]

@oliche Your May 15 review is still marked "changes requested", but I think it no longer applies.

That review was on a commit that tried a Hamming window on the FIR's sinc kernel, but I dropped it the same day. The kernel now uses a Kaiser window. The chunk-edge taper that I think you were referring to is unchanged: it is still the IBL raised cosine, zero at the ends. The FIR path needs no edge taper, because a margin of half the kernel length makes each chunk exact.

If that resolves your concern, could you please re-review? Would love to get this merged if there is still interest.

I've rebased onto main, including #4802. Here is the old-vs-new comparison you asked for, on 10 s of public Neuropixels 2.0 data:

import time
import numpy as np
import scipy.signal
import spikeinterface as si
import spikeinterface.preprocessing as spre

url = "s3://aind-open-data/ecephys_786867_2025-09-25_12-43-56/ecephys/ecephys_compressed/experiment1_Record Node 111#Neuropix-PXI-100.50213-1.zarr"
rec = si.read_zarr(url, storage_options={"anon": True}).frame_slice(600 * 30_000, 610 * 30_000)
rec = rec.save(format="memory", progress_bar=False)  # cache so S3 isn't timed

whole, chunked = {}, {}
for method, kw in [("fft", dict(dtype="float32")), ("fir", dict(dtype="float32"))]:
    ps = spre.phase_shift(rec, method=method, **kw)
    whole[method] = ps.get_traces(return_in_uV=True).astype(np.float64)  # one call, no seams
    t0 = time.perf_counter()
    saved = ps.save(format="memory", chunk_duration="1s", n_jobs=1, progress_bar=False)  # 9 seams
    print(f"{method}: {time.perf_counter() - t0:.2f} s")
    chunked[method] = saved.get_traces(return_in_uV=True).astype(np.float64)
    print(f"  chunked vs one call: max |diff| = {np.abs(chunked[method] - whole[method]).max():.3g} uV")

rms = lambda x: np.sqrt(np.mean(x**2))
inner = slice(3000, -3000)  # drop 100 ms at each end
for lo, hi in [(300, 6000), (300, 10000), (10000, 14900)]:
    sos = scipy.signal.butter(8, [lo, hi], btype="band", fs=30_000, output="sos")
    ref, diff = (scipy.signal.sosfiltfilt(sos, x, axis=0)[inner] for x in (whole["fft"], whole["fir"] - whole["fft"]))
    print(f"{lo}-{hi} Hz: rms(FIR - FFT) / rms(FFT) = {rms(diff) / rms(ref):.1e}")

fft: 2.82 s
  chunked vs one call: max |diff| = 0.481 uV
fir: 0.57 s
  chunked vs one call: max |diff| = 0 uV
300-6000 Hz: rms(FIR - FFT) / rms(FFT) = 3.4e-05
300-10000 Hz: rms(FIR - FFT) / rms(FFT) = 1.0e-04
10000-14900 Hz: rms(FIR - FFT) / rms(FFT) = 2.8e-01

• Seams: FIR output is identical whether computed in chunks or in one call.
• Accuracy: FIR and FFT differ by 34 ppm in 300–6000 Hz and 100 ppm in 300–10000 Hz. Above 10 kHz the FIR rolls off, where the probe's anti-aliasing filter has already removed the signal.
• Speed: 4.9× faster through save() with n_jobs=1, and 10× per 1 s get_traces call, both on one thread.

A median common reference downstream spreads the >10 kHz difference into every band (about 1.4% in 300–6000 Hz), because the median is nonlinear. A mean reference does not. That is the source of the ~1.6 µV per-spike difference in my earlier report.

The CI failure is unrelated to this PR.

[galenlynch]

method="fir" now works without numba, through a numpy path that matches the numba output exactly. On 1 s of 384 int16 channels it takes 143 ms, against 29 ms with numba and 272 ms for the FFT method.

[galenlynch]

One more fix: the numba kernels now run single-threaded, like the other preprocessing steps, with the same output and the same speed per thread: 31 ms per 1 s chunk of 384 channels, 9x faster than the FFT method.

They were parallel=True, which runs on whichever threading layer numba finds first, and two of the three break with chunk workers:

layer	pool_engine="thread"	pool_engine="process", fork
tbb (if installed)	ok	ok
omp	ok	workers die (seen after the parent ran the kernel once)
workqueue (numba's fallback)	aborts: it allows one caller at a time	ok

threadpool_limits also doesn't reach numba, so every job used every CPU.

Saving 30 s of 384 channels from RAM, in 1 s chunks. Median of 3 runs.

With several jobs the two match (parallel=True measured with tbb). A single job loses 29% on the phase shift alone (271 to 193 MS/s) and 5% in the full chain.

[galenlynch]

Two small changes, to match fir_filter in #4855:

• output_dtype is gone. dtype="float32" already gives the same output: with an int16 parent, the FIR reads int16 directly and writes float32. With both set, output_dtype silently overrode dtype.
• The kernels are built in fir_filter's tap order, y[n] = sum_j taps[j] * x[n + center - j] with its default center. The kernels are bit-identical to before, and a test checks the order against np.convolve. A follow-up can then hand them to fir_filter, or to a fused step, without reversing them.