- MNE-Python version: 0.23.2 - 3.7.2
- operating system: Win 10 and WSL2
Hi, as I discuss during the office hours meeting, I made a request for the tutorial.
I will try to explain it as short as possible. and also provide necessary examples for it.
Example use data from this public dataset. The minimal example which I started to create is uploaded on GD.
This is a graph of the pipeline I had in mind.
Questions are asked in the comments, next to the lines of code they refer to.
Willing to avoid writing a long essay on this thread by doubling the questions.
I put the code also here to make it accessible for those who don’t want to download files from GD.
import mne
import os.path as op
import numpy as np
from mne.datasets import fetch_fsaverage
from mne.datasets import sample
from mne.minimum_norm import apply_inverse_raw, read_inverse_operator, make_inverse_operator, apply_inverse, write_inverse_operator
import matplotlib.pyplot as plt
# Download fsaverage files
fs_dir = fetch_fsaverage(verbose=True)
subjects_dir = op.dirname(fs_dir)
# The files live in:
subject = 'fsaverage'
trans = 'fsaverage' # MNE has a built-in fsaverage transformation
src = op.join(fs_dir, 'bem', 'fsaverage-ico-5-src.fif')
bem = op.join(fs_dir, 'bem', 'fsaverage-5120-5120-5120-bem-sol.fif')
fname_raw = 'eeg_data/sub-xp201_task-MIpre_eeg.vhdr'
# load data as raw_brainvision
rawD = mne.io.read_raw_brainvision(fname_raw,preload=True)
notches = np.arange(50,100)
rawD.notch_filter(notches, phase='zero-double', fir_design='firwin2')
# Filter section for gradient noise from scanner artefacts (not-implemented)
# expect the data to be super noisy
# construct data structure by hand (see why on 39 line)
dataD = rawD.get_data()
dataD = np.concatenate((dataD[0:30], dataD[32:65]), axis=0) # this is done to prevent error from line 39
names = np.concatenate((rawD.ch_names[0:30], rawD.ch_names[32:65]))
data = dataD
plt.plot(data[0])
plt.show()
info = mne.create_info(names.tolist(), 100, ch_types=['eeg']*62)
# uncomment line 43 to perform read_raw_brainvision data structure
# this gives an error:
# ValueError: DigMontage is only a subset of info. There are 1 channel position not present in the DigMontage. The required channels are: ['ECG'].
raw = mne.io.RawArray(data, info)
#raw = rawD
montage = mne.channels.make_standard_montage('standard_1020')
raw.set_eeg_reference(projection=True)
raw.set_montage(montage, match_case=False)
# compute forward solution
fwd = mne.make_forward_solution(raw.info, trans=trans, src=src, bem=bem, meg=False, eeg=True, mindist=5.0, n_jobs=1)
mne.viz.plot_alignment(
raw.info, src=src, eeg=['original','projected'], trans=trans,
show_axes=True, mri_fiducials=True, dig='fiducials')
# Use fwd to compute the sensitivity map for illustration purposes
eeg_map = mne.sensitivity_map(fwd, ch_type='eeg', mode='fixed')
brain = eeg_map.plot(time_label='EEG sensitivity', subjects_dir=subjects_dir,
clim=dict(lims=[5, 50, 100]))
# compute covariance raw covariance
noise_cov = mne.compute_raw_covariance(raw,tmin=-0.01)
figs = mne.viz.plot_cov(noise_cov, raw.info)
# compute inverse operator
inverse_operator = make_inverse_operator(info, fwd, noise_cov,
loose=0.2, depth=0.8)
# save the inverser operator to disk
fname_inv = 'eeg_data/sub-xp201_task-MIpost_eeg-inv.fif'
write_inverse_operator(fname_inv, inverse_operator)
# how to apply in to this pipeline the sLoreta as in this example:
# https://mne.tools/stable/auto_examples/inverse/compute_mne_inverse_raw_in_label.html?highlight=loreta
# this currently gives an error: for j, label in enumerate(labels))
# ZeroDivisionError: integer division or modulo by zero
# It would be helpful to have an example how to prepare parcellation like in this example: 'Aud-lh'
label_name = 'Aud-lh'
data_path = sample.data_path()
fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name
label = mne.read_label(fname_label)
start, stop = raw.time_as_index([0, 15]) # read the first 15s of data
snr = 1.0 # use smaller SNR for raw data
lambda2 = 1.0 / snr ** 2
method = "sLORETA" # use sLORETA method (could also be MNE or dSPM)
# Compute inverse solution
stc = apply_inverse_raw(raw, inverse_operator, lambda2, method, label,
start, stop, pick_ori=None)
plt.plot(1e3 * stc.times, stc.data[::100, :].T)
plt.xlabel('time (ms)')
plt.ylabel('%s value' % method)
plt.show()
Please give me feedback if something is unclear. 
Best,
Lukas