cedalion.sigproc.quality

Signal quality metrics and channel pruning functionality.

Functions

`detect_baselineshift`(ts, outlier_mask)	Detect baselineshifts in fNIRSdata.
`detect_outliers`(ts, t_window_std[, ...])	Detect outliers in fNIRSdata based on standard deviation and gradient of signal.
`detect_outliers_grad`(ts[, iqr_threshold])	Detect outliers in fNIRSdata based on gradient of signal.
`detect_outliers_std`(ts, t_window[, ...])	Detect outliers in fNIRSdata based on standard deviation of signal.
`gvtd`(amplitudes[, stat_type, n_std])	Calculate GVTD metric based on Sherafati et al. [SSE+20].
`id_motion`(fNIRSdata[, t_motion, t_mask, ...])	Identify motion artifacts in an fNIRS input data array.
`id_motion_refine`(ma_mask, operator)	Refines motion artifact mask to simplify and quantify motion artifacts.
`mean_amp`(amplitudes, amp_range)	Calculate mean amplitudes and mask channels outside amplitude range.
`prune_ch`(amplitudes, masks, operator[, ...])	Prune channels from the the input data array using quality masks.
`psp`(amplitudes, window_length, psp_thresh[, ...])	Calculate the peak spectral power.
`repair_amp`(amp[, median_len, interp_nan])	Replace nonpositive amp values and optionally fill NaNs.
`sci`(amplitudes, window_length, sci_thresh[, ...])	Calculate the scalp-coupling index.
`sd_dist`(amplitudes, geo3D[, sd_range])	Calculate source-detector separations and mask channels outside a distance range.
`snr`(amplitudes[, snr_thresh])	Calculates signal-to-noise ratio for each channel and other dimension.
`stimulus_mask`(df_stim, mask)	Create a mask which events overlap with periods flagged as tainted in mask.

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.