Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

VIZ/DOC/ENH: add MNE rest example + smaller fixes #61

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion tutorials/plot_compute_evoked_inverse_solution.py
Original file line number Diff line number Diff line change
@@ -31,7 +31,7 @@
# These are the same as in :ref:`tut_plot_evoked`

hcp_evokeds = hcp.read_evokeds(onset='stim', subject=subject,
data_type=data_type, hcp_path=hcp_path)
data_type=data_type, hcp_path=hcp_path)
for evoked in hcp_evokeds:
if not evoked.comment == 'Wrkmem_LM-TIM-face_BT-diff_MODE-mag':
continue
199 changes: 199 additions & 0 deletions tutorials/plot_compute_resting_inverse_solution.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,199 @@
"""
.. _tut_compute_inverse_erf:

========================================
Compute inverse solution for evoked data
========================================

Here we'll use our knowledge from the other examples and tutorials
to compute an inverse solution and apply it on resting state data.
"""
# Author: Denis A. Engemann
# Luke Bloy <luke.bloy@gmail.com>
# Eric Larson <larson.eric.d@gmail.com>
# License: BSD 3 clause

import os.path as op
import mne
from mne.filter import next_fast_len
import hcp
from hcp import preprocessing as preproc

##############################################################################
# we assume our data is inside a designated folder under $HOME
storage_dir = op.expanduser('~/mne-hcp-data')
hcp_path = op.join(storage_dir, 'HCP')
recordings_path = op.join(storage_dir, 'hcp-meg')
subjects_dir = op.join(storage_dir, 'hcp-subjects')
subject = '105923' # our test subject
data_type = 'rest'
run_index = 0

##############################################################################
# We're reading the evoked data.
# These are the same as in :ref:`tut_plot_evoked`

raw = hcp.read_raw(subject=subject,
data_type=data_type, hcp_path=hcp_path,
run_index=run_index)
raw.load_data()
raw.crop(0, 250)
raw.resample(100.)
raw.info['bads'] = ['A147'] # we know this one is bad here.
preproc.set_eog_ecg_channels(raw)
preproc.apply_ref_correction(raw)

##############################################################################
# Remap BTi sensor coordinates to Neuromag for plotting.
info = raw.info.copy()
preproc.map_ch_coords_to_mne(raw) # we have to revert it later for inverse!


raw.info['projs'] = []
ecg_ave = mne.preprocessing.create_ecg_epochs(raw).average()
ecg_ave.plot_joint(times='peaks')
# We see no clear ECG artifact!

eog_ave = mne.preprocessing.create_eog_epochs(raw).average()
eog_ave.plot_joint(times=[0])
# But we see a clear EOG artifact. So we'll compute an SSP

ssp_eog, _ = mne.preprocessing.compute_proj_eog(
raw, n_grad=1, n_mag=1, average=True, reject=dict(mag=5e-12))
raw.add_proj(ssp_eog, remove_existing=True)

fig = mne.viz.plot_projs_topomap(raw.info['projs'],
info=raw.info)
fig.suptitle('HCP BTi/4D SSP')
fig.subplots_adjust(0.05, 0.05, 0.95, 0.85)

##############################################################################
# Explore data

n_fft = next_fast_len(int(round(4 * raw.info['sfreq'])))
print('Using n_fft=%d (%0.1f sec)' % (n_fft, n_fft / raw.info['sfreq']))
fig = raw.plot_psd(n_fft=n_fft, proj=True, fmax=50)
fig.subplots_adjust(0.1, 0.1, 0.95, 0.85)

# We see clear alpha and beta peaks!

##############################################################################
# We'll now use a convenience function to get our forward and source models
# instead of computing them by hand.

hcp.make_mne_anatomy(
subject, subjects_dir=storage_dir + '/hcp-subjects',
hcp_path=storage_dir + '/HCP',recordings_path=storage_dir + '/hcp-meg')

src_outputs = hcp.anatomy.compute_forward_stack(
subject=subject, subjects_dir=subjects_dir,
hcp_path=hcp_path, recordings_path=recordings_path,
# speed up computations here. Setting `add_dist` to True may improve the
# accuracy.
src_params=dict(add_dist=False),
info_from=dict(data_type=data_type, run_index=run_index))

fwd = src_outputs['fwd']

##############################################################################
# Now we can compute the noise covariance. For this purpose we will apply
# the same filtering as was used for the computations of the ERF in the first
# place. See also :ref:`tut_reproduce_erf`.

raw_noise = hcp.read_raw(subject=subject, hcp_path=hcp_path,
data_type='noise_empty_room')
raw_noise.load_data()

# apply ref channel correction, drop ref channels, add ssp.
preproc.apply_ref_correction(raw_noise)
raw_noise.add_proj(ssp_eog)


##############################################################################
# Note that using the empty room noise covariance will inflate the SNR of the
# evkoked and renders comparisons to `baseline` rather uninformative.
noise_cov = mne.compute_raw_covariance(raw_noise, method='oas')


##############################################################################
# Compute and apply inverse to PSD estimated using multitaper + Welch.
# Group into frequency bands, then normalize each source point and sensor
# independently. This makes the value of each sensor point and source location
# in each frequency band the percentage of the PSD accounted for by that band.

freq_bands = dict(
delta=(2, 4), theta=(5, 7), alpha=(8, 12), beta=(15, 29), gamma=(30, 45))
topos = dict()
stcs = dict()

##############################################################################
# redo original BTi sensor coordinates.
raw.info = info

inverse_operator = mne.minimum_norm.make_inverse_operator(
raw.info, forward=fwd, noise_cov=noise_cov, verbose=True)

stc_psd, ch_psd = mne.minimum_norm.compute_source_psd(
raw, inverse_operator, lambda2=1. / 9.,
n_fft=n_fft, dB=False, return_sensor=True, verbose=True,
fmax=50)

# remap sensors coords for power
preproc.map_ch_coords_to_mne(ch_psd)

stc_psd = stc_psd.to_original_src(
src_outputs['src_fsaverage'], subjects_dir=subjects_dir)

topo_norm = ch_psd.data.sum(axis=1, keepdims=True)
stc_norm = stc_psd.sum()
# Normalize each source point by the total power across freqs
for band, limits in freq_bands.items():
data = ch_psd.copy().crop(*limits).data.sum(axis=1, keepdims=True)
topos[band] = mne.EvokedArray(
100 * data / topo_norm, ch_psd.info)
stcs[band] = \
100 * stc_psd.copy().crop(*limits).sum() / stc_norm.data


###############################################################################
# Now we can make some plots of each frequency band. Note that the OPM head
# coverage is only over right motor cortex, so only localization
# of beta is likely to be worthwhile.
#
# Theta
# -----


def plot_band(band):
"""Plot channel and source band power."""
title = "%s\n(%d-%d Hz)" % ((band,) + freq_bands[band])
topos[band].plot_topomap(
times=0., scalings=1., cbar_fmt='%0.1f', vmin=0, cmap='inferno',
time_format=title)
brain = stcs[band].plot(
subject='fsaverage', subjects_dir=subjects_dir, views='cau', hemi='both',
time_label=title, title=title, colormap='inferno',
clim=dict(kind='percent', lims=(70, 85, 99)))
brain.show_view(dict(azimuth=0, elevation=0), roll=0)
return fig, brain


fig_theta, brain_theta = plot_band('theta')

###############################################################################
# Alpha
# -----

fig_alpha, brain_alpha = plot_band('alpha')

###############################################################################
# Beta
# ----

fig_beta, brain_beta = plot_band('beta')

###############################################################################
# Gamma
# -----

fig_gamma, brain_gamma = plot_band('gamma')