Calculating the Scalp Coupling Index
This notebook calculates the Scalp Coupling Index[1] metric for assessing the signal quality of a recording.
[1] L. Pollonini, C. Olds, H. Abaya, H. Bortfeld, M. S. Beauchamp, and J. S. Oghalai, “Auditory cortex activation to natural speech and simulated cochlear implant speech measured with functional near-infrared spectroscopy,” Hearing Research, vol. 309, pp. 84–93, Mar. 2014, doi: 10.1016/j.heares.2013.11.007.
[1]:
import cedalion
import cedalion.nirs
import cedalion.xrutils as xrutils
from cedalion.datasets import get_fingertapping
import numpy as np
import xarray as xr
import pint
import matplotlib.pyplot as p
import scipy.signal
import os.path
xr.set_options(display_expand_data=False)
[1]:
<xarray.core.options.set_options at 0x7f75f50bf510>
Loading raw CW-NIRS data from a SNIRF file
This notebook uses a finger-tapping dataset in BIDS layout provided by Rob Luke. It can can be downloaded via cedalion.datasets.
[2]:
elements = get_fingertapping()
amp = elements[0].data[0]
Calculating the SCI
From the paper:
Since the LED sources at 760 nm and 850 nm were co-located, an optical channel in good contact with the scalp exhibited a prominent synchronous cardiac pulsation in both photodetected signals. This observation was independent of the amplitude of the output voltage of the photodetector, which in turn depends on the inter-distance between sources and detector. For each channel, we filtered both photodetected signals between 0.5 and 2.5 Hz to preserve only the cardiac component and normalized the resulting signals to balance any difference between their amplitude. Then, we computed the cross-correlation and we extracted the value at a time lag of 0 to quantify the similarity between the filtered signals. In-phase and counter-phase identical waveforms yielded a zero-lag cross-correlation value of 1 and +1 respectively, whereas a null value derived from totally uncorrelated signals. Therefore, the zero-lag cross-correlation between photodetected signals of the same channel was used as a quantitative measure of the signal-to-noise ratio of the channel. We termed this value the scalp coupling index (SCI).
0. Utilities
[3]:
def plot_channel(array, channel, ylabel, xlabel="time", tmin=1000, tmax=1030):
f, ax = p.subplots(1,1, figsize=(12,4))
ax.plot(array.time, array.sel(channel=channel, wavelength=760), "r-")
ax.plot(array.time, array.sel(channel=channel, wavelength=850), "b-")
p.xlim(tmin, tmax)
p.xlabel(xlabel)
p.ylabel(ylabel)
1. Bandpass filter to extract the cardiac signal
[4]:
amp_filtered = amp.cd.freq_filter(0.5, 2.5, butter_order=4)
plot_channel(amp_filtered, "S5D7", "amplitude / V")
2. Normalize filtered amplitudes
Subtract the mean and normalize to each channels standard deviation.
[5]:
amp_filtered_normed = (amp_filtered - amp_filtered.mean("time")) / amp_filtered.std("time")
#amp_filtered_normed = (amp_filtered - amp_filtered.min("time")) / (amp_filtered.max("time") - amp_filtered.min("time"))
plot_channel(amp_filtered_normed, "S5D7", "normalized amplitude")
3. Moving windows
Calculate non-overlapping, moving windows of 5 seconds
[6]:
window_len_s = 5 # seconds
window_len_samples = int(np.ceil(window_len_s * amp_filtered_normed.cd.sampling_rate))
print(f"At a sampling rate of {amp_filtered_normed.cd.sampling_rate:.2f} Hz a {window_len_s} second window is {window_len_samples} samples long.")
At a sampling rate of 7.81 Hz a 5 second window is 40 samples long.
[7]:
# This creates a new DataArray with a new dimension "window", that is window_len_samples large.
# The time dimension will contain the time coordinate of the first sample in the window.
# Setting the stride size to the same value as the window length will result in non-overlapping windows.
windows = amp_filtered_normed.rolling(time=window_len_samples).construct("window", stride=window_len_samples)
display(windows)
f,ax = p.subplots(1,1, figsize=(12,2))
p.plot(amp_filtered_normed.time, np.ones(len(amp_filtered_normed.time)), "r|", label="amp_filtered_normed.time")
p.plot(windows.time, np.ones(len(windows.time)), "ks", label="windows.time")
p.xlim(0,40)
p.legend();
<xarray.DataArray (channel: 28, wavelength: 2, time: 581, window: 40)> Size: 10MB
[] nan nan nan nan nan nan nan ... -1.14 -0.4076 1.104 1.473 0.03794 -1.135
Coordinates:
* time (time) float64 5kB 0.0 5.12 10.24 ... 2.964e+03 2.97e+03
samples (time) int64 5kB 0 40 80 120 160 ... 23080 23120 23160 23200
* 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
Dimensions without coordinates: window
4. Calculate the correlation coefficient for each window
The cross-correlation of two time series \(X\) and \(Y\) at time lag \(\tau\) is:
At time lag \(\tau=0\) this reduces to:
This is here computed over the time window $[t_1, t_2] $ of length \(N\). The standard deviations \(\sigma_X\) and \(\sigma_X\) are calculated over the same time windows. The time series \(X\) and \(Y\) denote the two different wavelengths.
[8]:
sci = (windows - windows.mean("window")).prod("wavelength").sum("window") / window_len_samples
sci /= windows.std("window").prod("wavelength")
display(sci)
<xarray.DataArray (channel: 28, time: 581)> Size: 130kB
[] inf 0.6894 0.6215 0.4492 0.5985 0.673 ... 0.9552 0.9708 0.9479 0.9339 0.9169
Coordinates:
* time (time) float64 5kB 0.0 5.12 10.24 ... 2.959e+03 2.964e+03 2.97e+03
samples (time) int64 5kB 0 40 80 120 160 ... 23040 23080 23120 23160 23200
* 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'5. Illustrate heat maps of SCIs for the whole recording and all channels
[9]:
from matplotlib.colors import LinearSegmentedColormap
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"])))
[10]:
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);
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/micromamba/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)
6. Inspect time courses of good and bad channels
S1D1: SCI < 0.75 of the times
S3D3: SCI < 0.75 around t=2000s
S3D11: SCI < 0.75 around t=2000s
S6D5: SCI > 0.75 for all samples
[11]:
for channel in ["S1D1", "S3D3", "S3D11", "S6D5"]:
tmin, tmax = 2000, 2100
f,ax = p.subplots(2,1, figsize=(18,4), sharex=True)
m = (tmin <= amp.time) & (amp.time <= tmax)
ax[0].plot(amp.time[m], amp.sel(channel=channel, wavelength=760, time=m), "r-", alpha=.5)
ax[0].set_ylabel("amp. 760nm / V", color="r")
ax2 =ax[0].twinx()
ax2.plot(amp.time[m], amp.sel(channel=channel, wavelength=850, time = m), "b-", alpha=.5)
ax2.set_ylabel("amp. 850nm / V", color="b")
m = (tmin <= sci.time) & (sci.time <= tmax)
ax[1].scatter(sci.time[m], sci.sel(channel=channel, time=m), c=sci_cmap(sci.sel(channel=channel, time=m)))
ax[1].set_ylabel("SCI")
ax[1].set_xlabel("time / s")
ax[1].axhline(0.75, c="k", ls="--")
ax[1].set_ylim(0,1)
ax[0].grid(1)
ax[1].grid(axis="x")
f.suptitle(channel)
f.set_tight_layout(True)
/home/runner/micromamba/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)
7. Calculate a quality mask for each sample of the recording
[12]:
# sci.time coordinates contain the time of the first sample of each window.
# Use these as bin edges and calculate for each sample to which window it belongs.
# Subtract one from the indices returned by np.digitize as 0 denotes the underflow bin.
window_indices = np.digitize(amp.time, sci.time) - 1
# To obtain a quality mask for each sample we can threshold the sci array and then
# inflate it along the time dimension from n_windows values to n_samples:
qmask = (sci > 0.75)[:, window_indices]
qmask["time"] = amp.time # carry over time coordinates from original array
print(f'raw amplitude array has {len(amp.time)} values in the time dimension.')
print(f'SCI has {len(sci.time)} values in the time dimension.')
print(f'qmask has {len(qmask.time)} values in the time dimension.')
f,ax = p.subplots(1,1,figsize=(17,8))
m = ax.pcolormesh(
qmask.time,
np.arange(len(qmask.channel)),
qmask,
shading="nearest",
cmap=sci_binary_cmap,
edgecolors="w",
linewidths=0.5)
cb = p.colorbar(m, ax=ax)
cb.set_label("SCI > 0.75")
ax.set_xlabel("time / s")
p.tight_layout()
ax.yaxis.set_ticks(np.arange(len(qmask.channel)))
ax.yaxis.set_ticklabels(qmask.channel.values);
ax.set_xlim(500,520)
raw amplitude array has 23239 values in the time dimension.
SCI has 581 values in the time dimension.
qmask has 23239 values in the time dimension.
[12]:
(500.0, 520.0)
