mne_exercise.py

#!/usr/bin/python

import sys
import numpy as np
import mne
from mne.preprocessing import ICA

argc = len(sys.argv)

# lists of data sources
subjects = ('CC110033', 'CC110037', 'CC110045')
# different origins
archs = ('passive', 'rest', 'task')

id_subject = 0
id_arch = 0
subject = subjects[id_subject]
arch = archs[id_arch]

subjects_dir = 'subjects'

# load raw data from the first subject
raw_data = mne.io.Raw(subject + '/' + arch + '/' + arch + '_raw.fif')

# display the channels. expert can select a channel to 'blacklist' it
if argc > 1 and sys.argv[1] == '-p':
    #raw_data.plot(block=True, lowpass=40)
    raw_data.plot_psd()

# filter the bad channels
calibration_file = 'sss_cal.dat'
cross_talk_file = 'ct_sparse.fif'
preprocessed_data = mne.preprocessing.maxwell_filter(raw_data,
        calibration=calibration_file,
        cross_talk=cross_talk_file)

# display the channels after filtering
if argc > 1 and sys.argv[1] == '-p':
    preprocessed_data.plot(block=True, lowpass=40)
    preprocessed_data.plot_psd()

# keep frequencies from 1Hz to 45Hz (remove noise and signal of non-interest)
low_freq = 1 
high_freq = 45
bandpass_data = preprocessed_data.filter(low_freq, high_freq)

# checked that the filtered frequencies are dumped
if argc > 1 and sys.argv[1] == '-p':
    #bandpass_data.plot(block=True, lowpass=40)
    bandpass_data.plot_psd()

# compute ICA
picks_meg = mne.pick_types(bandpass_data.info, meg=True, eeg=False, eog=False,
        stim=False, exclude='bads')

method = 'fastica'
# a good estimate of the number of components is the rank of the Raw object
n_components = bandpass_data.estimate_rank()
decim = 3
random_state = 23
ica = ICA(n_components=n_components, method=method, random_state=random_state)

reject = dict(mag=5e-12, grad=4000e-13)
# goal: remove artifacts, so they should be easy to notice
ica.fit(bandpass_data, picks=picks_meg, decim=decim, reject=reject)

if argc > 1 and sys.argv[1] == '-p':
    ica.plot_components()


# pass on Epoch computing
# 1) extract events from stim channels
events = mne.find_events(bandpass_data, stim_channel='STI101')
# 2) create epochs from Raw + events
epochs = mne.Epochs(bandpass_data, events)
# 3) average the epochs to reduce noise
evoked = epochs.average()

evoked.shift_time(-0.05)
if argc > 1 and sys.argv[1] == '-p':
    evoked.plot()


bandpass_data.save('out.fif')
mne.gui.coregistration(subjects_dir=subjects_dir)

mne_exercise2.py

#!/usr/bin/python

import sys
import numpy as np
import mne
from mne.preprocessing import ICA

from mayavi import mlab  # noqa
from surfer import Brain  # noqa

argc = len(sys.argv)

# lists of data sources
subjects = ('CC110033', 'CC110037', 'CC110045')
# different origins
archs = ('passive', 'rest', 'task')

id_subject = 0
id_arch = 0
subject = subjects[id_subject]
arch = archs[id_arch]

subjects_dir = 'subjects'

# input Raw file
raw_fname = subject + '/' + arch + '/' + arch + '_raw.fif'
# Raw file after processing
#raw_fname = 'out.fif'

# The transformation file obtained by coregistration
trans_fname = 'CC110033-trans.fif'  

bem_fname = subjects_dir + '/' + subject + '/bem/' + 'CC110033-meg-bem.fif'

info = mne.io.read_info(raw_fname)

mne.viz.plot_bem(subject=subject, subjects_dir=subjects_dir,
                 brain_surfaces='white', orientation='coronal')

mne.viz.plot_alignment(info, trans_fname, subject=subject, dig=True,
                       meg=['helmet', 'sensors'], subjects_dir=subjects_dir,
                       surfaces=['head', 'brain'])

src = mne.setup_source_space(subject, spacing='oct6',
                             subjects_dir=subjects_dir, add_dist=False)
print(src)

brain = Brain(subject, 'lh', 'inflated', subjects_dir=subjects_dir)
surf = brain.geo['lh']

vertidx = np.where(src[0]['inuse'])[0]
#vertidx = src[0]['inuse'][0]

mlab.points3d(surf.x[vertidx], surf.y[vertidx],
              surf.z[vertidx], color=(1, 0, 0), scale_factor=1.5)


# compute forward solution
fw = mne.make_forward_solution(raw_fname, trans=trans_fname, src=src, bem=bem_fname, meg=True, eeg=False)

print(fw)

mne_readme

MEG Exercise September 2018

* concept of maxfilter and bad channels:

some calibration files are needed for maxwell filter (found in parietal repo).
visual inspection is enough to determine bad channels.

NB:

* patient is passive when doing a passive activity and rest he does nothing at all 

* HPI coil - Head Position Indicator placed at known locations on the scalp
of the subject, when energized, will generate a magnetic field that helps us to localize the position of
head in a three-dimensional space.

* ICA - computational method for separating a multivariate signal into additive subcomponents. This is
done by assuming that the subcomponents are non-Gaussian signals and that they are statistically
independent from each other

** ICA is not mandatory for evoked datasets.

Overview:
RAW -> Epochs -> Evoked

stimuli are stored in stim_channel (i.e. 'STI_101')

mne.find_events() can extract the events from the stim_channel as an array and we can
create Epochs from the couple (Raw, events).

in the plot of the avegared epoch, the main peak should be around 100ms.
use shift_time to correct axis if it's not the case.

before using mne.viz.plot_aligmnent() don't forget to create an interactive python environment:
$ ipython