cedalion.sigproc package

Submodules

cedalion.sigproc.epochs module

Extract epochs from a time series based on stimulus events.

cedalion.sigproc.epochs.to_epochs(
ts: cdt.NDTimeSeries,
df_stim: pd.DataFrame,
trial_types: list[str],
before: cdt.QTime,
after: cdt.QTime,
)[source]

Extract epochs from the time series based on stimulus events.

Parameters:

  • ts – the time series

  • df_stim – DataFrame containing stimulus events.

  • trial_types – List of trial types to include in the epochs.

  • before – Time before stimulus event to include in epoch.

  • after – Time after stimulus event to include in epoch.

Returns:

Array containing the extracted epochs.

Return type:

xarray.DataArray

cedalion.sigproc.frequency module

Frequency-related signal processing methods.

cedalion.sigproc.frequency.sampling_rate(
timeseries: Annotated[DataArray, DataArraySchema(dims='time', coords='time', 'time', 'samples')],
) → Quantity[source]

Estimate the sampling rate of the timeseries.

Note

This functions assumes uniform sampling.

Parameters:

timeseries (NDTimeSeries, (time,*)) – the input time series

Returns:

The sampling rate estimated by averaging time differences between samples.

cedalion.sigproc.frequency.freq_filter(
timeseries: Annotated[DataArray, DataArraySchema(dims='time', coords='time', 'time', 'samples')],
fmin: Annotated[Quantity, '[frequency]'],
fmax: Annotated[Quantity, '[frequency]'],
butter_order: int = 4,
) → Annotated[DataArray, DataArraySchema(dims='time', coords='time', 'time', 'samples')][source]

Apply a Butterworth bandpass frequency filter.

Parameters:

  • timeseries (NDTimeSeries, (time,*)) – the input time series

  • fmin (Quantity, [frequency]) – lower threshold of the pass band

  • fmax (Quantity, [frequency]) – higher threshold of the pass band

  • butter_order – order of the filter

Returns:

The frequency-filtered time series

cedalion.sigproc.motion_correct module

Various algorithms for motion correction of fNIRS data.

cedalion.sigproc.motion_correct.motion_correct_spline(
fNIRSdata: cdt.NDTimeSeries,
tIncCh: cdt.NDTimeSeries,
p: float,
) → cdt.NDTimeSeries[source]

Apply motion correction using spline interpolation to fNIRS data.

Based on Homer3 [1] v1.80.2 “hmrR_tInc_baselineshift_Ch_Nirs.m” Boston University Neurophotonics Center https://github.com/BUNPC/Homer3

Parameters:

  • fNIRSdata – The fNIRS data to be motion corrected.

  • tIncCh – The time series indicating the presence of motion artifacts.

  • p – smoothing factor

Returns:

The motion-corrected fNIRS data.

Return type:

dodSpline (cdt.NDTimeSeries)

cedalion.sigproc.motion_correct.compute_window(
SegLength: cdt.NDTimeSeries,
dtShort: Quantity,
dtLong: Quantity,
fs: Quantity,
)[source]

Computes the window size.

Window size is based on the segment length, short time interval, long time interval, and sampling frequency.

Parameters:

  • SegLength (cdt.NDTimeSeries) – The length of the segment.

  • dtShort (Quantity) – The short time interval.

  • dtLong (Quantity) – The long time interval.

  • fs (Quantity) – The sampling frequency.

Returns:

The computed window size.

Return type:

wind

cedalion.sigproc.motion_correct.motion_correct_splineSG(
fNIRSdata: cdt.NDTimeSeries,
p: float,
frame_size: Quantity = <Quantity(10,
'second')>,
)[source]

Apply motion correction using spline interpolation and Savitzky-Golay filter.

Parameters:

  • fNIRSdata (cdt.NDTimeSeries) – The fNIRS data to be motion corrected.

  • frame_size (Quantity) – The size of the sliding window in seconds for the Savitzky-Golay filter. Default is 10 seconds.

  • p – smoothing factor

Returns:

The motion-corrected fNIRS data after applying spline interpolation and Savitzky-Golay filter.

Return type:

dodSplineSG (cdt.NDTimeSeries)

cedalion.sigproc.motion_correct.motion_correct_PCA(
fNIRSdata: cdt.NDTimeSeries,
tInc: cdt.NDTimeSeries,
nSV: Quantity = 0.97,
)[source]

Apply motion correction using PCA filter identified as motion artefact segments.

Based on Homer3 [1] v1.80.2 “hmrR_MotionCorrectPCA.m” Boston University Neurophotonics Center https://github.com/BUNPC/Homer3

Inputs:

fNIRSdata: The fNIRS data to be motion corrected. tInc: The time series indicating the presence of motion artifacts. nSV (Quantity): Specifies the number of prinicpal components to remove from the

data. If nSV < 1 then the filter removes the first n components of the data that removes a fraction of the variance up to nSV.

Returns:

The motion-corrected fNIRS data. svs (np.array): the singular values of the PCA. nSV (Quantity): the number of principal components removed from the data.

Return type:

fNIRSdata_cleaned (cdt.NDTimeSeries)

cedalion.sigproc.motion_correct.motion_correct_PCA_recurse(
fNIRSdata: cdt.NDTimeSeries,
t_motion: Quantity = 0.5,
t_mask: Quantity = 1,
stdev_thresh: Quantity = 20,
amp_thresh: Quantity = 5,
nSV: Quantity = 0.97,
maxIter: Quantity = 5,
)[source]

Identify motion artefacts in input fNIRSdata.

If any active channel exhibits signal change greater than STDEVthresh or AMPthresh, then that segment of data is marked as a motion artefact. motion_correct_PCA is applied to all segments of data identified as a motion artefact. This is called until maxIter is reached or there are no motion artefacts identified.

Parameters:

  • fNIRSdata (cdt.NDTimeSeries) – The fNIRS data to be motion corrected.

  • t_motion – check for signal change indicative of a motion artefact over time range tMotion. (units of seconds)

  • t_mask (Quantity) – mark data +/- tMask seconds aroundthe identified motion artefact as a motion artefact.

  • stdev_thresh (Quantity) – if the signal d for any given active channel changes by more than stdev_thresh * stdev(d) over the time interval tMotion then this time point is marked as a motion artefact.

  • amp_thresh (Quantity) – if the signal d for any given active channel changes by more than amp_thresh over the time interval tMotion then this time point is marked as a motion artefact.

  • nSV – FIXME

  • maxIter – FIXME

Returns:

The motion-corrected fNIRS data. svs (np.array): the singular values of the PCA. nSV (int): the number of principal components removed from the data.

Return type:

fNIRSdata_cleaned (cdt.NDTimeSeries)

cedalion.sigproc.motion_correct.tddr(ts: cdt.NDTimeSeries)[source]

Implementation of the TDDR algorithm for motion correction.

Uses an iterative reweighting approach to reduce large fluctuations typically associated with motion artifacts. Adapted for cedalion from the python implementation at [Fis18], which is the reference implementation for the algorithm described in [FLVM19].

Parameters:

ts – The time series to be corrected. Should have dims channel and wavelength

Returns:

The corrected time series.

References

Paper: [FLVM19] Code: [Fis18]

cedalion.sigproc.motion_correct.pad_to_power_2(signal)[source]

Pad signal to next power of 2.

cedalion.sigproc.motion_correct.process_coefficients(coeffs, iqr_factor, signal_length)[source]

Deletes outlier coefficients based on IQR.

cedalion.sigproc.motion_correct.mad(x)[source]

Compute Median Absolute Deviation.

cedalion.sigproc.motion_correct.normalize_signal(signal, wavelet='db2')[source]

Normalize signal by its noise level using MAD of downsampled coefficients.

Implements Homer3’s NormalizationNoise function.

Parameters:

  • signal – 1D numpy array containing the signal to normalize

  • wavelet – wavelet to use (default: ‘db2’)

Returns:

normalized version of input signal norm_coef: normalization coefficient (multiply by this to denormalize)

Return type:

normalized_signal

References

[HDFB09]

cedalion.sigproc.motion_correct.motion_correct_wavelet(od, iqr=1.5, wavelet='db2', level=4)[source]

Wavelet-based motion correction, specializing in spike correction.

Implements the wavelet-based motion correction algorithm described in [MD12], closely following the MATLAB implementation found in Homer3 ([HDFB09])

Parameters:

  • od – The time series to be corrected. Should have dims channel and wavelength

  • iqr – The interquartile range factor for outlier detection. Set to -1 to disable. Increasing iqr will delete less coefficients.

  • wavelet – The wavelet to use for decomposition (default: ‘db2’)

  • level – The level of decomposition to use (default: 4)

Returns:

The corrected time series.

References

Original paper: [MD12] Implementation based on Homer3 v1.80.2 “hmrR_MotionCorrectWavelet.m” and its dependencies ([HDFB09]).

cedalion.sigproc.quality module

Signal quality metrics and channel pruning functionality.

cedalion.sigproc.quality.prune_ch(
amplitudes: cdt.NDTimeSeries,
masks: list[cdt.NDTimeSeries],
operator: str,
flag_drop: bool = True,
)[source]

Prune channels from the the input data array using quality masks.

Parameters:

  • amplitudes (NDTimeSeries) – input time series

  • masks (list[NDTimeSeries]) – list of boolean masks with coordinates comptabile to amplitudes

  • operator –

    operators for combination of masks before pruning data_array

    • ”all”: logical AND, keeps channel if it is good across all masks

    • ”any”: logical OR, keeps channel if it is good in any mask/metric

  • flag_drop – if True, channels are dropped from the data_array, otherwise they are set to NaN (default: True)

Returns:

A tuple (amplitudes_pruned, prune_list), where amplitudes_pruned is a the original time series channels pruned (dropped) according to quality masks. A list of removed channels is returned in prune_list.

cedalion.sigproc.quality.psp(
amplitudes: NDTimeSeries,
window_length: cdt.QTime,
psp_thresh: float,
cardiac_fmin: cdt.QFrequency = <Quantity(0.5,
'hertz')>,
cardiac_fmax: cdt.QFrequency = <Quantity(2.5,
'hertz')>,
)[source]

Calculate the peak spectral power.

The peak spectral power metric is based on Pollonini et al. [POA+14] / Pollonini et al. [PBO16].

Parameters:

  • amplitudes (NDTimeSeries, (channel, wavelength, time)) – input time series

  • window_length (Quantity, [time]) – size of the computation window

  • psp_thresh – if the calculated PSP metric falls below this threshold then the corresponding time window should be excluded.

  • cardiac_fmin – minimm frequency to extract cardiac component

  • cardiac_fmax – maximum frequency to extract cardiac component

Returns:

A tuple (psp, psp_mask), where psp is a DataArray with coords from the input NDTimeseries containing the peak spectral power. psp_mask is a boolean mask DataArray with coords from psp, true where psp_thresh is met.

cedalion.sigproc.quality.gvtd(amplitudes: cdt.NDTimeSeries, stat_type: str = 'default', n_std: int = 10)[source]

Calculate GVTD metric based on Sherafati et al. [SSE+20].

Parameters:

  • amplitudes (NDTimeSeries, (channel, wavelength, time)) – input time series

  • stat_type (string) – statistic of GVTD time trace to use to set the threshold (see _get_gvtd_threshold). Default = ‘default’

  • n_std (int) – number of standard deviations for consider above the statistic of interest.

Returns:

A DataArray with coords from the input NDTimeseries containing the GVTD metric.

cedalion.sigproc.quality.sci(
amplitudes: NDTimeSeries,
window_length: Quantity,
sci_thresh: float,
cardiac_fmin: cdt.QFrequency = <Quantity(0.5,
'hertz')>,
cardiac_fmax: cdt.QFrequency = <Quantity(2.5,
'hertz')>,
)[source]

Calculate the scalp-coupling index.

The scalp-coupling index metric is based on Pollonini et al. [POA+14] / Pollonini et al. [PBO16].

Parameters:

  • amplitudes (NDTimeSeries, (channel, wavelength, time)) – input time series

  • window_length (Quantity, [time]) – size of the computation window

  • sci_thresh – if the calculated SCI metric falls below this threshold then the corresponding time window should be excluded.

  • cardiac_fmin – minimm frequency to extract cardiac component

  • cardiac_fmax – maximum frequency to extract cardiac component

Returns:

A tuple (sci, sci_mask), where sci is a DataArray with coords from the input NDTimeseries containing the scalp-coupling index. Sci_mask is a boolean mask DataArray with coords from sci, true where sci_thresh is met.

cedalion.sigproc.quality.snr(amplitudes: cdt.NDTimeSeries, snr_thresh: float = 2.0)[source]

Calculates signal-to-noise ratio for each channel and other dimension.

SNR is the ratio of the average signal over time divided by its standard deviation.

Parameters:

  • amplitudes (NDTimeSeries, (time, *)) – the input time series

  • snr_thresh – threshold (unitless) below which a channel should be excluded.

Returns:

A tuple (snr, snr_mask) , where snr is a DataArray with coords from input NDTimeseries containing the ratio between mean and std of the amplitude signal for all channels. snr_mask is a boolean mask DataArray with the same coords as snr, true where snr_threshold is met.

References

Based on Homer3 v1.80.2 “hmR_PruneChannels.m” (Huppert et al. [HDFB09])

cedalion.sigproc.quality.mean_amp(amplitudes: cdt.NDTimeSeries, amp_range: tuple[Quantity, Quantity])[source]

Calculate mean amplitudes and mask channels outside amplitude range.

Parameters:

  • amplitudes (NDTimeSeries, (time, *)) – input time series

  • amp_range – if amplitudes.mean(“time”) < amp_threshs[0] or > amp_threshs[1] then it is excluded as an active channel in amp_mask

Returns:

A tuple (mean_amp, amp_mask), where mean_amp is DataArray with coords from input NDTimeseries containing the amplitudes averaged over time. The boolean DataArray amp_mask is true where amp_threshs are met

References

Based on Homer3 v1.80.2 “hmR_PruneChannels.m” (Huppert et al. [HDFB09])

cedalion.sigproc.quality.sd_dist(
amplitudes: cdt.NDTimeSeries,
geo3D: cdt.LabeledPointCloud,
sd_range: tuple[Quantity,
Quantity] = (<Quantity(0,
'centimeter')>,
<Quantity(4.5,
'centimeter')>),
)[source]

Calculate source-detector separations and mask channels outside a distance range.

Parameters:

  • amplitudes (NDTimeSeries, (channel, *)) – input time series

  • geo3D (LabeledPointCloud) – 3D optode coordinates

  • sd_range – if source-detector separation < sd_range[0] or > sd_range[1] then it is excluded as an active channelin sd_mask

Returns:

A tuple (sd_dist, sd_mask), where sd_dist contains the channel distances and sd_mask is a boolean NDTimeSeries, indicating where distances fall into sd_range.

References

Based on Homer3 v1.80.2 “hmR_PruneChannels.m” (Huppert et al. [HDFB09])

cedalion.sigproc.quality.id_motion(
fNIRSdata: cdt.NDTimeSeries,
t_motion: Quantity = <Quantity(0.5,
'second')>,
t_mask: Quantity = <Quantity(1.0,
'second')>,
stdev_thresh: float = 50.0,
amp_thresh: float = 5.0,
) → cdt.NDTimeSeries[source]

Identify motion artifacts in an fNIRS input data array.

If any active data channel exhibits a signal change greater than std_thresh or amp_thresh, then a segment of data around that time point is marked as a motion artifact.

Parameters:

  • fNIRSdata (NDTimeSeries, (time, channel, *)) – input time series

  • t_motion (Quantity, [time]) – time interval for motion artifact detection. Checks for signal change indicative of a motion artifact over time range t_motion.

  • t_mask (Quantity, [time]) – time range to mask around motion artifacts. Mark data over +/- t_mask around the identified motion artifact as a motion artifact.

  • stdev_thresh – threshold for std deviation of signal change. If the signal d for any given active channel changes by more than stdev_thresh * stdev(d) over the time interval tMotion, then this time point is marked as a motion artifact. The standard deviation is determined for each channel independently. Typical value ranges from 5 to 20. Use a value of 100 or greater if you wish for this condition to not find motion artifacts.

  • amp_thresh – threshold for amplitude of signal change. If the signal d for any given active channel changes by more that amp_thresh over the time interval t_motion, then this time point is marked as a motion artifact. Typical value ranges from 0.01 to 0.3 for optical density units. Use a value greater than 100 if you wish for this condition to not find motion artifacts.

Returns:

a DataArray that has at least the dimensions channel and time, dtype is boolean. At each time point, CLEAN indicates data included and TAINTED indicates motion artifacts.

References

Based on Homer3 v1.80.2 “hmR_MotionArtifact.m” and “hmR_MotionArtifactByChannel” (Huppert et al. [HDFB09]).

cedalion.sigproc.quality.id_motion_refine(ma_mask: cdt.NDTimeSeries, operator: str)[source]

Refines motion artifact mask to simplify and quantify motion artifacts.

Parameters:

  • ma_mask (NDTimeSeries, (time, channel, *)) – motion artifact mask as generated by id_motion().

  • operator –

    operation to apply to the mask. Available operators:

    • by_channel: collapses the mask along the amplitude/wavelength/chromo dimension to provide a single motion artifact mask per channel (default) over time

    • all: collapses the mask along all dimensions to provide a single motion artifact marker for all channels over time i.e. an artifact detected in any channel masks all channels.

Returns:

A tuple (ma_mask_new, ma_info), where ma_mask_new is the updated motion artifact mask, collapsed according to operator and ma_info is a pandas.DataFrame that contains 1) channels with motion artifacts, 2) # of artifacts detected per channel and 3) fraction of artifacts/total time.

cedalion.sigproc.quality.detect_outliers_std(ts: cdt.NDTimeSeries, t_window: cdt.QTime, iqr_threshold=2)[source]

Detect outliers in fNIRSdata based on standard deviation of signal.

:param ts NDTimeSeries: fNIRS timeseries data :type ts NDTimeSeries: time, channel, * :param : fNIRS timeseries data :type : time, channel, * :param t_window Quantity: time window over which to calculate std. deviations :param iqr_threshold: interquartile range threshold (detect outlier as any std.

deviation outside iqr_threshold * [25th percentile, 75th percentile])

Returns:

mask that is a DataArray containing TRUE anywhere the data is clean and FALSE anytime an outlier is detected based on the standard deviation

References

Based on Homer3 v1.80.2 “hmrR_tInc_baselineshift_Ch_Nirs.m” (Jahani et al. [JSBYucel18])

cedalion.sigproc.quality.detect_outliers_grad(ts: cdt.NDTimeSeries, iqr_threshold: float = 1.5)[source]

Detect outliers in fNIRSdata based on gradient of signal.

Parameters:

  • ts (NDTimeSeries, (time, channel, *)) – fNIRS timeseries data

  • iqr_threshold – interquartile range threshold (detect outlier as any gradient outside iqr_threshold * [25th percentile, 75th percentile])

Returns:

mask that is a DataArray containing TRUE anywhere the data is clean and FALSE anytime an outlier is detected

References

Based on Homer3 v1.80.2 “hmrR_tInc_baselineshift_Ch_Nirs.m” (Jahani et al. [JSBYucel18])

cedalion.sigproc.quality.detect_outliers(
ts: cdt.NDTimeSeries,
t_window_std: cdt.QTime,
iqr_threshold_std: float = 2,
iqr_threshold_grad: float = 1.5,
)[source]

Detect outliers in fNIRSdata based on standard deviation and gradient of signal.

Parameters:

  • ts (NDTimeSeries, (time, channel, *)) – fNIRS timeseries data

  • t_window_std (Quantity) – time window over which to calculate std. devs.

  • iqr_threshold_grad – interquartile range threshold (detect outlier as any gradient outside iqr_threshold * [25th percentile, 75th percentile])

  • iqr_threshold_std – interquartile range threshold (detect outlier as any standard deviation outside iqr_threshold * [25th percentile, 75th percentile])

Returns:

mask that is a DataArray containing TRUE anywhere the data is clean and FALSE anytime an outlier is detected

References

Based on Homer3 v1.80.2 “hmrR_tInc_baselineshift_Ch_Nirs.m” (Jahani et al. [JSBYucel18])

cedalion.sigproc.quality.detect_baselineshift(ts: cdt.NDTimeSeries, outlier_mask: cdt.NDTimeSeries)[source]

Detect baselineshifts in fNIRSdata.

Parameters:

  • ts (NDTimeSeries, (time, channel, *)) – fNIRS timeseries data

  • outlier_mask (NDTimeSeries) – mask containing FALSE anytime an outlier is detected in signal

Returns:

mask that is a DataArray containing TRUE anywhere the data is clean and FALSE anytime a baselineshift or outlier is detected.

References

Based on Homer3 v1.80.2 “hmrR_tInc_baselineshift_Ch_Nirs.m” (Jahani et al. [JSBYucel18])

cedalion.sigproc.quality.stimulus_mask(
df_stim: DataFrame,
mask: DataArray,
) → DataArray[source]

Create a mask which events overlap with periods flagged as tainted in mask.

Parameters:

  • df_stim – stimulus data frame

  • mask – signal quality mask. Must contain dimensions ‘channel’ and ‘time’

Returns:

A boolean mask with dimensions “stim”, “channel”. The stim dimension matches the stimulus dataframe. Stimuli are marked as TAINTED when there is any TAINTED flag in the mask between onset and onset+ duration.

cedalion.sigproc.quality.repair_amp(
amp: DataArray,
median_len=3,
interp_nan=True,
**kwargs,
)[source]

Replace nonpositive amp values and optionally fill NaNs.

TODO: Optimize handling of sequential nonpositive values.

Parameters:

  • amp – Amplitude data

  • median_len – Window size for the median filter

  • interp_nan – If True, interpolate NaNs in the data

  • **kwargs – Additional arguments for xarray interpolate_na function, such as method = “linear” (default), method = “nearest”, etc. See xarray documentation for more details.

cedalion.sigproc.tasks module

cedalion.sigproc.tasks.int2od(
rec: Recording,
ts_input: str | None = None,
ts_output: str = 'od',
)[source]

Calculate optical density from intensity amplitude data.

Parameters:

  • rec (Recording) – container of timeseries data

  • ts_input (str) – name of intensity timeseries. If None, this tasks operates on the last timeseries in rec.timeseries.

  • ts_output (str) – name of optical density timeseries.

cedalion.sigproc.tasks.od2conc(
rec: Recording,
dpf: dict[float, float],
spectrum: str = 'prahl',
ts_input: str | None = None,
ts_output: str = 'conc',
)[source]

Calculate hemoglobin concentrations from optical density data.

Parameters:

  • rec (Recording) – container of timeseries data

  • dpf (dict[float, float]) – differential path length factors

  • spectrum (str) – label of the extinction coefficients to use (default: “prahl”)

  • ts_input (str | None) – name of intensity timeseries. If None, this tasks operates on the last timeseries in rec.timeseries.

  • ts_output (str) – name of optical density timeseries (default: “conc”).

cedalion.sigproc.tasks.snr(
rec: Recording,
snr_thresh: float = 2.0,
ts_input: str | None = None,
aux_obj_output: str = 'snr',
mask_output: str = 'snr',
)[source]

Calculate signal-to-noise ratio (SNR) of timeseries data.

Parameters:

  • rec (Recording) – The recording object containing the data.

  • snr_thresh (float) – The SNR threshold.

  • ts_input (str | None, optional) – The input time series. Defaults to None.

  • aux_obj_output (str, optional) – The key for storing the SNR in the auxiliary object. Defaults to “snr”.

  • mask_output (str, optional) – The key for storing the mask in the recording object. Defaults to “snr”.

cedalion.sigproc.tasks.sd_dist(
rec: Recording,
sd_min: Annotated[Quantity, '[length]'],
sd_max: Annotated[Quantity, '[length]'],
ts_input: str | None = None,
aux_obj_output: str = 'sd_dist',
mask_output: str = 'sd_dist',
)[source]

Calculate source-detector separations and mask channels outside a range.

Parameters:

  • rec (Recording) – The recording object containing the data.

  • sd_min (Annotated[Quantity, "[length]"]) – The minimum source-detector separation.

  • sd_max (Annotated[Quantity, "[length]"]) – The maximum source-detector separation.

  • ts_input (str | None, optional) – The input time series. Defaults to None.

  • aux_obj_output (str, optional) – The key for storing the source-detector distances in the auxiliary object. Defaults to “sd_dist”.

  • mask_output (str, optional) – The key for storing the mask in the recording object. Defaults to “sd_dist”.

Module contents