ERP Construction.

Example on how to interpolate missing channels.

import os

import matplotlib.pyplot as plt
import mne

import numpy as np

from eegrasp import EEGrasp

current_dir = os.getcwd()
os.chdir(os.path.dirname(current_dir))
data_dir = './datasets'
#os.makedirs("./datasets", exist_ok=True)
#os.environ['MNE_EEGBCI_PATH'] = './datasets/'

subjects = np.arange(1, 10)
runs = [4, 8, 12]

# Download eegbci dataset through MNE
# Comment the following line if already downloaded

raw_fnames = [
    mne.datasets.eegbci.load_data(s, runs, path=data_dir, update_path=True)
    for s in subjects
]
raw_fnames = np.reshape(raw_fnames, -1)
raws = [mne.io.read_raw_edf(f, preload=True) for f in raw_fnames]
raw = mne.concatenate_raws(raws)
mne.datasets.eegbci.standardize(raw)
raw.annotations.rename(dict(T1='left', T2='right'))

montage = mne.channels.make_standard_montage('standard_1005')
raw.set_montage(montage)
eeg_pos = np.array(
    [pos for _, pos in raw.get_montage().get_positions()['ch_pos'].items()])
ch_names = montage.ch_names

Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S001/S001R04.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S001/S001R08.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S001/S001R12.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S002/S002R04.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S002/S002R08.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S002/S002R12.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S003/S003R04.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S003/S003R08.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S003/S003R12.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S004/S004R04.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S004/S004R08.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S004/S004R12.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S005/S005R04.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S005/S005R08.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S005/S005R12.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S006/S006R04.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S006/S006R08.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S006/S006R12.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S007/S007R04.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S007/S007R08.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S007/S007R12.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S008/S008R04.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S008/S008R08.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S008/S008R12.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S009/S009R04.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S009/S009R08.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...
Extracting EDF parameters from /home/docs/checkouts/readthedocs.org/user_builds/eegrasp/checkouts/latest/examples/datasets/MNE-eegbci-data/files/eegmmidb/1.0.0/S009/S009R12.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19679  =      0.000 ...   122.994 secs...

L_FREQ = 1  # Hz
H_FREQ = 30  # Hz
raw.filter(L_FREQ, H_FREQ, fir_design='firwin', skip_by_annotation='edge')
raw, ref_data = mne.set_eeg_reference(raw)

events, events_id = mne.events_from_annotations(raw)

Filtering raw data in 27 contiguous segments
Setting up band-pass filter from 1 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 529 samples (3.306 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  64 out of  64 | elapsed:    0.0s finished
EEG channel type selected for re-referencing
Applying average reference.
Applying a custom ('EEG',) reference.
Used Annotations descriptions: [np.str_('T0'), np.str_('left'), np.str_('right')]

Exclude bad channels

TMIN, TMAX = -1.0, 3.0
picks = mne.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=False,
                       exclude='bads')
epochs = mne.Epochs(raw, events, events_id, picks=picks, tmin=TMIN, tmax=TMAX,
                    baseline=(-1, 0), detrend=1)

Not setting metadata
810 matching events found
Applying baseline correction (mode: mean)
0 projection items activated

left = epochs['left']
erp_left = left.average()

right = epochs['right']
erp_right = right.average()

# Use only data on the Left condition to find
# the best distance (epsilon) value
data = erp_left.get_data()

# 1. Define index of the missing channel
MISSING_IDX = 5
# 2. Initialize instance of EEGrasp
eegsp = EEGrasp(data, eeg_pos, ch_names)
# 3. Compute the electrode distance matrix
dist_mat = eegsp.compute_distance(normalize=True)
# 4. Find the best parameter for the channel
results = eegsp.fit_sigma(missing_idx=MISSING_IDX, distances=dist_mat, epsilon=0.5,
                          min_sigma=0.01, max_sigma=0.5, step=0.05)

  0%|          | 0/10 [00:00<?, ?it/s]
 10%|█         | 1/10 [00:00<00:01,  8.01it/s]
 20%|██        | 2/10 [00:00<00:00,  8.14it/s]
 30%|███       | 3/10 [00:00<00:00,  8.22it/s]
 40%|████      | 4/10 [00:00<00:00,  8.24it/s]
 50%|█████     | 5/10 [00:00<00:00,  8.25it/s]
 60%|██████    | 6/10 [00:00<00:00,  8.25it/s]
 70%|███████   | 7/10 [00:00<00:00,  8.25it/s]
 80%|████████  | 8/10 [00:00<00:00,  8.25it/s]
 90%|█████████ | 9/10 [00:01<00:00,  8.24it/s]
100%|██████████| 10/10 [00:01<00:00,  8.24it/s]
100%|██████████| 10/10 [00:01<00:00,  8.23it/s]

error = results['error']
best_idx = np.argmin(error[~np.isnan(error)])
reconstructed_signal = results['signal'][MISSING_IDX, :]
best_sigma = results['best_sigma']
vdistances = results['sigma']

plt.subplot(211)
plt.plot(vdistances, error, color='black')
plt.scatter(vdistances, error, color='teal', marker='x', alpha=0.5)
plt.scatter(best_sigma, error[vdistances == best_sigma], color='red')
plt.xlabel(r'$\sigma$')
plt.ylabel(r'RMSE')
plt.title('Error')

plt.subplot(212)
plt.title('Best Reconstruction')
plt.plot(reconstructed_signal, label='Reconstructed Data')
plt.plot(data[MISSING_IDX, :], label='Original Data')
plt.xlabel('samples')
plt.ylabel('Voltage')
plt.legend()

plt.tight_layout()
plt.show()

original = erp_right.get_data()
new_data = original.copy()
# Delete information from the missing channel
new_data[MISSING_IDX, :] = np.nan

# Compute best graph for interpolation
eegsp.compute_graph(epsilon=0.5, sigma=best_sigma)
# Interpolate channel
interpolated = eegsp.interpolate_channel(data=new_data, missing_idx=MISSING_IDX)

plt.plot(interpolated[MISSING_IDX, :], label='Interpolated Data', color='purple')
plt.plot(original[MISSING_IDX, :], label='Original Data', color='teal')
plt.xlabel('Samples')
plt.ylabel('Voltage')
plt.legend()
plt.show()

Estimated memory usage: 847 MB