Channel Quality Assessment and Pruning

This notebook sketches how to prune bad channels

[1]:

import matplotlib.pyplot as p
import numpy as np
from matplotlib.colors import LinearSegmentedColormap

import cedalion
import cedalion.datasets as datasets
import cedalion.nirs
import cedalion.plots as plots
import cedalion.sigproc.quality as quality
import cedalion.xrutils as xrutils
from cedalion import units

Loading raw CW-NIRS data from a SNIRF file and converting it to OD and CONC

This notebook uses a finger-tapping dataset in BIDS layout provided by Rob Luke that is automatically fetched. You can also find it here.

[2]:

# get example finger tapping dataset

rec = datasets.get_fingertapping()
rec["od"] = cedalion.nirs.int2od(rec["amp"])

# Plot some data for visual validation
f,ax = p.subplots(1,1, figsize=(12,4))
ax.plot( rec["amp"].time, rec["amp"].sel(channel="S1D1", wavelength="850"), "r-", label="850nm")
ax.plot( rec["amp"].time, rec["amp"].sel(channel="S1D1", wavelength="760"), "r-", label="760nm")
p.legend()
ax.set_xlabel("time / s")
ax.set_ylabel("Signal intensity / a.u.")

display(rec["amp"])

<xarray.DataArray (channel: 28, wavelength: 2, time: 23239)> Size: 10MB
<Quantity([[[0.0913686 0.0909875 0.0910225 ... 0.0941083 0.0940129 0.0944882]
  [0.1856806 0.186377  0.1836514 ... 0.1856486 0.1850836 0.1842172]]

 [[0.227516  0.2297024 0.2261366 ... 0.2264519 0.2271665 0.226713 ]
  [0.6354927 0.637668  0.6298023 ... 0.6072068 0.6087293 0.6091066]]

 [[0.1064704 0.1066212 0.1053444 ... 0.121114  0.1205022 0.1205441]
  [0.2755033 0.2761615 0.2727006 ... 0.2911952 0.2900544 0.2909847]]

 ...

 [[0.2027881 0.1996586 0.2004866 ... 0.2318743 0.2311941 0.2330808]
  [0.4666358 0.4554404 0.4561614 ... 0.4809749 0.4812827 0.4862896]]

 [[0.4885007 0.4802285 0.4818338 ... 0.6109142 0.6108118 0.613845 ]
  [0.8457658 0.825988  0.8259648 ... 0.975894  0.9756599 0.9826459]]

 [[0.6304559 0.6284427 0.6287045 ... 0.6810626 0.6809573 0.6818709]
  [1.2285622 1.2205907 1.2190002 ... 1.2729124 1.2727222 1.2755645]]], 'volt')>
Coordinates:
  * time        (time) float64 186kB 0.0 0.128 0.256 ... 2.974e+03 2.974e+03
    samples     (time) int64 186kB 0 1 2 3 4 5 ... 23234 23235 23236 23237 23238
  * channel     (channel) object 224B 'S1D1' 'S1D2' 'S1D3' ... 'S8D8' 'S8D16'
    source      (channel) object 224B 'S1' 'S1' 'S1' 'S1' ... 'S8' 'S8' 'S8'
    detector    (channel) object 224B 'D1' 'D2' 'D3' 'D9' ... 'D7' 'D8' 'D16'
  * wavelength  (wavelength) float64 16B 760.0 850.0
Attributes:
    data_type_group:  unprocessed raw

../../_images/examples_signal_quality_21_data_quality_and_pruning_3_1.png

Calculating Signal Quality Metrics and applying Masks

To assess channel quality metrics such as SNR, channel distances, average amplitudes, sci, and others, we use small helper functions. As input, the quality functions should also expect thresholds for these metrics, so they can feed back both the calculated quality metrics as well as a mask. The masks can then be combined and applied - e.g. to prune channels with low SNR. The input and output arguments are based on xarray time series, quality parameters / instructions for thresholding. The returned mask is a boolean array in the shape and size of the input time series. It indicates where the threshold for our quality metric was passed (“True”) and is False otherwise. Masks can be combined with other masks, for instance to apply several metrics to assess or prune channels. At any point in time, the mask can be applied using the “apply_mask()” function available from cedalion’s the xrutils package.

If you are a user who is mainly interested in high-level application, you can skip to the Section “Channel Pruning using Quality Metrics and the Pruning Function” below. The “prune_ch()” function provides a higher abstraction layer to simply prune your data, using the same metrics and functions that are demonstrated below.

Channel Quality Metrics: SNR

[3]:

# Here we assess channel quality by SNR
snr_thresh = 16 # the SNR (std/mean) of a channel. Set high here for demonstration purposes

# SNR thresholding using the "snr" function of the quality subpackage
snr, snr_mask = quality.snr(rec["amp"], snr_thresh)

# apply mask function. In this example, we want all signals with an SNR below the threshold to be replaced with "nan".
# We do not want to collapse / combine any dimension of the mask (last argument: "none")
data_masked_snr_1, masked_elements_1 = xrutils.apply_mask(rec["amp"], snr_mask, "nan", "none")

# alternatively, we can "drop" all channels with an SNR below the threshold. Since the SNR of both wavelength might differ
# (pass the threshold for one wavelength, but not for the other), we collapse to the "channel" dimension.
data_masked_snr_2, masked_elements_2 = xrutils.apply_mask(rec["amp"], snr_mask, "drop", "channel")

# show some results
print(f"channels that were masked according to the SNR threshold: {masked_elements_2}")

# dropped:
data_masked_snr_2

channels that were masked according to the SNR threshold: ['S4D4' 'S5D7' 'S6D8' 'S8D8']

[3]:

<xarray.DataArray (channel: 24, wavelength: 2, time: 23239)> Size: 9MB
<Quantity([[[0.0913686 0.0909875 0.0910225 ... 0.0941083 0.0940129 0.0944882]
  [0.1856806 0.186377  0.1836514 ... 0.1856486 0.1850836 0.1842172]]

 [[0.227516  0.2297024 0.2261366 ... 0.2264519 0.2271665 0.226713 ]
  [0.6354927 0.637668  0.6298023 ... 0.6072068 0.6087293 0.6091066]]

 [[0.1064704 0.1066212 0.1053444 ... 0.121114  0.1205022 0.1205441]
  [0.2755033 0.2761615 0.2727006 ... 0.2911952 0.2900544 0.2909847]]

 ...

 [[0.187484  0.1868235 0.1866562 ... 0.1735965 0.1736705 0.1738339]
  [0.2424386 0.241503  0.2408491 ... 0.22303   0.2229887 0.2234081]]

 [[0.2027881 0.1996586 0.2004866 ... 0.2318743 0.2311941 0.2330808]
  [0.4666358 0.4554404 0.4561614 ... 0.4809749 0.4812827 0.4862896]]

 [[0.6304559 0.6284427 0.6287045 ... 0.6810626 0.6809573 0.6818709]
  [1.2285622 1.2205907 1.2190002 ... 1.2729124 1.2727222 1.2755645]]], 'volt')>
Coordinates:
  * time        (time) float64 186kB 0.0 0.128 0.256 ... 2.974e+03 2.974e+03
    samples     (time) int64 186kB 0 1 2 3 4 5 ... 23234 23235 23236 23237 23238
  * channel     (channel) object 192B 'S1D1' 'S1D2' 'S1D3' ... 'S8D7' 'S8D16'
    source      (channel) object 192B 'S1' 'S1' 'S1' 'S1' ... 'S7' 'S8' 'S8'
    detector    (channel) object 192B 'D1' 'D2' 'D3' 'D9' ... 'D15' 'D7' 'D16'
  * wavelength  (wavelength) float64 16B 760.0 850.0
Attributes:
    data_type_group:  unprocessed raw

[4]:

# we can plot the values per channel in a 2D montage
wl_idx = 0
# plots._simple_scalp_plot(rec['amp'], rec.geo3d, snr.isel(wavelength=wl_idx), ax, cmap=p.cm.jet, vmin=0, vmax=25, title='SNR ' + str(rec.wavelengths[wl_idx]))
# plots.scalp_plot(rec, snr.isel(wavelength=wl_idx).values, ax, colormap=p.cm.jet, vmin=0, vmax=25, title='SNR ' + str(rec.wavelengths[wl_idx]))

f, ax = p.subplots(1, 2, figsize=(10, 4))
for i, wl in enumerate(rec["amp"].wavelength.values):
    plots.scalp_plot(
        rec["amp"],
        rec.geo3d,
        snr.sel(wavelength=wl),
        ax[i],
        cmap="RdYlGn",
        vmin=0,
        vmax=50,
        title=f"{wl} nm",
        cb_label="SNR",
        channel_lw=2
    )
f.tight_layout()

../../_images/examples_signal_quality_21_data_quality_and_pruning_7_0.png

Channel Quality Metrics: Scalp Coupling Index

[5]:

# Here we can assess the scalp coupling index (SCI) of the channels
sci_threshold = 0.75
window_length = 5*units.s
sci, sci_mask = quality.sci(rec["amp"], window_length, sci_threshold)

[6]:

# we can plot a heat map to visualize the SCI across all time windows and channels
colors = ["black", "#DC3220", "#5D3A9B", "#0C7BDC"]
nodes = [0.0, 0.75, 0.751, 1.0]
sci_cmap = LinearSegmentedColormap.from_list("sci_cmap", list(zip(nodes,colors)))
sci_binary_cmap = LinearSegmentedColormap.from_list("sci_binary_cmap", list(zip([0,0.5,0.5,1],["#DC3220","#DC3220","#0C7BDC","#0C7BDC"])))

# plot the heatmap
f,ax = p.subplots(1,1,figsize=(17,8))

m = ax.pcolormesh(sci.time, np.arange(len(sci.channel)), sci, shading="nearest", cmap=sci_cmap)
cb = p.colorbar(m, ax=ax)
cb.set_label("SCI")
ax.set_xlabel("time / s")
p.tight_layout()
ax.yaxis.set_ticks(np.arange(len(sci.channel)))
ax.yaxis.set_ticklabels(sci.channel.values);

# plot the binary heatmap
f,ax = p.subplots(1,1,figsize=(17,8))

m = ax.pcolormesh(sci.time, np.arange(len(sci.channel)), sci>0.75, shading="nearest", cmap=sci_binary_cmap)
cb = p.colorbar(m, ax=ax)
p.tight_layout()
ax.yaxis.set_ticks(np.arange(len(sci.channel)))
ax.yaxis.set_ticklabels(sci.channel.values);
cb.set_label("SCI > 0.75")
ax.set_xlabel("time / s");

/home/runner/miniconda3/envs/cedalion/lib/python3.11/site-packages/xarray/core/variable.py:338: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  data = np.asarray(data)

../../_images/examples_signal_quality_21_data_quality_and_pruning_10_1.png

../../_images/examples_signal_quality_21_data_quality_and_pruning_10_2.png

[7]:

ch = "S3D3"
t1,t2 = 100, 1200
f,ax = p.subplots(2,1, figsize=(17,5), sharex=True)
ax[0].plot(rec["amp"].time, rec["amp"].sel(channel=ch, wavelength=760), "r-")
ax[0].set_ylabel("HbO", color="r")
ax02 = ax[0].twinx()
ax02.plot(rec["amp"].time, rec["amp"].sel(channel=ch, wavelength=850), "b-")
ax02.set_ylabel("HbR", color="b")
ax[1].plot(sci.time, sci.sel(channel=ch), "m.-")
ax[1].set_ylabel("sci", color="m")
ax12 = ax[1].twinx()
ax12.plot(sci_mask.time, sci_mask.sel(channel=ch), "g-")
ax12.set_ylabel("sci_mask", color="g")
ax[1].set_xlim(t1, t2)
ax[1].axhline(0.75, c="k", ls=":")
for i in sci.time.values:
    if i < t1 or i > t2:
        continue
    ax[0].axvline(i, c="k", alpha=.2)
    ax[1].axvline(i, c="k", alpha=.2)

../../_images/examples_signal_quality_21_data_quality_and_pruning_11_0.png

Channel Quality Metrics: Peak Spectral Power

[8]:

# We can also look at the peak spectral power which takes the peak power of the cross-correlation signal between the cardiac band of the two wavelengths
psp_threshold = 0.1
psp, psp_mask = quality.psp(rec["amp"], window_length, psp_threshold)

[9]:

# We can look at similar heatmaps across time and channels
# plot the heatmap
colors = ["black", "#DC3220", "#5D3A9B", "#0C7BDC"]
nodes = [0.0, 0.45, 0.46, 1.0]
psp_cmap = LinearSegmentedColormap.from_list("psp_cmap", list(zip(nodes,colors)))
psp_binary_cmap = LinearSegmentedColormap.from_list("psp_binary_cmap", list(zip([0,0.5,0.5,1],["#DC3220","#DC3220","#0C7BDC","#0C7BDC"])))

f,ax = p.subplots(1,1,figsize=(17,8))

m = ax.pcolormesh(psp.time, np.arange(len(psp.channel)), psp, shading="nearest", cmap=psp_cmap)
cb = p.colorbar(m, ax=ax)
cb.set_label("PSP")
ax.set_xlabel("time / s")
p.tight_layout()
ax.yaxis.set_ticks(np.arange(len(psp.channel)))
ax.yaxis.set_ticklabels(psp.channel.values);

# plot the binary heatmap
f,ax = p.subplots(1,1,figsize=(17,8))

m = ax.pcolormesh(psp.time, np.arange(len(psp.channel)), psp>psp_threshold, shading="nearest", cmap=psp_binary_cmap)
cb = p.colorbar(m, ax=ax)
p.tight_layout()
ax.yaxis.set_ticks(np.arange(len(psp.channel)))
ax.yaxis.set_ticklabels(psp.channel.values);
cb.set_label("PSP > 0.1")
ax.set_xlabel("time / s");

../../_images/examples_signal_quality_21_data_quality_and_pruning_14_0.png

../../_images/examples_signal_quality_21_data_quality_and_pruning_14_1.png

Combining SCI and PSP

We want both SCI and PSP to be above their respective thresholds for a window to be considered clean. We can then look at the percentage of time both metrics are above the thresholds.

[10]:

sci_x_psp_mask = sci_mask & psp_mask
perc_time_clean = sci_x_psp_mask.sum(dim="time") / len(sci.time)

# plot the percentage of clean time per channel
f,ax = p.subplots(1,1,figsize=(8,8))
# plots.scalp_plot(rec, perc_time_clean, ax, colormap=p.cm.jet, vmin=0, vmax=1, title='Percentage of clean time')
# plots._simple_scalp_plot(rec['amp'], rec.geo3d, perc_time_clean, ax, cmap=p.cm.jet, vmin=0, vmax=1, title='Percentage of clean time')

plots.scalp_plot(
    rec["amp"],
    rec.geo3d,
    perc_time_clean,
    ax,
    cmap="RdYlGn",
    vmin=0,
    vmax=1,
    title=None,
    cb_label="Percentage of clean time",
    channel_lw=2
)
f.tight_layout()

# we can also plot this as a binary heatmap
f,ax = p.subplots(1,1,figsize=(17,8))

m = ax.pcolormesh(sci_x_psp_mask.time, np.arange(len(sci_x_psp_mask.channel)), sci_x_psp_mask, shading="nearest", cmap=sci_binary_cmap)
cb = p.colorbar(m, ax=ax)
p.tight_layout()
ax.yaxis.set_ticks(np.arange(len(sci_x_psp_mask.channel)))
ax.yaxis.set_ticklabels(sci_x_psp_mask.channel.values);
cb.set_label("PSP > 0.1 and SCI > 0.75")
ax.set_xlabel("time / s");

../../_images/examples_signal_quality_21_data_quality_and_pruning_16_0.png

../../_images/examples_signal_quality_21_data_quality_and_pruning_16_1.png

Channel Quality Metrics: Channel Distance

[11]:

# Here we assess channel distances. We might want to exclude very short or very long channels
sd_threshs = [1, 4.5]*units.cm # defines the lower and upper bounds for the source-detector separation that we would like to keep

# Source Detector Separation thresholding
ch_dist, sd_mask = quality.sd_dist(rec["amp"], rec.geo3d, sd_threshs)
# print the channel distances
print(f"channel distances are: {ch_dist}")

# apply mask function. In this example, we want to "drop" all channels that do not fall inside sd_threshs
# i.e. drop channels shorter than 1cm and longer than 4.5cm. We want to collapse along the "channel" dimension.
data_masked_sd, masked_elements = xrutils.apply_mask(rec["amp"], sd_mask, "drop", "channel")

# display the resultings
print(f"channels that were masked according to the SD Distance thresholds: {masked_elements}")
data_masked_sd

channel distances are: <xarray.DataArray (channel: 28)> Size: 224B
<Quantity([0.039 0.039 0.041 0.008 0.037 0.038 0.037 0.007 0.04  0.037 0.008 0.041
 0.034 0.008 0.039 0.039 0.041 0.008 0.037 0.037 0.037 0.008 0.04  0.037
 0.007 0.041 0.033 0.008], 'meter')>
Coordinates:
  * channel   (channel) object 224B 'S1D1' 'S1D2' 'S1D3' ... 'S8D8' 'S8D16'
    source    (channel) object 224B 'S1' 'S1' 'S1' 'S1' ... 'S7' 'S8' 'S8' 'S8'
    detector  (channel) object 224B 'D1' 'D2' 'D3' 'D9' ... 'D7' 'D8' 'D16'
channels that were masked according to the SD Distance thresholds: ['S1D9' 'S2D10' 'S3D11' 'S4D12' 'S5D13' 'S6D14' 'S7D15' 'S8D16']

[11]:

<xarray.DataArray (channel: 20, wavelength: 2, time: 23239)> Size: 7MB
<Quantity([[[0.0913686 0.0909875 0.0910225 ... 0.0941083 0.0940129 0.0944882]
  [0.1856806 0.186377  0.1836514 ... 0.1856486 0.1850836 0.1842172]]

 [[0.227516  0.2297024 0.2261366 ... 0.2264519 0.2271665 0.226713 ]
  [0.6354927 0.637668  0.6298023 ... 0.6072068 0.6087293 0.6091066]]

 [[0.1064704 0.1066212 0.1053444 ... 0.121114  0.1205022 0.1205441]
  [0.2755033 0.2761615 0.2727006 ... 0.2911952 0.2900544 0.2909847]]

 ...

 [[0.2225884 0.2187791 0.2195495 ... 0.2564863 0.2551258 0.2560233]
  [0.3994258 0.3917637 0.389261  ... 0.4304597 0.430814  0.4331249]]

 [[0.2027881 0.1996586 0.2004866 ... 0.2318743 0.2311941 0.2330808]
  [0.4666358 0.4554404 0.4561614 ... 0.4809749 0.4812827 0.4862896]]

 [[0.4885007 0.4802285 0.4818338 ... 0.6109142 0.6108118 0.613845 ]
  [0.8457658 0.825988  0.8259648 ... 0.975894  0.9756599 0.9826459]]], 'volt')>
Coordinates:
  * time        (time) float64 186kB 0.0 0.128 0.256 ... 2.974e+03 2.974e+03
    samples     (time) int64 186kB 0 1 2 3 4 5 ... 23234 23235 23236 23237 23238
  * channel     (channel) object 160B 'S1D1' 'S1D2' 'S1D3' ... 'S8D7' 'S8D8'
    source      (channel) object 160B 'S1' 'S1' 'S1' 'S2' ... 'S7' 'S8' 'S8'
    detector    (channel) object 160B 'D1' 'D2' 'D3' 'D1' ... 'D7' 'D7' 'D8'
  * wavelength  (wavelength) float64 16B 760.0 850.0
Attributes:
    data_type_group:  unprocessed raw

Channel Quality Metrics: Mean Amplitudes

[12]:

# Here we assess average channel amplitudes. We might want to exclude very small or large signals
amp_threshs = [0.1, 3]*units.volt # define whether a channel's amplitude is within a certain range

# Amplitude thresholding
mean_amp, amp_mask = quality.mean_amp(rec["amp"], amp_threshs)

# apply mask function. In this example, we want drop all channels that do not fall inside the amplitude thresholds.
# We collapse to the "channel" dimension.
data_masked_amp, masked_elements = xrutils.apply_mask(rec["amp"], amp_mask, "drop", "channel")

# display the results
print(f"channels that were masked according to the amplitude threshold: {masked_elements}")
data_masked_amp

channels that were masked according to the amplitude threshold: ['S1D1' 'S1D9' 'S3D2' 'S6D8' 'S7D6']

[12]:

<xarray.DataArray (channel: 23, wavelength: 2, time: 23239)> Size: 9MB
<Quantity([[[0.227516  0.2297024 0.2261366 ... 0.2264519 0.2271665 0.226713 ]
  [0.6354927 0.637668  0.6298023 ... 0.6072068 0.6087293 0.6091066]]

 [[0.1064704 0.1066212 0.1053444 ... 0.121114  0.1205022 0.1205441]
  [0.2755033 0.2761615 0.2727006 ... 0.2911952 0.2900544 0.2909847]]

 [[0.5512474 0.5510672 0.5476283 ... 0.6179242 0.6188702 0.6187721]
  [1.125532  1.1238331 1.1119423 ... 1.1817728 1.1819598 1.1832658]]

 ...

 [[0.2027881 0.1996586 0.2004866 ... 0.2318743 0.2311941 0.2330808]
  [0.4666358 0.4554404 0.4561614 ... 0.4809749 0.4812827 0.4862896]]

 [[0.4885007 0.4802285 0.4818338 ... 0.6109142 0.6108118 0.613845 ]
  [0.8457658 0.825988  0.8259648 ... 0.975894  0.9756599 0.9826459]]

 [[0.6304559 0.6284427 0.6287045 ... 0.6810626 0.6809573 0.6818709]
  [1.2285622 1.2205907 1.2190002 ... 1.2729124 1.2727222 1.2755645]]], 'volt')>
Coordinates:
  * time        (time) float64 186kB 0.0 0.128 0.256 ... 2.974e+03 2.974e+03
    samples     (time) int64 186kB 0 1 2 3 4 5 ... 23234 23235 23236 23237 23238
  * channel     (channel) object 184B 'S1D2' 'S1D3' 'S2D1' ... 'S8D8' 'S8D16'
    source      (channel) object 184B 'S1' 'S1' 'S2' 'S2' ... 'S8' 'S8' 'S8'
    detector    (channel) object 184B 'D2' 'D3' 'D1' 'D3' ... 'D7' 'D8' 'D16'
  * wavelength  (wavelength) float64 16B 760.0 850.0
Attributes:
    data_type_group:  unprocessed raw

Channel Pruning using Quality Metrics and the Pruning Function

To prune channels according to quality criteria, we do not have to manually go through the steps above. Instead, we can create quality masks for the metrics that we are interested in and hand them to a dedicated channel pruning function. The prune function expects a list of quality masks alongside a logical operator that defines how these masks should be combined.

[13]:

# as above we use three metrics and define thresholds accordingly
snr_thresh = 16 # the SNR (std/mean) of a channel.
sd_threshs = [1, 4.5]*units.cm # defines the lower and upper bounds for the source-detector separation that we would like to keep
amp_threshs = [0.1, 3]*units.volt # define whether a channel's amplitude is within a certain range

# then we calculate the masks for each metric: SNR, SD distance and mean amplitude
_, snr_mask = quality.snr(rec["amp"], snr_thresh)
_, sd_mask = quality.sd_dist(rec["amp"], rec.geo3d, sd_threshs)
_, amp_mask = quality.mean_amp(rec["amp"], amp_threshs)
# you can also include other masks, e.g. the SCI mask

# put all masks in a list
masks = [snr_mask, sd_mask, amp_mask]

# prune channels using the masks and the operator "all", which will keep only channels that pass all three metrics
amp_pruned, drop_list = quality.prune_ch(rec["amp"], masks, "all")

# print list of dropped channels
print(f"List of pruned channels: {drop_list}")
# display the new data xarray
amp_pruned

List of pruned channels: ['S1D1' 'S1D9' 'S2D10' 'S3D2' 'S3D11' 'S4D4' 'S4D12' 'S5D7' 'S5D13' 'S6D8'
 'S6D14' 'S7D6' 'S7D15' 'S8D8' 'S8D16']

[13]:

<xarray.DataArray (channel: 13, wavelength: 2, time: 23239)> Size: 5MB
<Quantity([[[0.227516  0.2297024 0.2261366 ... 0.2264519 0.2271665 0.226713 ]
  [0.6354927 0.637668  0.6298023 ... 0.6072068 0.6087293 0.6091066]]

 [[0.1064704 0.1066212 0.1053444 ... 0.121114  0.1205022 0.1205441]
  [0.2755033 0.2761615 0.2727006 ... 0.2911952 0.2900544 0.2909847]]

 [[0.5512474 0.5510672 0.5476283 ... 0.6179242 0.6188702 0.6187721]
  [1.125532  1.1238331 1.1119423 ... 1.1817728 1.1819598 1.1832658]]

 ...

 [[0.3463254 0.3424951 0.3408207 ... 0.3929267 0.3941368 0.3945422]
  [0.6978315 0.6875081 0.6857653 ... 0.7259991 0.7271688 0.7292138]]

 [[0.2225884 0.2187791 0.2195495 ... 0.2564863 0.2551258 0.2560233]
  [0.3994258 0.3917637 0.389261  ... 0.4304597 0.430814  0.4331249]]

 [[0.2027881 0.1996586 0.2004866 ... 0.2318743 0.2311941 0.2330808]
  [0.4666358 0.4554404 0.4561614 ... 0.4809749 0.4812827 0.4862896]]], 'volt')>
Coordinates:
  * time        (time) float64 186kB 0.0 0.128 0.256 ... 2.974e+03 2.974e+03
    samples     (time) int64 186kB 0 1 2 3 4 5 ... 23234 23235 23236 23237 23238
  * channel     (channel) object 104B 'S1D2' 'S1D3' 'S2D1' ... 'S7D7' 'S8D7'
    source      (channel) object 104B 'S1' 'S1' 'S2' 'S2' ... 'S6' 'S7' 'S8'
    detector    (channel) object 104B 'D2' 'D3' 'D1' 'D3' ... 'D7' 'D7' 'D7'
  * wavelength  (wavelength) float64 16B 760.0 850.0
Attributes:
    data_type_group:  unprocessed raw

Timeseries Quality Metric: Global Variance of the Temporal Derivative

[14]:

# we can plot the timeseries of the GVTD to evaluate motion in the data
gvtd, gvtd_mask = quality.gvtd(rec["amp"], stat_type="histogram_mode")

f, ax = p.subplots(2,1, figsize=(12,4), sharex=True)
ax[0].plot(gvtd.time, gvtd)
ax[1].plot(gvtd.time, gvtd_mask)
ax[0].set_xlabel("time / s")
ax[0].set_ylabel("GVTD / $s^{-1}$")
ax[1].set_ylabel("GVTD mask")
ax[1].set_yticks([0,1])
ax[1].set_yticklabels(["TAINTED","CLEAN"])
ax[0].set_xlim(0, 250)
ax[0].grid()
ax[1].grid()
f.suptitle("Global Variance of Temporal Derivative")
p.tight_layout()

../../_images/examples_signal_quality_21_data_quality_and_pruning_24_0.png