Spike detectors

Spike detectors#

Interictal epileptiform discharge (spike) detectors: Janca (Hilbert-envelope distribution modelling; eeg_forge formulation with presets 'spike' and 'ripple' – the latter not validated on real ripples – and MATLAB-v24 port) and Barkmeier (multichannel half-wave morphology). The detectors are raw (finite input only); GapAwareSpikeDetector wraps any of them for data with NaN/inf gaps. Algorithms, every parameter with its default and source, the filter verification and the comparison with the reference implementations are documented in brainmaze_eeg/spikes/README.md and in the module docstrings below.

Interictal epileptiform discharge (spike) detectors.

Detectors#

  • detect_spikes_janca() – Hilbert-envelope distribution-modelling detector (Janca et al. 2015), eeg_forge formulation with verified filters; presets 'spike' (default) and 'ripple' (80-250 Hz, not validated on real ripples) of one implementation. Layout (n_samples,) or (n_channels, n_samples).

  • SpikeDetectorHilbert (alias spike_detector_hilbert_v24) – port of the MATLAB spike_detector_hilbert_v24 with its full output. Layout (n_samples, n_channels) (MATLAB convention).

  • detect_spikes_barkmeier() – amplitude/slope/ duration half-wave detector (Barkmeier et al. 2012). Layout (n_samples,) or (n_channels, n_samples).

Two layers:

  • raw detectors (above): the algorithm only; finite input required (NaN/inf raise ValueError);

  • GapAwareSpikeDetector: wraps a detector object (JancaDetector, BarkmeierDetector, SpikeDetectorHilbert or any object with a detect(x, fs) method), fills NaN/inf gaps, runs it, removes detections in or near gaps and reports the valid time per channel.

The 2-D entry points raise ValueError on an array with more channels than samples (probably transposed).

See brainmaze_eeg/spikes/README.md for the algorithms, parameters, filter verification and the comparison with the reference implementations.

Janca#

Janca (Hilbert-envelope) interictal spike detectors#

Interictal epileptiform discharge (IED) detection by modelling the distribution of the band-passed signal’s Hilbert envelope as log-normal and flagging envelope maxima that exceed an adaptive threshold derived from that distribution:

Janca, R., Jezdik, P., Cmejla, R., Tomasek, M., Worrell, G.A., Stead, M., Wagenaar, J., Jefferys, J.G.R., Krsek, P., Komarek, V., Jiruska, P., Marusic, P. (2015). Detection of Interictal Epileptiform Discharges Using Signal Envelope Distribution Modelling: Application to Epileptic and Non-Epileptic Intracranial Recordings. Brain Topography 28(1), 172-183. https://doi.org/10.1007/s10548-014-0379-1

Two implementations are provided:

detect_spikes_janca() (primary)

The fast, streamlined formulation of the algorithm from eeg_forge by xnejed07 (eeg_forge/detection/spike_detection_janca.py, xnejed07/eeg_forge), reproduced faithfully (detection indices are identical on the parity data, see the README) with its numerical defects fixed and every parameter configurable. Multichannel layout (n_channels, n_samples). One implementation serves several frequency bands through named presets (see below).

SpikeDetectorHilbert (MATLAB-v24 compatible)

A port of the published MATLAB spike_detector_hilbert_v24.m with its full output (per-detection CDF/PDF weights, multichannel discharge grouping, ambiguous k2 class, segment buffering). Multichannel layout (n_samples, n_channels) (MATLAB’s [samples, channels]).

Presets of detect_spikes_janca() / JancaDetector#

A preset is a complete, named set of parameter values (JANCA_PRESETS); any parameter passed explicitly overrides the preset’s value, e.g. detect_spikes_janca(x, fs, preset='ripple', threshold=4). There is one code path; a preset only supplies defaults. janca_params() returns the resolved, validated values.

'spike' (default)

The eeg_forge reference values: band 10-60 Hz, analysis rate target_fs 200 Hz, 50 Hz notch, 5 s window, threshold 3.65, 0.1 s minimum distance. Source: eeg_forge spike_detection_Janca (commit 54c3704) and Janca et al. (2015).

'ripple'

Ripple band 80-250 Hz (the clinical ripple band, e.g. Zijlmans, M., Jiruska, P., Zelmann, R., Leijten, F.S.S., Jefferys, J.G.R., Gotman, J. (2012). High-frequency oscillations as a new biomarker in epilepsy. Annals of Neurology 71(2), 169-178. https://doi.org/10.1002/ana.22548) analysed at target_fs = 1000 Hz [ours: four times the band top, so that after decimation the band top sits at half the analysis Nyquist; with the default decimation='integer' the analysis rate is fs / floor(fs / 1000), e.g. 1000 Hz for 2-5 kHz input, 1024 Hz for 2048 Hz input]. Input of 501-1999 Hz is not decimated (integer needs fs >= 2 * target_fs) and is analysed at its own rate, where 250 Hz can lie close to Nyquist (0.98 of it at 512 Hz, 0.998 at 501 Hz): the realised edges are still -6.02 dB, but there is almost no room above the band. Input of 500 Hz or less raises. Every other value (threshold, window, minimum distance, notch) is the 'spike' preset’s, carried over untuned. Applying the Janca envelope model to the ripple band is our choice; this preset has not been validated on real ripples (only its filters, resampling and validation are verified by the test-suite). Note that mains harmonics inside 80-250 Hz are not notched by default (notch_harmonics=1 notches only powerline); pass e.g. notch_harmonics=5 to notch 100-250 Hz (50 Hz mains) as well.

Algorithm of detect_spikes_janca()#

For each channel, independently (all filters zero-phase, sosfiltfilt, at the input rate fs):

  1. Band-pass Butterworth, order filter_order (3), edges band (10, 60) Hz. band is the design edge, i.e. the -3 dB point of the single-pass filter; the forward-backward (zero-phase) application squares the response, so the realised response is -6.02 dB at both edges (as in eeg_forge and v24, which also run their designs forward-backward).

  2. Band-stop Butterworth, order notch_order (3), powerline +/- notch_width/2 (50 +/- 2.5 Hz; -3 dB single pass, -6 dB zero-phase at those edges), optionally also at harmonics (notch_harmonics).

  3. Resample with scipy.signal.resample_poly() (its anti-alias FIR runs on the already band-limited signal). decimation='integer' (default, reference): if target_fs is set and fs >= 2 * target_fs, decimate by the integer factor q = floor(fs / target_fs); the analysis rate fs_a = fs / q is then generally not target_fs (500 Hz -> 250 Hz; 512 -> 256; 2048 -> 204.8; 256 or 399 Hz -> not decimated). decimation='exact' (ours): resample to target_fs whenever fs > target_fs by the rational factor up / down (see janca_resampling(); within a relative 1e-6 of target_fs). The resampler’s anti-alias filter attenuates the top ~15 % below the analysis Nyquist, so a resampled configuration is only accepted if that filter loses at most MAX_RESAMPLER_LOSS_DB (0.1 dB) at band[1] (about band[1] <= 0.85 * fs_a / 2); the realised high edge is then -6.02 to -6.12 dB (-6.12 dB at the limit, band[1] = 0.858 * fs_a / 2).

  4. Envelope e = |hilbert(x)|. When the analysis length has a prime factor > 1000 the FFT is padded to scipy.fft.next_fast_len() (pocketfft is 3-11x slower on such lengths; measured). Padding changes the envelope near both ends of the record (the zeros wrap the end into the start): on band-passed noise, > 1 % relative change in the first ~0.2-0.27 s and the last ~0.05-0.2 s. Detections there can differ (on a 6.8 h record cut to a prime analysis length: 5 of ~1500, all in the last 1.6 s; none changed on a 1 h prime-length record), as they would for a record a few samples longer. All other lengths are transformed unpadded, exactly as in the reference.

  5. Sliding statistics of L = log(e + eps) over a centred window of W = int(window_s * fs_a) samples (made odd), mode='reflect' at the ends:

    mu[n] = mean_{k in win(n)} L[k]
sd[n] = sqrt( mean_{k in win(n)} (L[k] - mu[k])**2 )

    Note that sd subtracts the per-sample local mean mu[k] inside the window, not the window-centre mean mu[n]; this is the reference definition and is kept as is (it is not the textbook moving standard deviation; for a stationary background the two agree to first order).

  6. Log-normal mode and median: mode = exp(mu - sd**2), median = exp(mu); threshold T = threshold * (mode + median).

  7. Detections are the maxima of e found by scipy.signal.find_peaks() with height=T and distance=int(min_distance_s * fs_a) samples; returned as sample indices of the input signal (index_a * q, resolution q input samples; with decimation='exact': round(index_a * down / up), i.e. mapped back with the realised ratio, so the rational approximation causes no timing drift).

Channels are processed one at a time, so the working memory is a few times one channel, not a few times the whole montage.

Known differences from the eeg_forge reference (all deliberate fixes)#

  • Filters in ``sos`` form. The reference designs b, a transfer functions; at high sampling rates these lose precision (max complex response error of the 50 Hz band-stop 5e-7 at 2 kHz, 5e-5 at 5 kHz, 2e-3 at 8 kHz, 7e-3 at 10 kHz, 6e-2 at 16 kHz) and the band-stop becomes unstable at 32 kHz (pole radius 1.0006) -> NaN -> silently no detections. The sos design has the same response where the b, a one is accurate: on the parity data the detection indices are identical at 200-10000 Hz; at 16 kHz one detection of 23 moves by one analysis sample; at 32 kHz the reference finds nothing.

  • Scale-invariant epsilon. The reference adds an absolute 1e-6 to the envelope before the log. For data in volts (envelope ~1e-5..1e-8) that constant dominates the background, the threshold no longer adapts, and the detector returns nothing. Here eps = eps_rel * median(e) (eps_rel = 1e-6), so the result does not depend on the unit; at uV scale detections are unchanged.

  • Validation. Every parameter is checked (type, finiteness, range) when it is resolved (janca_params(), also at JancaDetector construction); band edges, notch frequencies, the analysis Nyquist, the resampler loss at the band edge, the window length and the record length are checked against fs at call time (ValueError). A power-line notch that does not fit below Nyquist is skipped with a UserWarning (powerline=None disables it silently).

  • Multichannel input (n_channels, n_samples) (the reference is 1-D only).

  • Non-finite input raises ValueError (the reference silently returns zero detections for a channel with a single NaN). For signals with gaps use GapAwareSpikeDetector around JancaDetector.

  • min_distance_s * fs_a < 1 is clamped to 1 sample (the reference would raise).

  • Optional rational resampling (decimation='exact'); the default keeps the reference’s integer decimation so that results are identical to it.

  • Hilbert FFT length padded to a fast length only for lengths with a prime factor > 1000 (step 4); the parity fixtures and the 6.8 h parity recording are not affected.

brainmaze_eeg.spikes.janca.JANCA_PRESETS = mappingproxy({'spike': mappingproxy({'band': (10.0, 60.0), 'filter_order': 3, 'powerline': 50.0, 'notch_width': 5.0, 'notch_order': 3, 'notch_harmonics': 1, 'target_fs': 200.0, 'decimation': 'integer', 'window_s': 5.0, 'threshold': 3.65, 'min_distance_s': 0.1, 'eps_rel': 1e-06}), 'ripple': mappingproxy({'band': (80.0, 250.0), 'filter_order': 3, 'powerline': 50.0, 'notch_width': 5.0, 'notch_order': 3, 'notch_harmonics': 1, 'target_fs': 1000.0, 'decimation': 'integer', 'window_s': 5.0, 'threshold': 3.65, 'min_distance_s': 0.1, 'eps_rel': 1e-06})})#

Named parameter sets of detect_spikes_janca() (read-only; see the module docstring for their sources). 'ripple' differs from 'spike' only in band and target_fs and is not validated on real ripples.

class brainmaze_eeg.spikes.janca.JancaDetector(preset='spike', **overrides)#

detect_spikes_janca() as a detector object (the protocol of GapAwareSpikeDetector).

JancaDetector(preset, **overrides).detect(x, fs) equals detect_spikes_janca(x, fs, preset=preset, **overrides) for 2-D x (n_channels, n_samples): a list with one int64 array of sample indices per channel. All parameters are resolved and validated at construction (janca_params(): names, types, finiteness, ranges); the checks that need the sampling rate (Nyquist, resampler loss, window and record length) run in detect().

Examples:

JancaDetector()                                  # spikes, eeg_forge values
JancaDetector(powerline=60)                      # spikes, North-American mains
JancaDetector('ripple')                          # 80-250 Hz at ~1 kHz (unvalidated)
JancaDetector('ripple', notch_harmonics=5, threshold=4.0)
detect(x, fs)#

Detection sample indices per channel of x (n_channels, n_samples).

brainmaze_eeg.spikes.janca.MAX_RESAMPLER_LOSS_DB = 0.1#

Maximum loss (dB) of the resampler’s anti-alias filter allowed at band[1].

class brainmaze_eeg.spikes.janca.SpikeDetectorHilbert(**kwargs)#

Janca envelope-distribution IED detector, port of MATLAB spike_detector_hilbert_v24.

Use this when you need v24’s full output (detection CDF/PDF weights, multichannel discharge grouping, the ambiguous k2 class) or comparability with MATLAB v24 results; otherwise prefer detect_spikes_janca(). Input layout is (n_samples, n_channels) (MATLAB convention), the opposite of the rest of the package.

Pipeline (v24): resample to decimation Hz -> power-line notch comb -> 1 Hz high-pass (Butterworth order 2) -> per segment: band-pass bandwidth -> Hilbert envelope -> per-window (winsize/noverlap) log-normal MLE, smoothed and interpolated -> threshold k1*(mode+median) - k3*(mean-mode) -> local maxima, poly-spike union -> output.

Filters (all zero-phase, sos; verified by the test-suite):

  • Power-line notch comb: 2nd-order IIR notches at main_hum_freq and harmonics up to 1.1 * bandwidth[1], unit gain away from the notch. v24 uses a fixed pole radius 0.985, which gives a ~1 Hz wide notch at 200 Hz but widens proportionally with the rate (~24 Hz at 5 kHz with decimation=0). Here the radius is 1 - 0.015 * 200 / fs so the width stays ~1 Hz (identical to v24 at 200 Hz).

  • High-pass 1 Hz Butterworth order 2.

  • Band-pass f_type:

    1. Chebyshev-II (default): minimum-order low-pass and high-pass meeting cheb_rp dB (6) max loss at the band edges and cheb_rs dB (60) stop-band attenuation cheb_transition_hz = (5, 10) Hz beyond the low/high edge, i.e. the v24 spec (normalised 0.05/0.1 at 200 Hz) expressed in Hz so it is valid at any rate. Effective zero-phase response: about -12 dB at the band edges (v24’s 6 dB single-pass spec), < -120 dB in the stop bands. Fixed defect: the earlier port passed the pass-band edge as cheby2’s Wn (which is the stop-band edge) and discarded the order-design Wn, shrinking the effective band to ~18-48 Hz (15 Hz attenuated to 0.03, 50 Hz to 0.15 of the input power).

    2. Butterworth order 4 high-pass + order 4 low-pass (-6 dB at the edges).

    3. FIR (firwin, fs/2 taps, odd) high-pass + low-pass (-12 dB at the edges).

    Change from v24: v24 replaced Chebyshev by Butterworth (with a warning) whenever decimation was not 200 Hz, because its Chebyshev spec was normalised to 200 Hz. Here the spec is in Hz and valid at any analysis rate (verified at 200 Hz-32 kHz), so f_type is always used as given. Pass f_type=2 to reproduce v24 at other rates.

Resampling (decimation > 0 and different from fs): rational polyphase (scipy.signal.resample_poly()) by the smallest-denominator ratio up / down within 1e-6 (relative) of decimation / fs, so any real input rate works (TDT 24414.0625 Hz, 511.99 Hz). The analysis then runs at the realised rate fs * up / down (filters designed and positions converted at that rate, so there is no timing drift). The band must lie below the input Nyquist as well, and the resampler’s anti-alias filter may lose at most MAX_RESAMPLER_LOSS_DB at bandwidth[1].

Parameters (defaults follow v24)#

bandwidth(float, float)

Band-pass edges [low, high] in Hz. Default [10, 60].

k1float

Threshold multiplier for obvious spikes. Default 3.65.

k2float

Threshold multiplier for ambiguous spikes, 0 < k2 <= k1 (MATLAB v24 help: “k1 >= k2”). A local maximum above the k2 threshold but not above k1 is reported as ambiguous (con 0.5) only if an obvious detection on any channel lies within the preceding 10 ms ([i - 10 ms, i]; v24 tests the single sample i - 10 ms, which we read as a typo for this window). Default equals k1 (ambiguous class disabled). With k2 < k1 the result of a channel depends on the other channels, so channel_independent is False and GapAwareSpikeDetector passes the whole montage. (An earlier version enforced k2 >= k1, which inverted the v24 constraint; the ambiguous class could then never fire.) Not comparable to MATLAB v24 when k2 < k1: v23/v25 test the literal single sample, and on real 15-channel iEEG the window used here gives 4-60x more ambiguous detections (a symmetric +-10 ms window 12-50 % more still). See the README (“Ambiguous class is not comparable to v24”). Defaults are unaffected.

k3float

Threshold tilt term. Default 0.

main_hum_freqfloat or None

Mains frequency to notch (Hz). Default 50 (use 60 for North-American data). None disables the notch.

decimationfloat

Target analysis sampling rate (Hz). Default 200. Set 0 to keep the input rate.

bufferingfloat

Core window length for batch processing (seconds). Default 300.

winsize, noverlapfloat

Envelope-model window and overlap in seconds. Defaults 5 and 4. A winsize longer than the analysis record issues a UserWarning (one fit for the whole record).

polyspike_union_timefloat

Poly-spike union interval (seconds). Default 0.12.

discharge_tolfloat

Grouping tolerance for multichannel events (seconds). Default 0.005.

f_typeint

Band-pass family: 1 = Chebyshev-II (default), 2 = Butterworth, 3 = FIR.

cheb_rp, cheb_rsfloat

Chebyshev-II max pass-band loss / min stop-band attenuation, single pass, dB. Defaults 6 and 60 (v24).

cheb_transition_hz(float, float)

Chebyshev-II transition widths (Hz) below the low and above the high edge. Default (5, 10) (= v24 at 200 Hz).

betafloat

Low edge (Hz) of beta rejection. inf (default) disables it; any finite value raises NotImplementedError.

Input must be finite (NaN/inf raise ValueError). For data with gaps use GapAwareSpikeDetector(SpikeDetectorHilbert(...)): it calls detect(), which takes (n_channels, n_samples) and returns detection sample indices per channel.

Buffering scheme#

The record is partitioned into contiguous core windows that tile [0, N) exactly. Each core is analysed inside a block extended by margin = 3 * winsize samples on each side, and only detections inside the core are kept, so every detection is produced exactly once (verified by a whole-vs-buffered test).

Not implemented: beta/mu-activity rejection and the ti_switch == 2 timing mode.

property channel_independent#

ambiguous detections are accepted only next to an obvious detection on any channel.

Type:

True unless the ambiguous class is enabled (k2 < k1)

design_filters(fs)#

Design the filters applied at the analysis rate fs (i.e. after decimation).

Returns:

'notches': list of (centre_hz, sos); 'highpass': 1 Hz sos; 'bandpass': list of filters applied in sequence, each sos (IIR) or ('fir', taps); 'f_type': family actually used.

Return type:

dict

detect(x, fs)#

Detector protocol of GapAwareSpikeDetector.

Parameters:

  • x (np.ndarray, shape (n_channels, n_samples)) – Channels first (unlike run()).

  • fs (float)

Returns:

Per channel, sorted int64 sample indices (input rate) of the detections (obvious and, with k2 < k1, ambiguous ones), round(pos * fs).

Return type:

list of np.ndarray

run(d, fs)#

Detect IEDs in d sampled at fs.

Parameters:

  • d (np.ndarray) – Signal, shape (n_samples,) or (n_samples, n_channels).

  • fs (float) – Input sampling frequency (Hz).

Returns:

  • out (dict) – Per-detection arrays: pos (s), dur (s), chan (int), con (1 obvious / 0.5 ambiguous), weight (envelope CDF), pdf.

  • discharges (dict) – Per multichannel event: MV [n_events, n_chan] spike type, MA max envelope above background, MP start position (s), MD duration (s), MW CDF weight, MPDF pdf.

  • d_decim (np.ndarray) – Decimated, hum-notched, 1 Hz high-passed signal (n_samples_dec, n_chan).

  • envelope (np.ndarray) – Hilbert envelope of the band-passed signal (n_samples_dec, n_chan).

  • background (np.ndarray) – Threshold curves (n_samples_dec, n_chan, 2) for k1 and k2.

  • envelope_pdf (np.ndarray) – Log-normal PDF of the envelope (n_samples_dec, n_chan).

brainmaze_eeg.spikes.janca.design_janca_filters(fs, band=(10.0, 60.0), filter_order=3, powerline=50.0, notch_width=5.0, notch_order=3, notch_harmonics=1)#

Design the filters detect_spikes_janca() applies at the input rate fs.

Parameters are those of detect_spikes_janca(). band and the band-stop edges are design edges (-3 dB single pass); the detector applies the filters forward-backward, so the realised response is -6.02 dB there.

Returns:

'bandpass': sos of the band-pass; 'notches': list of (centre_hz, sos) band-stops actually applied; 'skipped_notches': centres (Hz) that did not fit below Nyquist and were skipped (a UserWarning is issued).

Return type:

dict

Raises:

ValueError – Invalid band, order, notch width or harmonic count.

brainmaze_eeg.spikes.janca.detect_spikes_janca(x, fs, *, preset='spike', band=<preset>, filter_order=<preset>, powerline=<preset>, notch_width=<preset>, notch_order=<preset>, notch_harmonics=<preset>, target_fs=<preset>, decimation=<preset>, window_s=<preset>, threshold=<preset>, min_distance_s=<preset>, eps_rel=<preset>, return_details=False)#

Janca envelope-distribution detector (eeg_forge formulation, fixed), for spikes or, with preset='ripple', the ripple band.

See the module docstring for the algorithm, the presets and the differences from the reference. Every parameter left at <preset> takes the value of preset; any value given explicitly overrides it. The defaults below are those of preset='spike'.

Parameters:

  • x (np.ndarray) – Signal, (n_samples,) or (n_channels, n_samples). Any amplitude unit (the detector is scale-invariant). Must be finite: NaN/inf raise ValueError (wrap the detector in GapAwareSpikeDetector for data with gaps).

  • fs (float) – Sampling frequency of x (Hz).

  • preset ({'spike', 'ripple'}) – Named parameter set (JANCA_PRESETS). 'spike' (default): the eeg_forge reference values. 'ripple': band 80-250 Hz, target_fs 1000 Hz, everything else as 'spike'; not validated on real ripples.

  • band ((float, float)) – Band-pass design edges in Hz (spike: 10, 60; ripple: 80, 250). Single-pass -3 dB; realised (zero-phase) -6.02 dB at both edges. Must satisfy 0 < low < high, high < Nyquist of fs, and, when resampling, lose at most MAX_RESAMPLER_LOSS_DB in the resampler (about high <= 0.85 * fs_a / 2).

  • filter_order (int) – Butterworth prototype order of the band-pass (reference: 3; 1-20).

  • powerline (float or None) – Power-line frequency in Hz (reference: 50). Set 60 for North-American data. None disables the notch.

  • notch_width (float) – Total width (Hz) of the band-stop, centred on powerline (reference: 5, i.e. +/-2.5 Hz design edges; -6 dB at the edges after zero-phase filtering). Must be < 2 * powerline.

  • notch_order (int) – Butterworth prototype order of the band-stop (reference: 3; 1-20).

  • notch_harmonics (int) – Number of notches at k * powerline, k = 1..notch_harmonics (reference: 1). Notches not fitting below Nyquist are skipped with a warning.

  • target_fs (float or None) – Decimation target (Hz) (spike: 200; ripple: 1000). None keeps the input rate.

  • decimation ({'integer', 'exact'}) – 'integer' (default, reference): decimate by the integer factor q = floor(fs / target_fs) only when fs >= 2 * target_fs; the analysis rate fs / q is then generally not target_fs (512 -> 256 Hz, 2048 -> 204.8 Hz, 250..399 Hz -> not decimated). 'exact' (ours): rational polyphase resampling to target_fs (within 1e-6 relative) whenever fs > target_fs, so the analysis rate (and with it the effective time resolution of the envelope model) is the same for every input rate. Detections are mapped back to input samples by rounding, with the realised ratio.

  • window_s (float) – Length (s) of the sliding window of the log-envelope statistics (reference w: 5). Must span at least 3 analysis samples. A window longer than the analysis record issues a UserWarning: the statistics then cover the whole (reflected) record, not a local window.

  • threshold (float) – Threshold multiplier on mode + median of the local log-normal model (reference thr: 3.65, the paper’s k1). Finite, > 0.

  • min_distance_s (float) – Minimum distance between detections in seconds (reference: 0.1). Of two close maxima the larger is kept. Finite, >= 0.

  • eps_rel (float) – Envelope offset before the log, relative to the channel’s median envelope (our choice, replaces the reference’s absolute 1e-6; see module docstring). Finite, >= 0.

  • return_details (bool) – Also return per-channel diagnostics (see Returns).

Returns:

  • detections (np.ndarray or list of np.ndarray) – For 1-D input: int64 array of detection sample indices into x (input rate), sorted. For 2-D input: a list with one such array per channel (row of x). Convert to seconds with detections / fs.

  • details (dict or list of dict) – Only with return_details=True (one dict per channel for 2-D input): fs_analysis (Hz), up/down resampling factors (fs_analysis = fs * up / down), envelope and threshold (at the analysis rate; sample i corresponds to input sample i * down / up), filters (output of design_janca_filters()), preset and params (the resolved parameter values).

Raises:

ValueError, TypeError – Invalid parameters (see janca_params()), a band edge at/above the Nyquist frequency or too close to it after resampling, a window shorter than 3 analysis samples, a record too short for the zero-phase filters, a 2-D array with more rows than columns (probably transposed), or a non-finite value (NaN/inf) in x.

brainmaze_eeg.spikes.janca.janca_decimation_factor(fs, target_fs=200.0)#

Integer decimation factor of the reference rule (decimation='integer').

floor(fs / target_fs) when target_fs is set and fs >= 2 * target_fs, otherwise 1 (no decimation). With the default target_fs=200 this is the eeg_forge rule (decimate only when fs >= 400). The resulting analysis rate fs / q is generally not target_fs (512 Hz -> 256 Hz, 2048 Hz -> 204.8 Hz, 399 Hz -> 399 Hz).

brainmaze_eeg.spikes.janca.janca_params(preset='spike', **overrides)#

Resolved, validated parameters of detect_spikes_janca().

Parameters:
Returns:

One value per parameter (band … eps_rel), ready for detect_spikes_janca(x, fs, **params).

Return type:

dict

Raises:

ValueError, TypeError – Unknown preset or parameter, or a value of the wrong type / out of range (NaN and inf are rejected everywhere except target_fs=None / powerline=None, which mean “off”). Checks that need the sampling rate happen in detect_spikes_janca().

brainmaze_eeg.spikes.janca.janca_resampling(fs, target_fs=200.0, decimation='integer')#

Resampling applied by detect_spikes_janca() after filtering.

Parameters:

  • fs (float) – Input sampling rate (Hz).

  • target_fs (float or None) – Decimation target (Hz); None keeps the input rate.

  • decimation ({'integer', 'exact'}) – 'integer' (reference): decimate by q = floor(fs / target_fs) when fs >= 2 * target_fs (see janca_decimation_factor()). 'exact': rational resampling to target_fs whenever fs > target_fs (never upsamples). up / down is the smallest-denominator continued-fraction convergent of target_fs / fs within a relative 1e-6, so any real rate works (TDT 24414.0625 Hz, rates estimated from timestamps such as 511.99 Hz).

Returns:

(up, down, fs_analysis) – scipy.signal.resample_poly() factors and the realised analysis rate fs * up / down (equal to target_fs to within 1e-6 relative with 'exact').

Return type:

(int, int, float)

brainmaze_eeg.spikes.janca.resampler_gain_db(freq, fs, up, down)#

Gain (dB) at freq Hz of the anti-alias FIR that scipy.signal.resample_poly() designs by default for up / down at input rate fs (Kaiser window, beta 5, cutoff at the lower of the two Nyquist frequencies, 20 * max(up, down) + 1 taps; scipy’s own design rule). 0 when no resampling happens.

Barkmeier#

Barkmeier interictal spike detector#

Multichannel amplitude/slope/duration half-wave spike detector after:

Barkmeier, D.T., Shah, A.K., Flanagan, D., Atkinson, M.D., Agarwal, R., Fuerst, D.R., Jafari-Khouzani, K., Loeb, J.A. (2012). High inter-reviewer variability of spike detection on intracranial EEG addressed by an automated multi-channel algorithm. Clinical Neurophysiology 123(6), 1088-1095. https://doi.org/10.1016/j.clinph.2011.09.023 (PMC3277646)

Written from the paper’s Methods (“Spike detection algorithm”); no third-party code. Each step below is marked [paper] when it is specified by the paper and [ours] when the paper leaves it open and this implementation makes a documented choice.

Algorithm#

  1. [paper] One-minute blocks. The record is processed in successive blocks of block_s seconds (60). The scaling factor, the artifact-channel rule and the candidate threshold are computed per block. [ours] Filtering is done once on the whole record (no block-edge transients); a trailing remainder shorter than half a block is merged into the previous block; block_s=None treats the whole record as one block.

  2. [paper] Artifact channels – OFF by default here (opt-in). The paper excludes, in each block, a channel whose average slope is more than 10 standard deviations from the mean slope of the channels. [ours] Average slope = mean |dx/dt| of the input signal over the (valid samples of the) block. No formulation of this rule we tested keeps every real detection and removes realistic artifacts, so the default is artifact_sd=None, artifact_ratio=None (no channel is ever excluded) and two formulations are offered as opt-in. Measured on 1/f background with 60 s blocks at 256 and 1000 Hz (brainmaze-work/scratch/eeg-spikes/r3/v1_*.py; table in the README):

    • Literal paper rule (mean/SD over all channels): the tested channel is part of the SD, so its z-score is bounded by sqrt(n_channels - 1); it cannot fire below 102 channels (flagged nothing in any test).

    • Leave-one-out mean/SD (round 1 of this module): the reference SD of similar channels is tiny, so a genuinely spiking channel is excluded (8 channels, 300 uV IEDs at 1/s on one: 120 of 598 kept at 256 Hz; 1000 uV at 3/s: 0 of 1794 kept), and 3-4 channel noise montages have 3-12 % false flags.

    • Spatial robust rule (artifact_sd, e.g. 10; round 2 of this module): centre = median slope of the usable channels, spread = max(1.4826 * MAD, artifact_rel_floor * median), flagged when |slope - median| > artifact_sd * spread. With the floor 0.2 a channel is flagged whenever its slope exceeds 3x the median channel’s, whatever the cause: in a montage of 12 contacts at 1x and 4 at 3.5-6x amplitude (grey vs white matter, no artifact) the 4 large channels are flagged in every block and a spiking large channel loses all its detections (0 of 221, 181, 154 and 145 kept at 3.5x, 4x, 5x, 6x). It also removes strong IED bursts (2000 uV at 3/s: 541 -> 4 at 256 Hz). Use it only for montages of similar contacts.

    • Self-referenced rule (artifact_ratio, e.g. 3; ours): each channel’s slope is divided by its own median slope over all blocks, and a channel-block is flagged when that ratio exceeds artifact_ratio times the median ratio of the montage in that block. Heterogeneous montages and steady spiking are never flagged (all cases above kept), but it needs >= 3 blocks, misses an artifact present in more than half of the blocks, and also removes strong IED bursts confined to a few blocks (2000 uV at 3/s in 3 of 10 blocks at 256 Hz: 541 -> 4).

    What either rule catches: broadband noise (white noise at 3x the background rms: slope ratio 5.5-6.2) and mains pick-up (300 uV, 60 Hz). What neither catches (mean slope barely changes, ratio 0.03-1.9): slow drifts, EMG bursts, electrode pops, flat stretches with jumps – although EMG, pops and jumps cause false detections (17-94 in a 3-minute stretch). On real 15-channel iEEG (2 x 1 h, 256 Hz, contact slopes 0.14-4.1x the median) no rule flagged anything. Only usable channels count (>= 3 needed, positive median); excluded channel-blocks get no detections and do not enter the scaling, and a UserWarning names them. Both rules may be combined (a channel-block is excluded if either flags it).

  3. [paper] Candidates. Band-pass 20-50 Hz (narrow_band); candidates are local maxima of the rectified narrow-band signal above a threshold of std_coeff (4) standard deviations. [ours: interpretation] The paper’s wording (“four standard deviations of the channel mean amplitude”) is ambiguous; the threshold here is mean + std_coeff * std of the rectified narrow-band signal in the block (on Gaussian noise about 3.2 SD of the narrow-band signal itself). [ours] Filter type/order are not given in the paper: Butterworth, prototype order narrow_order = 2 (the paper’s broad-band design), zero-phase.

  4. [paper] Morphology band and block scaling. Band-pass 1-35 Hz, 2nd-order Butterworth (broad_band, broad_order). All channels of a block are multiplied by one factor that brings the median (across channels) of the channel mean rectified amplitudes to scale (70 uV). [ours] zero-phase application (sosfiltfilt): the paper only says “second order digital Butterworth”. The band edges are design edges (-3 dB single pass); applied forward-backward the attenuation doubles, so the realised response is -6 dB at 1 and 35 Hz.

  5. [paper] For each candidate, the broad-band peak (largest |x| within +/-2 ms [ours]) and the flanking opposite extrema within trough_search (50 ms [ours]) define two half-waves. A spike is accepted iff total amplitude of both half-waves > 600, each half-wave slope > 7 uV/ms and each half-wave duration > 10 ms, in the block-scaled domain (DEFAULT_THRESHOLDS).

  6. [ours] Accepted detections on a channel closer than trough_search are merged (largest total amplitude kept; one discharge produces several narrow-band maxima); then the optional refractory (s) is applied to the merged list (default 0, the paper defines none).

Filters (verified by the test-suite at 200-32000 Hz)#

signal

type

order

zero-phase response

narrow

Butterworth band-pass

2

-6 dB at 20 / 50 Hz, ~0 dB at 30-35 Hz, -15 dB at 60 Hz

broad

Butterworth band-pass

2

-6 dB at 1 / 35 Hz, ~0 dB at 5-10 Hz

Differences from the earlier version of this module#

  • Broad band was 1-80 Hz (2nd-order high-pass + 4th-order low-pass) and its docstring claimed “per the paper”; the paper specifies 1-35 Hz 2nd-order Butterworth. The 80 Hz band admits sharper noise transients: on 1/f noise the false-positive rate drops ~5x with the paper’s band (see the README for numbers).

  • No blocks: one scaling factor and one threshold per channel for the whole record. Now per block_s (paper: one minute).

  • No artifact-channel rule. Now available as two opt-in formulations (see step 1); off by default because none keeps every real detection.

  • The refractory period was applied before the 50 ms merge, so a merge could pick a detection the refractory had already used to suppress its neighbour. Now merge first.

  • A single NaN anywhere silently made the median scaling factor NaN for every channel. Non-finite input now raises ValueError; data with gaps go through GapAwareSpikeDetector (with BarkmeierDetector), which fills the gaps, passes the gap mask as valid so that block statistics use only real samples, and drops detections in/near gaps.

  • A transposed (n_samples, n_channels) array was silently accepted. A 2-D input with more rows than columns, or shorter than 1 s, now raises ValueError.

Note on false positives#

Block scaling normalises the median channel to scale whatever the channel count, so the fixed thresholds are relative to the typical channel of the montage. On pure 1/f background the detector fires at a low but non-zero rate (0.015-0.027/s per channel with the paper’s 1-35 Hz band; 8 channels of 30 uV 1/f noise) independently of the number of channels; the multichannel design makes the result comparable between channels, it does not by itself remove noise detections. The earlier statement that this was a single-channel artefact was wrong.

class brainmaze_eeg.spikes.barkmeier.BarkmeierDetector(**params)#

detect_spikes_barkmeier() as a detector object (the protocol of GapAwareSpikeDetector).

BarkmeierDetector(**params).detect(x, fs) runs detect_spikes_barkmeier() on the whole montage x (n_channels, n_samples) and returns a list with one list of detection dicts per channel (same dicts as detect_spikes_barkmeier(), sorted by peak_index). Parameters are those of detect_spikes_barkmeier() except valid and return_info; they are validated at construction (names, types, finiteness, ranges); the checks that need fs (band vs Nyquist, trough_search of at least 2 samples) run in detect().

detect(x, fs, valid=None)#

Per-channel lists of detection dicts for x (n_channels, n_samples).

brainmaze_eeg.spikes.barkmeier.design_barkmeier_filters(fs, narrow_band=(20.0, 50.0), broad_band=(1.0, 35.0), narrow_order=2, broad_order=2)#

Design the band-pass filters used by detect_spikes_barkmeier().

Returns:

{'narrow': sos, 'broad': sos} (Butterworth band-passes; applied zero-phase).

Return type:

dict

Raises:

ValueError – Unless 0 < low < high < fs/2 for both bands and the orders are positive integers.

brainmaze_eeg.spikes.barkmeier.detect_spikes_barkmeier(sig, fs, scale=70.0, std_coeff=4.0, trough_search=0.05, thresholds=None, narrow_band=(20.0, 50.0), broad_band=(1.0, 35.0), refractory=0.0, *, narrow_order=2, broad_order=2, block_s=60.0, artifact_sd=None, artifact_rel_floor=0.2, artifact_ratio=None, valid=None, return_info=False)#

Detect interictal spikes with the Barkmeier (2012) multichannel half-wave criteria.

See the module docstring for the algorithm and which choices are the paper’s.

Parameters:

  • sig (np.ndarray) – iEEG in uV, (n_samples,) or (n_channels, n_samples) – pass the whole montage: scaling and the artifact rule are across channels. Must be finite: NaN/inf raise ValueError (use GapAwareSpikeDetector for data with gaps).

  • fs (float) – Sampling frequency in Hz.

  • scale (float) – Block-scaling target for the median channel mean rectified amplitude (paper: 70 uV). Finite, > 0.

  • std_coeff (float) – Candidate threshold in SDs of the rectified narrow-band signal (paper: 4). Finite, >= 0.

  • trough_search (float) – Half-window (s) each side of the peak in which to find the flanking troughs, and merge distance of detections (ours: 0.05). Finite, > 0 and at least 2 samples.

  • thresholds (dict, optional) – {'total_amp', 'slope', 'half_dur'} in the block-scaled domain; partial dicts are completed from DEFAULT_THRESHOLDS (paper: 600 uV, 7 uV/ms = 7000 uV/s, 10 ms). Each finite, >= 0.

  • narrow_band ((float, float)) – Candidate band design edges (paper: 20-50 Hz); -6 dB realised (zero-phase).

  • broad_band ((float, float)) – Morphology/scaling band design edges (paper: 1-35 Hz); -6 dB realised.

  • refractory (float) – Minimum time (s) between accepted spikes on a channel, after merging (default 0). Finite, >= 0.

  • narrow_order (int) – Butterworth prototype orders (broad: paper 2; narrow: ours 2). 1-20.

  • broad_order (int) – Butterworth prototype orders (broad: paper 2; narrow: ours 2). 1-20.

  • block_s (float or None) – Block length in seconds (paper: 60), finite and >= 1. None: one block for the whole record.

  • artifact_sd (float or None) – Opt-in spatial robust artifact-channel rule (paper: 10 SD; default None = off), finite and > 0. Flags any channel whose block slope exceeds 1 + artifact_sd * artifact_rel_floor (3x) times the median channel’s whatever the cause, so large normal channels of a heterogeneous montage lose all detections, and frequent large spikes alone can exceed the limit (the slope ratio depends on the sampling rate: 1000 uV IEDs at 3/s, 256 Hz: 510 of 1235 kept); use only for montages of similar contacts. See the module docstring, step 1.

  • artifact_rel_floor (float) – Floor of the spatial rule’s spread, relative to the median channel slope (ours: 0.2). Finite, >= 0.

  • artifact_ratio (float or None) – Opt-in self-referenced artifact rule (ours; e.g. 3; default None = off), finite and > 1: flags a channel-block whose slope, relative to the channel’s own median over blocks, exceeds artifact_ratio times the montage’s median relative slope in that block. Safe for heterogeneous montages; needs >= 3 blocks; misses artifacts present in most blocks and removes strong IED bursts confined to a few blocks. See the module docstring, step 1.

  • valid (np.ndarray of bool, optional) – Same shape as sig; samples that are real data (default: all). Only valid samples enter the per-block statistics (scaling factor, candidate threshold, artifact slope); a channel with no valid sample in a block is excluded from that block. Used by GapAwareSpikeDetector so that filled gaps do not bias the statistics; detections are not filtered by it.

  • return_info (bool) – Also return a dict with per-block diagnostics (see Returns).

Returns:

  • detections (list of dict) – One dict per spike, sorted by (channel, time): channel, peak_index, peak_time, block, peak_amp, left_amp, left_dur, left_slope, right_amp, right_dur, right_slope, total_amp. Index/time refer to sig (samples / seconds). Amplitudes and slopes are in the block-scaled domain (multiply by 1/info['scale_factor'][block] for input units); durations are seconds.

  • info (dict) – Only with return_info=True: blocks (n_blocks, 2) [start, stop) samples; scale_factor (n_blocks,); artifact (n_blocks, n_channels) bool; channel_slope and candidate_threshold (n_blocks, n_channels) in input units; filters.

Raises:

ValueError, TypeError – Invalid parameter (type, NaN/inf, range), band at/above Nyquist, input not 1-D/2-D, a 2-D input with more rows than columns (probably transposed), a record shorter than 1 s, a non-finite value, or a valid mask of the wrong shape.

Gap handling#

Gap-aware wrapper for the spike detectors#

The detectors in this package are raw: they run their algorithm on a finite signal and raise ValueError on NaN/inf. GapAwareSpikeDetector is the layer for signals with missing data. Around any detector that follows the protocol below, it

  1. finds the gaps of every channel of the original signal: runs of NaN/inf and, by default, runs of an exactly constant value lasting at least flat_as_gap_s (0.1 s; see “Missing data stored as a constant” below);

  2. fills them – gaps up to short_gap_s (0.02 s) by linear interpolation, longer ones with fill ('mirror', 'pink' or 'linear', see brainmaze_eeg.spikes._gaps) – so that filters and background estimates run on a signal without steps or silent stretches;

  3. runs the detector on the filled signal (channels that are entirely missing are left out, so they cannot bias cross-channel statistics such as Barkmeier’s median scaling; detectors that declare channel_independent = True are fed one channel at a time to bound memory);

  4. removes every detection inside a gap or within edge_margin_s (0.2 s) of it – for short, interpolated gaps as well as long ones;

  5. reports the gaps and the valid time per channel, so rates can be normalised by the time in which a detection could have been reported (a RuntimeWarning names channels with less than half of the record valid).

On a gap-free signal the result is identical to calling the detector directly.

Missing data stored as a constant#

Many recordings do not store missing data as NaN: MEF/EDF exports and acquisition systems often write a constant (0, the last value, or a fixed code) during dropouts. A raw detector sees steps at both ends of such a flat run and fires trains of false detections there (on the 6.8 h eeg_forge parity recording, constant runs inside its data_present == 0 intervals produced ~70-85 of the 1494 detections). Convert missing data to NaN when you know where it is (e.g. x[~data_present] = np.nan). As a safety net the wrapper treats any run of exactly equal consecutive samples lasting at least flat_as_gap_s seconds as a gap (default 0.1 s; None disables it). Evidence for the default: in that recording every constant run of >= 0.1 s is missing data, and constant runs within real signal last at most 4 ms (2 samples). Saturated (clipped) stretches are caught too, which is intended. A channel that is constant throughout becomes an all-missing channel (flagged, no detections).

Short gaps and the margin (trade-off)#

Linear interpolation is right for very short dropouts but carries no band power, so with many longer interpolated gaps a background-modelling detector (Janca) sees a lower background and a lower threshold everywhere: on 30 min of real Fz-Cz with 100 ms dropouts every 1 s, the old default short_gap_s=0.1 gave 173 detections in valid time vs 108 without gaps (+67). With short_gap_s=0.02 (gaps > 20 ms mirrored), the same probe gives 110 vs 108, and every tested dropout pattern (4-200 ms, every 0.1-2 s) stays within a few detections of the gap-free run (brainmaze-work/scratch/eeg-spikes/r2/). Even 1-2 sample gaps move nearby detections by > 20 ms, so the margin applies to every gap. The price is valid time: with edge_margin_s=0.2 a gap every 0.4 s or less leaves no valid time at all. That is reported (info['valid_s'], warning) rather than hidden; lower edge_margin_s only if you accept detections influenced by the fill.

Detector protocol#

An object with a method detect(x, fs) where x is a finite float64 array of shape (n_channels, n_samples), returning a list of length n_channels; element c is either

  • a 1-D integer array of detection sample indices into x[c], or

  • a list of dicts, each with an integer 'peak_index' (sample index); a 'channel' key, if present, is rewritten to the channel index of the caller’s array. Such detectors set the class attribute output = 'records' so that empty channels keep the list type (output = 'indices' otherwise).

If the object has accepts_valid = True, detect is called as detect(x, fs, valid=mask) with a boolean (n_channels, n_samples) mask that is False on filled samples, so the detector can exclude them from its statistics. If it has channel_independent = True (its result for one channel never depends on the others), the wrapper calls it once per channel with a (1, n_samples) array.

Implementations: JancaDetector (also with preset='ripple'), BarkmeierDetector and SpikeDetectorHilbert (its detect method).

Example

>>> from brainmaze_eeg.spikes import GapAwareSpikeDetector, JancaDetector
>>> det = GapAwareSpikeDetector(JancaDetector(powerline=60))
>>> spikes, info = det.detect(x, fs, return_info=True)    # x: (n_channels, n_samples)
>>> rate_per_min = [len(s) / (v / 60) for s, v in zip(spikes, info['valid_s'])]
class brainmaze_eeg.spikes.gap_aware.GapAwareSpikeDetector(detector, detector_kwargs=None, *, short_gap_s=0.02, fill='mirror', edge_margin_s=0.2, flat_as_gap_s=0.1, seed=0, fill_kwargs=None)#

Run a spike detector on a signal with gaps (NaN/inf, and constant runs); see the module docstring.

Parameters:

  • detector (object or class) – A detector instance following the protocol (e.g. JancaDetector(), JancaDetector('ripple'), BarkmeierDetector()), or a detector class, instantiated with detector_kwargs.

  • detector_kwargs (dict, optional) – Constructor arguments when detector is a class.

  • short_gap_s (float) – Gaps up to this length (s) are linearly interpolated (default 0.02; see the module docstring for why not longer). Finite, >= 0.

  • fill ({'mirror', 'pink', 'linear'}) – Fill of longer gaps (default 'mirror'; see the README for the measurements behind this choice). Always passed explicitly to the fill function.

  • edge_margin_s (float) – Detections within this distance (s) of a gap, or inside it, are removed (default 0.2: covers the mirror fill’s influence on peak selection, measured at 0.10-0.16 s from the edge, and min_distance_s + half a spike). Finite, >= 0.

  • flat_as_gap_s (float or None) – Runs of exactly equal consecutive samples lasting at least this long (s) are treated as gaps (default 0.1; None disables). See “Missing data stored as a constant”.

  • seed (None, int, sequence of int, np.random.SeedSequence or np.random.Generator) – Seed of the 'pink' fill noise (default 0: reproducible). Validated at construction.

  • fill_kwargs (dict, optional) – Further options of brainmaze_eeg.spikes._gaps.fill_gaps() (context_s > 0, taper_s >= 0, beta; all finite). Validated at construction. The wrapper’s context_s default is 10 s (passed explicitly; _gaps.fill_gaps itself defaults to 30 s like brainmaze_utils.gaps).

detect(x, fs, return_info=False, return_mask=False)#

Detect spikes in x with gaps handled.

Parameters:

  • x (np.ndarray) – (n_samples,) or (n_channels, n_samples); NaN/inf (and constant runs of at least flat_as_gap_s) mark missing data. Not modified; any real dtype (float32 is converted one channel at a time).

  • fs (float) – Sampling frequency (Hz).

  • return_info (bool) – Also return a dict (see Returns).

  • return_mask (bool) – Include info['gap_mask'] (bool, shape of x; True on missing samples).

Returns:

  • detections – The detector’s output with in/near-gap detections removed: for 2-D input a list with one entry per channel (int64 sample-index array, or list of dicts), for 1-D input that single entry. Entirely missing channels give an empty entry.

  • info (dict) – Only with return_info=True. Per channel: gaps ((n_gaps, 2) [start, stop) samples, NaN/inf and flat runs merged), flat_runs (the constant runs treated as gaps), gap_intervals_s (gaps in seconds), all_nan (bool: no usable sample), n_removed (detections dropped in/near gaps), valid_s (seconds in which a detection can be reported: record length minus gaps widened by edge_margin_s; normalise rates by this), valid_fraction; and gap_mask with return_mask=True.

Gap (NaN/inf) helpers used by GapAwareSpikeDetector.

Three steps around an unchanged detector:

  1. find_gaps() – runs of non-finite samples (NaN or +/-inf) of the original signal, as [start, stop) sample indices.

  2. fill_gaps() – make the signal finite so filters/FFT/Hilbert can run:

    • gaps up to max_interp_s (default 0.1 s here; the wrapper passes its own short_gap_s, 0.02 s) are linearly interpolated;

    • longer gaps, method='pink': 1/f noise scaled to the robust amplitude (MAD) of the neighbouring context_s seconds, offset to the local level (bridged linearly between the two sides) and cross-faded with a raised-cosine taper of taper_s seconds into the mirror image of the neighbouring signal at each edge (continuous at the edges);

    • longer gaps, method='mirror': the neighbouring signal mirrored into the gap from both sides, the two images cross-faded over the whole gap (keeps the local spectrum). Real events next to the gap are copied into it. Detections inside the gap are removed by the wrapper, but the mirrored copies still compete with real maxima in peak selection (find_peaks(distance=...)) just outside the gap: on spike-dense real data the extra/missing detections sit 0.10-0.16 s from the edge, which is why the wrapper’s edge_margin_s default is 0.2 s (independent review of PR #67, R7);

    • method='linear': straight line for every gap (not recommended for gaps > ~0.1 s before a background-modelling detector such as Janca).

    Each gap draws its noise from its own generator, seeded from seed and the gap’s position and neighbouring data, so fills are reproducible, independent across channels and independent of the other gaps.

  3. mask_in_gaps() / drop_in_gaps() – remove detections inside a gap or within margin_s of it. Detections on filled samples are never real, whatever the fill.

Relation to brainmaze_utils.gaps#

The canonical gap helpers of the BrainMaze family are being added to brainmaze-utils (brainmaze_utils.gaps, PR bnelair/brainmaze-utils#26, to ship in brainmaze-utils 3.0.0). This module is a deliberately thin, self-contained stand-in with the same names, signatures and semantics as that module’s final API: find_gaps, gap_intervals, fill_gaps(x, fs, *, max_interp_s, method, context_s, taper_s, beta, seed, axis, all_nan, copy), pink_noise, and mask_in_gaps/drop_in_gaps(det, gaps, fs, *, units, gap_units, margin_s=0.1, end=None) with explicit, separate units for the detections and the gaps (both required). The one difference: there is no 'spectral' fill here (the default there), so this module’s default method is 'mirror'. GapAwareSpikeDetector always passes method explicitly, so switching to from brainmaze_utils.gaps import ... is a one-line change (follow-up, once brainmaze-utils 3.0.0 is released; then benchmark 'spectral' as the wrapper default).

brainmaze_eeg.spikes._gaps.drop_in_gaps(det, gaps, fs, *, units, gap_units, margin_s=0.1, end=None)#

det without the detections flagged by mask_in_gaps() (original dtype), or (det, end) when end is given.

brainmaze_eeg.spikes._gaps.fill_gaps(x, fs, *, max_interp_s=0.1, method='mirror', context_s=30.0, taper_s=0.5, beta=1.0, seed=0, axis=-1, all_nan='keep', copy=True)#

Fill the non-finite gaps of a signal (see module docstring).

Parameters:

  • x (array_like) – Signal, 1-D (n_samples,) or N-D with time along axis; every other index is filled independently (a row of an N-D array is filled exactly as the same 1-D signal).

  • fs (float) – Sampling frequency (Hz).

  • max_interp_s (float) – Gaps up to this length (s) are linearly interpolated (default 0.1 s).

  • method ({'mirror', 'pink', 'linear'}) – Fill of longer gaps (default 'mirror'; brainmaze_utils.gaps additionally has 'spectral', its default).

  • context_s (float) – Seconds of valid data on each side used for amplitude, level and mirroring (default 30, as in brainmaze_utils.gaps; the wrapper passes 10 explicitly, its tuned value).

  • taper_s (float) – 'pink': raised-cosine cross-fade length (s) into the mirrored signal at each edge, capped at half the gap (default 0.5).

  • beta (float) – 'pink': spectral exponent, P(f) ~ 1/f**beta (default 1).

  • seed (None, int, sequence of int, np.random.SeedSequence or np.random.Generator) – Noise seed (default 0: reproducible). Each gap gets its own stream derived from the seed, the gap’s position and a hash of its neighbouring data, so channels with the same seed still get independent noise.

  • axis (int) – Time axis (default -1).

  • all_nan ({'keep', 'zero', 'raise'}) – A channel without any finite sample: left unchanged ('keep'), set to zeros ('zero') – both with one warning listing the channels – or ValueError ('raise').

  • copy (bool) – True (default): x is not modified. False: a writable floating-point ndarray is filled in place and returned.

Returns:

Same shape as x; floating input keeps its dtype, anything else becomes float64.

Return type:

np.ndarray

brainmaze_eeg.spikes._gaps.find_gaps(x)#

Runs of non-finite samples (NaN, +inf, -inf) in a 1-D signal.

Returns:

[start, stop) sample indices of each run (stop exclusive), in order.

Return type:

np.ndarray, shape (n_gaps, 2), int64

brainmaze_eeg.spikes._gaps.gap_intervals(x, fs)#

Gaps of a 1-D signal as [start, stop) times in seconds, shape (n_gaps, 2).

brainmaze_eeg.spikes._gaps.mask_in_gaps(det, gaps, fs, *, units, gap_units, margin_s=0.1, end=None)#

Boolean mask of detections overlapping a gap widened by margin_s on each side.

Parameters:

  • det (array_like) – Detection positions, or interval starts, in units.

  • gaps (array_like, shape (n_gaps, 2)) – [start, stop) of each gap in gap_units: find_gaps() gives samples, gap_intervals() gives seconds.

  • fs (float) – Sampling frequency (Hz), required.

  • units ({'seconds', 'samples'}) – Units of det (and end); required, never inferred.

  • gap_units ({'seconds', 'samples'}) – Units of gaps; required, never inferred. Integer-typed values with 'seconds' raise ValueError (pass whole seconds as floats); with 'samples', float values must be integral within 1e-6 (t * fs is fine) or ValueError is raised (a units mix-up must not silently mask nothing).

  • margin_s (float) – Exclusion margin in seconds (default 0.1).

  • end (array_like, optional) – Interval ends (same shape as det, end >= det).

Returns:

True where [det, end] overlaps [gap_start - margin, gap_stop + margin).

Return type:

np.ndarray of bool

brainmaze_eeg.spikes._gaps.pink_noise(n, beta=1.0, fmin_bins=1, rng=None)#

Zero-mean, unit-variance 1/f**beta noise of length n (spectral synthesis). Frequency bins below fmin_bins are zeroed.

Filter design#

Shared filter design and verification helpers for the spike detectors.

Every filter used by a detector is designed here, in second-order-section (sos) form, from parameters given in Hz (never normalised frequencies), and validated against the Nyquist frequency of the rate at which it will run. The sos form is used throughout because the transfer-function (b, a) form of narrow or low-cutoff IIR filters becomes numerically unstable at high sampling rates (e.g. a 3rd-order 47.5-52.5 Hz Butterworth band-stop has a pole outside the unit circle at 32 kHz in b, a form, and its frequency response is already wrong by >1e-3 at 10 kHz), which silently turns the output into NaN or garbage.

All detectors apply their filters forward-backward (scipy.signal.sosfiltfilt()), so the effective magnitude response is |H(f)|**2 (attenuation in dB doubles, phase is zero). zero_phase_response_db() returns exactly that and is what the test-suite uses to verify that every filter does what its parameters say.

brainmaze_eeg.spikes._filters.butter_bandpass(band, fs, order)#

Butterworth band-pass, sos. order is the prototype order (as scipy/MATLAB).

Single-pass response is -3 dB at both edges; forward-backward (zero-phase) -6 dB.

brainmaze_eeg.spikes._filters.butter_bandstop(band, fs, order)#

Butterworth band-stop, sos; -3 dB single-pass (-6 dB zero-phase) at the edges.

brainmaze_eeg.spikes._filters.butter_highpass(cutoff, fs, order)#

Butterworth high-pass, sos; -3 dB single-pass (-6 dB zero-phase) at cutoff.

brainmaze_eeg.spikes._filters.butter_lowpass(cutoff, fs, order)#

Butterworth low-pass, sos; -3 dB single-pass (-6 dB zero-phase) at cutoff.

brainmaze_eeg.spikes._filters.cheby2_highpass(passband_edge, stopband_edge, fs, rp, rs)#

Chebyshev-II high-pass counterpart of cheby2_lowpass() (stopband < passband).

brainmaze_eeg.spikes._filters.cheby2_lowpass(passband_edge, stopband_edge, fs, rp, rs)#

Minimum-order Chebyshev-II low-pass meeting: at most rp dB loss at passband_edge and at least rs dB attenuation from stopband_edge (single pass).

The order and the stop-band edge passed to scipy.signal.cheby2() are the ones returned by scipy.signal.cheb2ord() (cheby2’s Wn is the stop-band edge, not the pass-band edge – passing the pass-band edge shifts the whole filter inwards).

brainmaze_eeg.spikes._filters.check_band(band, fs, name='band')#

Validate a (low, high) band in Hz against the Nyquist frequency of fs.

Raises:

ValueError – Unless 0 < low < high < fs/2.

brainmaze_eeg.spikes._filters.iir_notch_sos(freq, fs, pole_radius)#

Second-order notch with zeros on the unit circle at freq and poles at radius pole_radius (same angle), as one sos row. The -3 dB width (single pass) is about (1 - pole_radius) * fs / pi Hz.

brainmaze_eeg.spikes._filters.zero_phase_response_db(sos, freqs, fs)#

Effective magnitude response (dB) of sos applied forward-backward (scipy.signal.sosfiltfilt()), i.e. 20*log10(|H(f)|**2), at freqs (Hz).

Several cascaded filters can be passed as a list of sos arrays; their responses multiply.