Signal#

Tools for digital signal processing: decimation/resampling, FFT filtering, low-frequency filtering utilizing downsampling, buffering etc.

Conventions

Multichannel arrays are (n_signals, n_samples): time runs along the last axis.
Sampling frequencies and cutoffs are in Hz, durations in seconds.
Functions do not modify their inputs.
NaN marks missing data (gaps); each function documents how it treats NaN.

Which functions filter?

Function	Anti-aliasing / filtering	NaN handling
`decimate()`	yes: 16th-order Butterworth (SOS, zero-phase), default cutoff `fs_new / 3`	interpolated for filtering, re-masked in the output
`nandecimate()`	yes, weak: 30-tap FIR (legacy; prefer `decimate()`)	mean-filled for filtering, re-masked in the output
`resample()`	no - interpolation only, by design; low-pass first when downsampling	never interpolated across; an output is NaN if a bracketing input sample is. When downsampling, a gap shorter than one output period can vanish
`unify_sampling_frequency()`	yes (calls `decimate()`)	as `decimate()`
`downsample_min_max()`	no (min/max envelope, for display)	NaN ignored within a window
`LowFrequencyFilter`	is itself a zero-phase low-/high-pass	short gaps bridged, long gaps split into separately filtered segments; NaN out exactly where NaN in (`inf` raises)
`fft_filter()`	is itself a brick-wall low-/high-pass	as `LowFrequencyFilter`

class brainmaze_utils.signal.LowFrequencyFilter(fs=None, cutoff=None, n_decimate=1, n_order=None, dec_cutoff=0.3, filter_type='lp', ftype='fir', max_gap_fill=None)#

Zero-phase low-pass or high-pass filter for very low cutoff frequencies relative to the sampling rate (e.g. a 0.5 Hz high-pass on 8 kHz data), implemented as a decimate -> filter -> upsample cascade.

A cutoff of 1e-4 * fs cannot be realised well by a direct filter (an FIR would need ~1e5 taps, an IIR is ill-conditioned). Instead the signal is

halved in rate n_decimate times: each step applies a zero-phase anti-aliasing low-pass (normalised cutoff dec_cutoff x Nyquist of that stage) and keeps every 2nd sample;
low-passed at cutoff at the low rate fs / 2**n_decimate (zero-phase);
brought back to fs by n_decimate steps of zero-insertion followed by the same anti-aliasing low-pass (gain 2).

The result of 1-3 is the low-frequency content of the signal. filter_type='lp' returns it; filter_type='hp' returns x minus it.

Parameters:

fs (float) – Sampling frequency of the signals that will be filtered, in Hz.
cutoff (float) – Cutoff frequency in Hz. Must lie below the Nyquist frequency of the low rate, cutoff < fs / 2**n_decimate / 2, and should lie well inside the band kept by the anti-aliasing stages (< dec_cutoff * fs / 2**(n_decimate + 1)).
n_decimate (int) – Number of halving steps (>= 0). The cutoff filter runs at fs / 2**n_decimate. Choose it so that this rate is ~20-500 x cutoff, e.g. n_decimate=5 (250 Hz) for 0.5 Hz on 8 kHz data.
n_order (int, optional) – Filter order: number of taps for ftype='fir' (default 101), Butterworth order for ftype='iir' (default 3). Used for both the anti-aliasing and the cutoff filter. An FIR needs about n_order >= 1.5 * fs_low / cutoff taps to resolve the cutoff (fs_low = fs / 2**n_decimate), otherwise its transition band is much wider than cutoff and the effective cutoff is higher: the default 101 taps for 0.5 Hz at fs_low = 250 Hz pass 0.94 at 0.5 Hz and reach -6 dB only at ~1.7 Hz. A UserWarning is issued when the FIR’s low-frequency gain at cutoff exceeds 0.3 (nominal 0.25).
dec_cutoff (float) – Normalised (to the Nyquist of each stage, 0 < dec_cutoff < 1) cutoff of the anti-aliasing low-pass used by the halving/upsampling steps. Default 0.3.
filter_type ({'lp', 'hp'}) –
Which side of cutoff is returned:
- 'lp' returns the low-frequency content below cutoff.
- 'hp' returns x minus that low-frequency content, i.e. the content above cutoff (use this to remove slow drift / DC).
max_gap_fill (float, optional) – NaN gaps (interior runs of NaN) shorter than this many seconds are bridged by linear interpolation before filtering; longer gaps split the signal into segments that are filtered separately. Default 0.5 / cutoff (1 s for 0.5 Hz). 0 always splits. See Gaps in the notes.
ftype ({'fir', 'iir'}) – 'fir': scipy.signal.firwin() windowed-sinc filters (cutoff is the design -6 dB point). 'iir': Butterworth filters in second-order sections (cutoff is the design -3 dB point). All filters are applied forward and backward, so the phase is zero and the magnitude response is squared: at cutoff the low-frequency gain is ~0.5 (IIR) or ~0.25 (FIR, only if it has enough taps, see n_order).

Examples

import numpy as np
from brainmaze_utils.signal import LowFrequencyFilter

fs = 8000
x = np.random.randn(60 * fs) + 1500.0 + np.linspace(0, 300, 60 * fs)  # DC + drift

# remove DC and slow drift below 0.5 Hz (high-pass), keep everything above
hp = LowFrequencyFilter(fs=fs, cutoff=0.5, n_decimate=5, ftype='iir', n_order=3, filter_type='hp')
x_hp = hp(x)

# keep only the slow content below 0.5 Hz
lp = LowFrequencyFilter(fs=fs, cutoff=0.5, n_decimate=5, ftype='iir', n_order=3, filter_type='lp')
x_slow = lp(x)

Notes

Input. x may be 1-D (n_samples,) or N-D with time along the last axis, e.g. (n_channels, n_samples). NaN marks gaps (see below); inf raises ValueError. The input is not modified.

Gaps (NaN). The output is NaN exactly where the input is NaN, and finite everywhere else, so the filter can follow a NaN-aware decimate() in a downsampling cascade.

Gaps shorter than max_gap_fill (default half a cutoff period, e.g. dropped samples) are bridged by linear interpolation, filtered through, and re-masked. The bridge does not check the two sides for a level jump: a jump across a bridged gap (amplifier re-zero, reconnection after a dropout) is filtered as a genuine step. E.g. a 500 uV DC step across a 0.1 s gap at 0.5 Hz high-pass gives a transient above 10 uV for about +-1.5 s around the gap (the same as the response to a real step). If such jumps are expected, lower max_gap_fill below the gap length (0 always splits) so that each side is filtered on its own.
Longer gaps split the signal: each valid segment is filtered on its own, with the same edge handling as the record edges (below), so there is no jump at a gap edge, and nothing on one side of a long gap affects the other side. A segment shorter than about 1 / cutoff s (e.g. a few samples between two long gaps) carries no information at the cutoff frequency: its output is finite but is not a meaningful low-/high-pass estimate (a 3-sample segment can be off by about 10 for SD-20 data). Treat such segments as unusable.

Measured on 1/f noise (SD 20) with offset, drift and a 10 Hz tone, 0.5 Hz IIR high-pass, maximum error within 10 s of a gap compared with the same signal without the gap: 1 ms gap 0.008, 10 ms 0.14, 100 ms 0.95 (bridged); 1 s 5.4, 10 s 7.5 (split; bridging would give 6.0 and 8.1). The error next to a split gap is the record-edge error described below. Leading and trailing NaN are simply excluded. A multichannel array whose channels share the same gaps is processed in one pass. Every split segment is extended by 3 / cutoff s at both ends, so many long gaps cost time: 10 min at 8 kHz (0.5 Hz IIR) takes 1.3 s without gaps, 1.5 s with 500 dropouts, 2.8 s with 50 two-second gaps and 11 s with 300 split segments.

Record edges. A filter with a 0.5 Hz cutoff has a transient lasting seconds, so what is assumed about the signal beyond the two ends of the record determines the first and last few seconds of the output. This class

removes a least-squares straight line (offset + drift) from each signal and adds it back to the low-frequency output (a zero-phase low-pass passes a straight line unchanged, so for 'hp' offset and linear drift are removed exactly, edges included);
extends each end by max(3 / cutoff s, the legacy pad) samples: a line fitted to the outermost 2 / cutoff s is extrapolated and the residual is mirrored about it, so the extension is continuous in value and slope with the data;
runs the cascade on the extended signal and crops the extension.

For a signal with a 2000 (uV) offset, a 40 uV/s drift and three 10 uV tones (3, 10, 40 Hz) at 8 kHz, a 0.5 Hz 'hp' (n_decimate=5, IIR order 3) gives a maximum edge error of ~0.75 uV over the first/last second (previously ~1900 uV); offset and drift no longer contribute at all (the output equals that for the zero-mean signal). The remaining edge error is the irreducible effect of not knowing the signal beyond the record, of the order of A * cutoff / (pi * f) for a component of amplitude A at frequency f; on 1/f-like EEG it scales with the signal power near and below cutoff. Far from the edges (more than ~3 / cutoff s), the output equals that of the cascade applied to an infinitely long signal.

Design choice: extensions were compared on synthetic data (scripts in the PR). On 1/f-like signals the mirrored extension gives the smallest edge error (zero extension of the detrended signal is ~18 % worse); for a signal whose content lies entirely well above cutoff (pure tones), zero extension would be ~2x better (0.36 vs 0.77 uV in the example above) but only without any offset or drift.

Note

Changed after v2.0.0: The signal used to be zero-padded by only 2 * n_order * 2**n_decimate samples (e.g. 24 ms for a 0.5 Hz filter at 8 kHz). Any offset therefore became a step at both edges (output off by about half the offset, e.g. ~1000 uV for a 2000 uV offset, over the first/last seconds), and the filter transient reached far into the record. The frequency response in the interior is unchanged. IIR filters are now applied in second-order sections (numerically identical response; the b_*/a_* attributes are kept for reference). N-D input is supported. NaN gaps are handled (see Gaps); v2.0.0 returned an all-NaN output for any input containing a NaN. inf raises ValueError. Short FIRs that cannot resolve the cutoff now emit a UserWarning.

decimate(X)#: Anti-aliasing low-pass (zero-phase) and downsampling by 2 along the last axis.

design_filters()#: Design the anti-aliasing (*_dec) and cutoff (*_filt) filters and the edge-extension length n_pad (samples at fs).

filter_signal(X)#

Low-frequency content of X (below cutoff), same shape as X.

See the class notes for the edge and NaN handling: every run of valid samples between NaN gaps is filtered on its own, and NaN samples stay NaN. Raises ValueError if X contains inf.

upsample(X)#: Upsample by 2 along the last axis: zero-insertion, then the anti-aliasing low-pass (zero-phase) with gain 2.

brainmaze_utils.signal.PSD(x: ndarray, fs: float, nperseg=None, noverlap=0, nfft=None)#

Estimates PSD of an input signal or signals using Welch’s method. If nperseg is None, the spectrum is estimated from the whole signal in a single window.

Parameters:

x (np.ndarray) – A single signal with a shape (n_samples) or set of signals with a shape (n_signals, n_shapes)
fs (float) – Sampling frequency
nperseg (int) – Number of samples for a segment
noverlap (int) – Number of overlap samples.

Returns:

freq (np.ndarray) – Frequency axis for estimated PSD
psd (np.ndarray) – Power spectral density estimate

brainmaze_utils.signal.buffer(x: ndarray, fs: float = 1, segm_size: float = None, overlap: float = 0, drop: bool = True)#

Cut a 1-D signal into (optionally overlapping) fixed-length segments.

Segment k starts at sample round(k * (segm_size - overlap) * fs) - start positions are computed from exact times, so there is no cumulative drift for non-integer fs * (segm_size - overlap). Every segment has round(segm_size * fs) samples.

Parameters:

x (numpy.ndarray) – 1-D signal (n_samples,). For multichannel data call per channel.
fs (float) – Sampling frequency in Hz.
segm_size (float, optional) – Segment length in seconds. If None, x is returned unchanged.
overlap (float) – Overlap between consecutive segments in seconds, 0 <= overlap < segm_size.
drop (bool) – If True (default) a trailing incomplete segment is dropped; otherwise one final segment covering the remaining samples is zero-padded (False for boolean input) to full length.

Returns:

(n_segments, round(segm_size * fs)) with the dtype of x. NaNs are copied as-is. If no segment fits, the result has shape (0, n_segm).

Return type:

numpy.ndarray

Raises:

ValueError – If x is not 1-D, or overlap is negative or >= segm_size (which previously caused an infinite loop), or the hop segm_size - overlap is shorter than one sample.

brainmaze_utils.signal.decimate(x, fs, fs_new, cutoff=None, datarate=False)#

Resample signal(s) to fs_new with a zero-phase anti-aliasing low-pass filter. NaN-aware. Despite the name it also upsamples (fs_new > fs).

Processing steps, per signal:

NaN samples are temporarily filled by linear interpolation between the neighbouring valid samples (leading/trailing NaNs take the nearest valid value) so the filter does not see steps at gap edges.
Zero-phase low-pass at cutoff: 16th-order Butterworth in second-order sections, applied forward and backward (scipy.signal.sosfiltfilt()); the magnitude response is therefore that of a 32nd-order filter with -6 dB at cutoff. When upsampling it is skipped if cutoff >= 0.45 * fs: the input holds nothing above fs / 2, and the interpolator (step 3) already band-limits to 0.45 * fs.
Output sample k is the filtered signal at time exactly k / fs_new s; the output has n_new = round(n_samples * fs_new / fs) samples.
- Integer ratios fs / fs_new pick every q-th sample (exact).
- Any other ratio, including upsampling and rates with decimals such as 30000.5 or 511.9999 Hz, evaluates the band-limited signal at the exact times with a Kaiser-windowed sinc interpolator (2 x 5 to 2 x 40 taps, designed for 120 dB). The position of every output sample is computed directly from k * fs / fs_new, so there is no timing drift, however long the record (timing error < 1e-6 samples). Interpolation error for content inside the band is ~1e-7 of its amplitude (measured on tones; the integer path is exact).
The NaN mask is re-applied: an output sample is NaN if any input sample within +-0.5 / fs_new s of it, or one of the two input samples bracketing it, was NaN. Signals that are entirely NaN stay entirely NaN.

Parameters:

x (numpy.ndarray) – (n_samples,) or (n_signals, n_samples) - samples along the last axis. Not modified.
fs (float) – Sampling frequency of x in Hz (> 0).
fs_new (float) – Target sampling frequency in Hz (> 0). May be lower, equal or higher than fs.
cutoff (float, optional) – Anti-aliasing cutoff in Hz. Default fs_new / 3 (two thirds of the new Nyquist frequency). When downsampling it must satisfy 0 < cutoff < fs_new / 2; a higher cutoff would let content above the new Nyquist frequency alias into the output. When upsampling it must be > 0; if it is >= 0.45 * fs no filter is applied.
datarate (bool) – If True, also return get_datarate() of the input.

Returns:

Resampled signal(s) with the same number of dimensions as x (float64), or (resampled, datarate) if datarate=True.

Return type:

numpy.ndarray or tuple

Raises:

ValueError – If fs or fs_new is not positive, if cutoff is outside the range given above, or if x has more than 2 dimensions.

Notes

fs_new == fs still applies the low-pass at cutoff (as in v2.0.0).
Upsampling: with the default cutoff the low-pass is applied only for ratios fs_new / fs < 1.35 (where fs_new / 3 < 0.45 * fs); v2.0.0 applied it up to 1.5, returned NaN between ~1.45 and 1.5 and raised above. The interpolator preserves content up to 0.45 * fs (or 1.25 * cutoff, if lower); content between 0.45 * fs and fs / 2 is attenuated. When the Butterworth is skipped (cutoff >= 0.45 * fs) such content is not removed but leaves an image at fs - f: a tone at 0.48 * fs comes out with gain ~0.93 plus a ~7 % image.
When upsampling, round(n * fs_new / fs) output samples can include up to 3 samples after the last input sample (e.g. 250->1000 Hz). They are extrapolated (odd extension) and are the least reliable part of the output: on a unit-amplitude tone their error is ~0.1 at 20 Hz and reaches ~1 for tones near 0.4 * fs. Discard them if the end of the record matters.
Samples close to (but outside) a NaN gap are computed from the interpolated fill and the filter’s impulse response, so a few output samples next to a gap may carry a small transient; they are not masked.
The record edges are extended by odd reflection over ~6 * fs / cutoff samples before filtering, which keeps edge transients small (~1e-3 of the amplitude for an in-band tone) but not zero.

Note

Changed after v2.0.0: Filter is now applied in second-order sections; the previous (b, a) 16th-order design was numerically unstable for fs / fs_new >= ~10 and returned all-NaN output (e.g. 3000->250, 1000->50, 32000->1000 Hz). NaNs are now re-applied to the output instead of being silently filled. The final resampling step no longer uses the FFT (scipy.signal.resample(), which assumes a periodic signal and wraps the end of the record into its start); it picks samples (integer ratios) or interpolates at exact times. Upsampling by ratios >= 1.5 used to raise in the filter design and now works. A cutoff at or above the new Nyquist frequency when downsampling now raises ValueError (it used to alias silently).

brainmaze_utils.signal.detrend(y, x=None, y2=None, method='lstsq')#

Remove a linear trend from a 1-D signal.

Parameters:

y (numpy.ndarray) – 1-D signal to detrend.
x (numpy.ndarray, optional) – Abscissa of y (same length). Default numpy.linspace(0, 1, len(y)).
y2 (numpy.ndarray, optional) – Second signal (same length) from which the trend fitted on y is subtracted as well.
method ({'lstsq', 'endpoints'}) –
- 'lstsq' (default): ordinary least-squares line fitted to all finite samples of y.
- 'endpoints': the line through the first and last sample (the behaviour of this function before the fix); sensitive to noise on the two endpoints.

Returns:

y - trend, or (y - trend, y2 - trend) if y2 is given. NaNs in the inputs stay NaN; they are ignored when fitting.

Return type:

numpy.ndarray or tuple of numpy.ndarray

Notes

Note

Changed after v2.0.0: Default changed from the endpoint line to a least-squares fit. With the endpoint line a single noisy first/last sample tilted the whole trend (e.g. a +10 outlier at y[0] left a residual slope of ~9.7 on a slope-5 ramp). Use method='endpoints' to reproduce the old results.

brainmaze_utils.signal.downsample_min_max(signal: ndarray, original_fs: float, final_fs: float) → tuple[ndarray, float]#

Downsamples an iEEG signal using the min-max method, preserving the temporal order of min and max values within each downsampling window.

The method processes the input signal in non-overlapping windows. For each window, it finds the minimum and maximum values and their original temporal order. These two values (min and max) are then placed in the output signal in their temporal order of appearance within the window.

Parameters:

signal (np.ndarray) – The input iEEG signal. Can be 1D (samples,) or 2D (channels, samples).
original_fs (float) – The original sampling rate of the signal in Hz.
final_fs (float) – The desired final sampling rate of the output points in Hz. Since each original window produces two points (a min and a max), the number of original signal windows processed per second is final_fs / 2.

Returns:

np.ndarray: The downsampled iEEG signal. If the input was 1D,
the output is 1D. If 2D, output is 2D.
float: The actual final sampling rate of the output signal in Hz.
This will be close to the requested final_fs but may differ slightly due to integer window sizes.

Return type:

tuple[np.ndarray, float]

Raises:

TypeError – If the input signal is not a NumPy array.
ValueError – If signal dimensions are incorrect, sampling rates are not positive, or if final_fs implies a window_size less than 1.

Note

The signal is processed in full windows. Any remaining samples at the end of the signal that do not form a complete window are ignored.

NaN samples are ignored when searching for the min/max of a window; a window that is entirely NaN produces two NaN output points. (Previously a single NaN made both points of its window NaN.)

brainmaze_utils.signal.fft_filter(X: ndarray, fs: float, cutoff: float, type: str = 'lp', edges=None, max_gap_fill=None)#

Ideal (brick-wall) FFT filter. NaN-aware.

The signal is transformed with an FFT along the last axis, every frequency bin on the rejected side of cutoff is set to zero, and the result is transformed back. Bin frequencies are the exact DFT frequencies k * fs / n_samples (numpy.fft.fftfreq()); positive and negative frequencies are treated symmetrically, so the output is real.

'lp' keeps bins with |f| <= cutoff (a bin exactly at cutoff is kept).
'hp' keeps bins with |f| > cutoff (DC is always removed by 'hp').

Parameters:

X (numpy.ndarray) – Signal (n_samples,) or a stack of signals (..., n_samples); filtered along the last axis. NaN marks gaps (see Notes). Not modified.
fs (float) – Sampling frequency in Hz.
cutoff (float) – Cutoff frequency in Hz (>= 0). A cutoff at or above Nyquist makes 'lp' the identity and 'hp' return zeros.
type (str) – 'lp' or 'hp'.
edges ({None, 'periodic', 'extend'}) –
What is assumed beyond the ends of the record (and of each segment between gaps):
- 'periodic': nothing is added; the FFT treats the record as one period, so a mismatch between its start and end rings (Gibbs) into both edges. This is the v2.0.0 behaviour.
- 'extend': the LowFrequencyFilter edge handling - a least-squares line is removed and added back to the low-pass part, and each end is extended by 5 / cutoff s (a local line plus the mirrored residual), so offset, drift and gap edges produce no jump. Falls back to 'periodic' for cutoff == 0 (DC removal has no edge transient).
- None (default): 'periodic' if X contains no NaN (unchanged from v2.0.0), 'extend' if it does.
Warning

With None the edge handling depends on whether X contains NaN: a single NaN anywhere switches the whole output from 'periodic' to 'extend', which changes samples far from the gap (e.g. up to ~1 for SD-172 1/f data, 1 Hz high-pass, 30-60 s from the NaN). Pass edges explicitly ('extend' is the better choice for EEG-like data) whenever the result must not depend on the presence of gaps.
max_gap_fill (float, optional) – Gaps shorter than this many seconds are bridged by linear interpolation; longer gaps split the signal into separately filtered segments. Default 0.5 / cutoff (0 for cutoff == 0). A level jump across a bridged gap is filtered like a real step; use a smaller value to split instead.

Returns:

Same shape as X (float64). NaN exactly where X is NaN.

Return type:

numpy.ndarray

Raises:

ValueError – If X contains inf, or for an invalid type, fs, cutoff or edges.

Notes

A brick-wall filter rings (Gibbs phenomenon) around sharp transients. With edges='periodic' no windowing or padding is applied.

Gaps (NaN). As in LowFrequencyFilter: interior gaps shorter than max_gap_fill are bridged linearly, filtered through and re-masked; at longer gaps the signal is split and each segment is filtered as if it were a complete record, with the chosen edges handling. Leading/trailing NaN are excluded. A level jump across a bridged gap is filtered as a real step, and a segment of only a few samples between long gaps gives finite but meaningless values (see LowFrequencyFilter). (An FFT of the whole array would spread a single NaN over the entire output.)

Note

Changed after v2.0.0: NaN input used to return an all-NaN output; gaps are now handled as above. edges and max_gap_fill are new; the default result for NaN-free input is that of v2.0.0 apart from the frequency-axis fix (see the Changes page).

brainmaze_utils.signal.find_peaks(y)#

Finds the peaks in a given signal.

Parameters:

y (numpy.ndarray) – The input signal in which to find peaks.

Returns:

position (numpy.ndarray) – The positions of the peaks in the input signal.
value (numpy.ndarray) – The values of the peaks in the input signal.

brainmaze_utils.signal.get_datarate(x)#

Fraction of valid (non-NaN) samples per signal.

Parameters:

x (numpy.ndarray or list of numpy.ndarray) – A single signal (n_samples,), a stack of signals (n_signals, n_samples) (samples along the last axis), or a list of 1-D signals of possibly different lengths.

Returns:

Values in [0, 1]: 1 - n_nan / n_samples.

1-D input -> float
N-D input -> numpy.ndarray of shape x.shape[:-1]
list input -> list of floats (one per element)

Return type:

float, numpy.ndarray or list

brainmaze_utils.signal.nandecimate(x, fs, fs_new, cutoff=None, datarate=False)#

Downsample signal(s) with a short FIR anti-aliasing filter. NaN-aware.

Legacy, lightweight variant of decimate(): NaNs are filled with the channel mean, the signal is low-passed with a 30-tap FIR (firwin, applied forward-backward), resampled with scipy.signal.resample() (FFT), and the NaN mask is re-applied (an output sample is NaN if any input sample within +-0.5 / fs_new s of it was NaN).

Warning

A 30-tap FIR has a wide transition band, so anti-aliasing is weak for large ratios fs / fs_new. Prefer decimate() for analysis.

Parameters:

x (numpy.ndarray) – (n_samples,) or (n_signals, n_samples) - samples along the last axis. Not modified.
fs (float) – Sampling frequency of x in Hz.
fs_new (float) – Target sampling frequency in Hz.
cutoff (float, optional) – FIR cutoff in Hz. Default fs_new / 3.
datarate (bool) – If True, also return get_datarate() of the input.

Returns:

Decimated signal(s) (round(n_samples * fs_new / fs) samples), or (decimated, datarate) if datarate=True.

Return type:

numpy.ndarray or tuple

Notes

Note

Changed after v2.0.0: NaNs were previously filled with the NaN fraction of the channel (a value in [0, 1]) instead of the channel mean, producing large spurious transients next to every gap for signals with a DC offset.

brainmaze_utils.signal.resample(x, fsamp_orig, fsamp_new)#

Resample a signal to a new sampling frequency by linear interpolation. NaN-aware.

Warning

No anti-aliasing filter is applied - by design. resample only interpolates. When downsampling, any content above the new Nyquist frequency fsamp_new / 2 folds back (aliases) into the output: e.g. a 400 Hz tone resampled 1000 -> 250 Hz appears as a full-amplitude 100 Hz tone. Band-limiting the signal is deliberately left to the caller, who knows which band matters. Before downsampling, low-pass below fsamp_new / 2 (with margin), or use decimate(), which filters and downsamples in one step.

import scipy.signal as ss
from brainmaze_utils.signal import resample, decimate

fs, fs_new = 1000, 250
# option 1: your own zero-phase low-pass, then resample
sos = ss.butter(8, 0.8 * fs_new / 2, 'lp', fs=fs, output='sos')
y = resample(ss.sosfiltfilt(sos, x), fs, fs_new)  # x must be NaN-free here
# option 2: decimate (anti-aliasing filter included, NaN-aware)
y = decimate(x, fs, fs_new)

Upsampling (fsamp_new > fsamp_orig) does not need a filter; linear interpolation slightly attenuates content close to the original Nyquist.

Sample i of the input is at time i / fsamp_orig and sample k of the output at k / fsamp_new (both starting at 0); the output has round(n_samples * fsamp_new / fsamp_orig) samples. Each output value is the linear interpolation between the two input samples that bracket its time.

Parameters:

x (numpy.ndarray) – (n_samples,) or (..., n_samples) - resampled along the last axis. Not modified.
fsamp_orig (float) – Sampling frequency of x in Hz.
fsamp_new (float) – Target sampling frequency in Hz.

Returns:

Resampled signal (float64). If fsamp_orig == fsamp_new a float copy of x is returned.

Return type:

numpy.ndarray

Notes

NaN handling: an output sample is NaN if either bracketing input sample is NaN (an output time that coincides with an input sample uses only that sample). Values are never interpolated across a gap. When downsampling, short gaps can disappear: a gap that falls entirely between two output instants marks no output sample (e.g. a 3-sample gap, 1000 -> 250 Hz, gives no NaN at all). Use decimate() if the gap mask must be preserved (it flags every output sample within half an output period of a NaN).
Empty input (or an output of 0 samples) returns an empty array of shape x.shape[:-1] + (n_new,).
Output times beyond the last input sample (possible when upsampling) take the last input value.

Note

Changed after v2.0.0: Time axes were linspace(0, 1, N) / linspace(0, 1, N_new), which pins both endpoints and stretched time: the effective output rate was (N_new - 1) / (N - 1) * fs rather than fsamp_new. Also fixed np.NaN (removed in NumPy 2) and division by zero for constant signals.

brainmaze_utils.signal.unify_sampling_frequency(x: list, sampling_frequency: list, fs_new=None) → tuple#

Bring a list of signals to a common sampling frequency using decimate().

If all frequencies are equal and fs_new is None: nothing is done.
If frequencies differ and fs_new is None: all signals are decimated to the lowest frequency present.
If fs_new is given: every signal whose frequency differs from fs_new is resampled to fs_new with decimate(), which downsamples or upsamples (e.g. a 150 Hz channel brought to fs_new=200).

Signals already at the target frequency are passed through unchanged (they are not low-passed).

Parameters:

x (list of numpy.ndarray) – Signals (1-D, or 2-D with samples along the last axis).
sampling_frequency (list or numpy.ndarray) – Sampling frequency in Hz of each signal.
fs_new (float, optional) – Target sampling frequency in Hz. Signals below it are upsampled, signals above it are low-passed (fs_new / 3) and downsampled.

Returns:

(list of numpy.ndarray, fs_new). A new list is returned; the input list and its arrays are not modified.

Return type:

tuple

Signal

Contents

Signal#