Gaps (NaN handling)

Gaps (NaN handling)#

Gap (NaN / inf) handling for detectors: find gaps, fill them before detection, and remove detections that fall in or near them afterwards.

Why fill at all#

Most filters (IIR filtfilt, FFT, Hilbert) spread a single NaN or inf over the whole signal, so a detector run on raw gappy data silently returns nothing. Background-modelling detectors (Janca, RMS/HFO detectors, …) additionally estimate a running background from the signal: a gap filled with zeros or a straight line lowers that background for several seconds around the gap and produces bursts of false detections there; a fill with too much power in the detector’s band raises it and hides real events next to the gap; and a hard step or kink at a gap edge rings through every filter. The fill therefore has to look like the neighbouring background (same level, same power in every band) and must join the data the way the data joins itself.

Fill methods#

  • Short gaps (<= max_interp_s, default 0.1 s), and every gap with method='linear': linear interpolation between the two edge samples (constant extension for a gap at the start or end). All short gaps are filled in one vectorised step, so packet-loss-like recordings with thousands of tiny gaps are cheap.

  • Long gaps, method='spectral' (default): noise with the power spectrum of the neighbouring data, conditioned on the samples at the gap edges (see “Edges”). The spectrum is a Welch estimate (Hann window, half-overlapping segments) from up to context_s of data on each side. The segment is the gap length rounded up to a power of two, but at least ~1 s (so the delta band is in the spectrum) and at most ~8 s (0.12 Hz resolution for long gaps, so a 0.75 Hz slow oscillation is not smeared into 1-4 Hz) and at most half the context. Only segments made of original samples are used when there are at least three (interpolated short gaps and earlier fills have no high-frequency power); with dense packet loss the spectrum above 2 fs / k comes from the longest original-only segments (length k). The estimate is robust to spikes and artifacts in the context: a segment is dropped when its log power exceeds the median over segments by more than max(3 MAD, log 2) in at least two of four log-spaced bands (spikes, pops and movement are broadband), or by more than max(5 MAD, log 16) in any one band; the remaining segments are averaged. Physiological bursts (spindles, alpha/beta bursts) raise one band by less than that and are kept. Every frequency bin of the fill gets exactly the power that the context has there (random phase), so the fill matches the neighbours band by band (the 1/f slope, alpha or other peaks, the white noise floor) rather than only in total RMS. It rides on the local level (median of the 0.5 s next to each edge, linearly bridged across the gap).

  • Long gaps, method='pink': 1/f^``beta`` noise (from 1 / gap up) scaled to the robust (MAD) amplitude of the context, on the same level bridge and with the same edge conditioning. The spectral shape is fixed, so its band powers generally differ from the data (see below), and so does its roughness at the junctions; kept for comparison and for data that really is 1/f^beta.

  • Long gaps, method='mirror': the neighbouring signal mirrored into the gap from both sides, the two images cross-faded over the whole gap with an amplitude correction that accounts for the measured correlation of the two images, plus a 10 ms point reflection blend at each edge. Keeps the local spectrum, but copies real neighbouring events (e.g. a spike next to the gap) into the gap.

Edges#

For 'spectral' and 'pink' the noise is a conditional simulation of the Gaussian process with that spectrum, given ~20 ms (4-64 samples) of data on each side of the gap (kriging with the autocovariance of the same spectrum). The fill therefore continues the data exactly as far as the data’s own autocorrelation reaches: a smooth signal keeps its value and slope, an oscillation keeps its phase into the gap, and white noise carries nothing over (no low-frequency bump from one noisy edge sample). The variance is right everywhere, so the band power neither dips nor bulges at the edges, and the junction is statistically like the data’s own sample-to-sample steps. Edge residuals larger than 4 robust SD of the neighbouring data (e.g. an artifact cut by the dropout) are clipped before conditioning, so an artifact at the edge is not continued into the gap. The same clip applies to a large real event (e.g. a 300 uV K-complex or a 400 uV spike) that the gap cuts: the fill then starts from the clipped value, which leaves a step at the junction (K-complex: junction step 42x the median absolute sample difference, 10-60 Hz envelope 12.6x vs 2.2 gap-free; without the clip 1.2-1.9x). The step is confined to +-0.1 s of the edge, which the default drop_in_gaps margin (margin_s=0.1) covers; do not use a smaller margin if such events occur. taper_s (default 0.5 s) bounds how far into the gap the conditioning acts (the whole gap if it is shorter than 2 * taper_s; faded out over the second half of taper_s). taper_s=0 gives unconditioned noise with a step at each edge.

Units (post-filters)#

mask_in_gaps / drop_in_gaps take two required keywords: units for the detections (and end) and gap_units for the gaps, each 'seconds' or 'samples'; margin_s is always in seconds and fs is always required. Gaps from find_gaps() are sample indices, from gap_intervals() seconds; detections are whatever your detector returns (e.g. eeg_forge’s Janca returns sample indices). Because each side is stated separately, detections in samples work against gaps in seconds (and vice versa) without conversion, and a units mix-up cannot hide behind one shared keyword. Integer-typed values with 'seconds' raise (sample indices declared as seconds; pass genuine whole seconds as floats), and non-integral values with 'samples' raise (floats within 1e-6 of an integer, such as t * fs, are accepted).

Randomness#

seed (default 0) may be an int, a sequence of ints, None (fresh entropy), a numpy.random.SeedSequence or a numpy.random.Generator. The noise of every long gap is drawn from its own stream derived from the seed, the gap’s length and a hash of the context data next to it, quantised to 1/1000 of its robust SD around its median. Consequently:

  • a channel’s fill does not depend on other channels, their gaps or their order, nor on the other gaps of the same channel (as long as they are outside its context);

  • the same channel gives the same fill in a 1-D call and inside an N-D call;

  • the fill is the same across machines and numpy/scipy versions, and after upstream processing that differs at the rounding (ulp) level; for data rescaled (uV vs V) or and for the same gap at another array offset (segment-wise processing), as long as the context is the same. Data that differ by more than about 1/1000 SD give a different, independent realisation. A float32 cast of real data is not guaranteed to give the same fill (the rounding of ~6e-8 relative can cross a 1/1000 SD quantisation boundary: 5 of 20 real segments changed, 20 of 20 were unchanged for uV -> V or 1e-9 relative noise); fill in one dtype throughout if bit-reproducibility across dtypes matters;

  • different channels sharing a gap (recording-wide dropouts) get independent noise, also when they are filled one at a time in separate calls with the same seed, so the fills are not spuriously identical. (Two channels whose context data are equal up to scale get identical fills.) This does not make montages safe, see “Montages” below.

Montages (bipolar, CAR): derive first, then fill#

Derive the montage (bipolar, common average, …) from the raw channels first and fill the derived signal afterwards. Filling the referential channels and then deriving subtracts independent noise realisations from channels that are strongly correlated in the real data (the real common part cancels, the independent fills add), which inflates the gap in the derived signal 5-22x in RMS (two referential channels with correlation 0.962: 5.2x; 0.998: 21.6x) and raises the Janca threshold 1.22x (p90 1.38) up to 2.5 s outside the gap (0.995 correlation). Deriving first and then filling gives 1.01-1.04 (gap RMS 1.01-1.02). Gaps of the derived signal are the union of the channels’ gaps (non-finite samples propagate through the subtraction), so find the gaps on the derived signal.

Defaults and the evidence behind them#

max_interp_s=0.1: in the Janca detection benchmark (1 h scalp EEG, 500 Hz) gaps <= 0.05 s made no difference for any method, while straight lines over 0.5 s gaps already produced false detections; 0.1 s also matches the default post-filter margin.

method='spectral': measured on real Fz-Cz EEG (500 Hz, 20-40 gaps per length) and on synthetic backgrounds (probes and tables in brainmaze-utils PR #26):

  • band RMS of the fill / band RMS of the neighbouring data, median over gaps, 0.5 Hz to 1 kHz: 'spectral' 0.96-1.07 for 2-60 s gaps and 0.84-0.99 for 1 s gaps (real), 0.88-1.04 for 1-60 s gaps at 5 kHz; pooled over 39 gaps of stationary synthetic data 0.81-1.19 also for 0.15-0.5 s gaps. In 0.2-0.5 s gaps the delta band (0.5-4 Hz) of real EEG is under-filled (0.57-0.76; frequencies below ~1/gap cannot be represented). 'pink': 1.4-3.2 in 4-200 Hz on real data and about 5-7 in 80-1000 Hz at 5 kHz, because its 1/f shape is fixed;

  • slow oscillation (0.75 Hz) + white noise, 30 s gaps: 0.6-1 Hz 0.94, 1-2 Hz 1.5, 2-4 Hz 0.99 (with the former ~1 s Welch segment: 0.38 / 6.7 / 7.7);

  • noisy (white-dominated) edges, 5 s gaps: 1-4 Hz power in the first 0.5 s of the fill 1.65x the original, the same as unconditioned noise (the former one-sample point reflection: 181x); 1-4 Hz difference in the valid data 0.1-1 s before the gap 0.50 band-SD median (point reflection: 6.1);

  • a 300 uV artifact cut by the dropout: max |fill - level| in the first 0.5 s 3.7 SD (point reflection: 19.9 SD, 'mirror': 10.7 SD); junctions on real EEG are statistically like the data’s own (2nd difference p90 3.2-3.8 vs 3.4 for the data);

  • an RMS (80-500 Hz, 10 s window) background threshold in the valid data 0.1-4 s outside a gap, filled / gap-free: 'spectral' 1.00 and 'mirror' 0.99-1.00 for 0.5, 2 and 10 s gaps; 'pink' 1.69 / 2.67 / 2.94, i.e. an RMS-threshold detector misses real events next to a pink-filled gap;

  • eeg_forge’s Janca threshold 0.1-2.5 s outside a gap, filled / gap-free (0.2-60 s gaps): 'spectral' median 0.99-1.01, p90 <= 1.09; 'mirror' p90 <= 1.04; 'pink' p90 up to 1.23;

  • Janca detections (1 h, 23 gaps per length, 0.5-60 s, transients injected 0.15-1.2 s outside both edges, drop_in_gaps margin 0.1 s): no false detections near the gaps for 'spectral' or 'mirror' ('linear': 12-415). With strong transients both find 100 %. With weak ones (near the threshold) the sensitivity at 0.15 s from the edge is 0.74-0.91 gap-free, 0.72-0.85 with 'spectral' and 0.41-0.78 with 'mirror' (it copies the transients into the gap, raising the background).

'spectral' is the default because it is the only fill that keeps the background of the neighbouring data in every band without copying neighbouring events (spikes, artifacts) into the gap; its Welch estimate also rejects such events in the context. context_s=30: room for the 8 s Welch segments of long gaps (>= 3 per side; with 10 s the slow-oscillation test above gave 1-2 Hz 3.3x); shorter gaps use about 16 gap lengths (at least three ~1 s segments) per side, so the cost scales with the gap. taper_s=0.5: covers the autocorrelation of typical EEG; the conditioning is exact within it.

Release#

This module is new in brainmaze-utils 3.0.0 (a major release because numerical results of filled signals change); the version itself is bumped by the release workflow.

Limitations#

  • At very high packet loss (>= 75 %, islands of valid data shorter than 16 samples) the fill underestimates the spectrum: 8-30 Hz 0.26-0.83 and 30-240 Hz 0.44-1.03 of the original (10-50 % loss: 0.99-1.00).

  • Dense large slow-oscillation or K-complex trains next to short gaps: their cycles are rejected as artifacts by the Welch rule, so the fill is slightly too weak (1 s gaps: 11-15 Hz 0.83-0.89, 30-60 Hz 0.90-0.92 vs 0.96-1.01 without rejection; 5-30 s gaps 0.93-1.02).

  • The fill is stationary noise: it does not reproduce oscillatory bursts, spikes, sleep spindles or other non-stationary structure, nor the phase relations between channels (each channel is filled independently; derive bipolar/CAR montages before filling). Do not compute features that depend on the content of the gap (event rates, coherence, phase) without excluding the gaps.

  • Bursty rhythms are filled at their average power: with the two-band rejection rule the burst band of the fill is 0.94-1.04 of the context for spindles and alpha/beta bursts (it was 0.62-0.91 with a one-band rule). The price: a single narrow-band artifact below 16x is kept (0.3 s, 600 uV movement artifacts every 5 s next to a 1 s gap: 1-4 Hz 1.24x of the clean background instead of 1.10x; spikes in the context are still rejected: 10-60 Hz 0.98-1.04, 1.7x without rejection).

  • For gaps >= ~4 s the Welch segments are 8 s long; artifacts that recur more often than that are in every segment and cannot be rejected: the movement artifacts above give 1-4 Hz 3.5x next to a 5 s gap (10-60 Hz unaffected, 0.99).

  • The spectrum is estimated from at most context_s on each side; if the state changes across the gap the fill uses the average of both sides. With fewer than 16 valid context samples the fill falls back to 'pink' scaled by the MAD of what is there.

  • Frequencies below about 1 / segment (0.12-1 Hz) are represented only by the level bridge and the edge conditioning; the fill has no infra-slow drift. A strong narrow peak (a 0.75 Hz slow oscillation) still leaks ~1.5x into the neighbouring band.

  • Remove line noise before filling: in the shortest long gaps (<= ~0.15 s) the line in the fill is not phase-locked to the data and passes a notch (40 uV 50 Hz, eeg_forge Janca 10-60 Hz band with a 50 Hz notch: 2.4x in 0.13 s gaps, 0.94-1.04 in 0.2-1 s gaps, where the conditioning keeps the line’s phase).

  • With dense packet loss (no ~1 s stretch of original data) the spectrum below 2 fs / k (see above) comes from segments that include interpolated samples and is biased low.

  • A detector’s background estimate straddling a gap is still partly made of fill, so the post-filter margin should cover the detector’s own edge sensitivity (default 0.1 s; increase it for detectors with long filters or windows).

  • Short gaps are linear: continuous in value but not in slope. A slope-matching (cubic) interpolation would extrapolate sample-to-sample noise over the gap; the kink is covered by the post-filter margin.

  • 'mirror' corrects its cross-fade amplitude with the broadband correlation of the two mirror images; a narrow-band coherent component such as line noise can still come out ~1.1-1.3x too strong in short gaps (it was up to 1.5x without the correction).

  • Cost: a long gap at 32 kHz takes ~0.5 s per minute of gap (10 min: ~5 s) and allocates ~5.6-5.8x the gap length in float64 temporaries (30 min at 32 kHz: ~2.5 GiB); medium gaps (0.1-0.5 s) at 32 kHz take ~25-55 ms each (1500 of 0.12 s in 10 min: ~40 s). Fill very long gaps at a lower rate, or split the recording.

  • A channel that is entirely NaN/inf cannot be filled; see all_nan.

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

Remove the detections flagged by mask_in_gaps() (same arguments).

Returns:

The kept detections (and interval ends), with their original dtype.

Return type:

np.ndarray, or tuple (det, end) when end is given

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

Fill NaN/inf gaps so a detector can run on the signal (see the module docstring for the pipeline, the methods and the evidence behind the defaults).

Parameters:

  • x (array_like) – Signal, 1-D (n_samples,) or N-D with time along axis. Non-finite samples (NaN, +inf, -inf) are gaps.

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

  • max_interp_s (float) – Gaps up to this length (seconds) are linearly interpolated; longer gaps use method. Default 0.1 s.

  • method ({'spectral', 'pink', 'mirror', 'linear'}) – Fill for long gaps. 'spectral' (default): noise with the band powers (robust Welch PSD) of the neighbouring data. 'pink': 1/f^beta noise at the MAD amplitude of the neighbouring data. 'mirror': neighbouring signal mirrored in from both sides (copies neighbouring events into the gap). 'linear': straight line (not suitable for gaps longer than ~0.1 s before background-modelling detectors).

  • context_s (float) – Maximum seconds of valid data used on each side of a long gap (spectrum, level, amplitude, mirror images). Default 30 s. 'spectral' uses less for shorter gaps (about 16 gap lengths, at least three ~1 s Welch segments).

  • taper_s (float) – 'spectral'/'pink': how far (seconds) into the gap the fill is conditioned on the data at each edge (the whole gap if it is shorter than 2 * taper_s; faded out over the second half). Default 0.5 s. 0 disables the conditioning (the noise then starts with a step; not recommended).

  • beta (float) – Spectral exponent for 'pink' (and for the 'spectral' fallback when fewer than 16 valid context samples exist).

  • seed (None, int, sequence of int, SeedSequence or Generator) – Root of the per-gap random streams (see “Randomness” in the module docstring). Default 0: reproducible. A gap’s noise depends only on the seed, the gap’s length and the data next to it (quantised to 1/1000 of its robust SD), not on its position in the array.

  • axis (int) – Time axis for N-D input; every other index is a channel, filled independently.

  • all_nan ({'keep', 'zero', 'raise'}) – A channel without any finite sample: 'keep' (default) leaves it unchanged (non-finite) and warns once, listing the channel indices; 'zero' fills it with zeros and warns; 'raise' raises ValueError.

  • copy (bool) – If False and x is a writable floating-point ndarray, it is filled in place and returned (no extra memory). Otherwise a filled copy is returned.

Returns:

Same shape as x; same dtype for floating-point input (float32 stays float32), float64 otherwise. Finite everywhere except all-NaN channels with all_nan='keep'. Valid samples are never changed.

Return type:

np.ndarray

Notes

Always remove detections in gaps afterwards with drop_in_gaps() / mask_in_gaps(), using gaps found on the original signal.

brainmaze_utils.gaps.find_gaps(x)#

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

Parameters:

x (array_like, shape (n_samples,))

Returns:

[start, stop) sample indices of each run (stop exclusive), in order. Use them in mask_in_gaps() / drop_in_gaps() with gap_units='samples'.

Return type:

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

brainmaze_utils.gaps.gap_intervals(x, fs)#

Gaps of a 1-D signal as [start, stop) times in seconds, shape (n_gaps, 2), float64. Use them with gap_units='seconds' in mask_in_gaps() / drop_in_gaps().

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

Boolean mask of detections that are in or near a gap (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. Need not be sorted; may overlap.

  • fs (float) – Sampling frequency (Hz). Always required (margin_s is in seconds).

  • units ({'seconds', 'samples'}) – Units of det (and end). Required, no default.

  • gap_units ({'seconds', 'samples'}) – Units of gaps. Required, no default. The two are separate on purpose: a mix-up (e.g. sample-index detections against gaps converted to seconds) would silently mask nothing, so you state each one. 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.

  • margin_s (float) – Exclusion margin around every gap, in seconds, >= 0. Default 0.1 s.

  • end (array_like, optional) – Interval ends (same units and shape, end >= det) for interval detections; an interval is masked if any part of [det, end] overlaps a widened gap.

Return type:

np.ndarray of bool, shape of np.atleast_1d(det); True = drop the detection.

brainmaze_utils.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).

Parameters:

  • n (int) – Number of samples.

  • beta (float) – Spectral exponent of the power spectrum, P(f) ~ 1/f**beta (1 = pink, 0 = white, 2 = brown).

  • fmin_bins (int) – Lowest non-zero frequency bin kept. Bins below are zeroed so that the realisation does not wander on time scales longer than the segment itself.

  • rng (np.random.Generator, int or None) – Generator or seed (passed to numpy.random.default_rng()).

Return type:

np.ndarray, shape (n,), float64