Dipole modeling and topomaps

simulate_dipoles_to_evoked.py

from os import path as op
import numpy as np
# from scipy import linalg

import mne

data_path = mne.datasets.sample.data_path()
subjects_dir = op.join(data_path, 'subjects')
fname_ave = op.join(data_path, 'MEG', 'sample', 'sample_audvis-ave.fif')
fname_cov = op.join(data_path, 'MEG', 'sample', 'sample_audvis-cov.fif')
fname_bem = op.join(subjects_dir, 'sample', 'bem',
                    'sample-5120-5120-5120-bem-sol.fif')
fname_trans = op.join(data_path, 'MEG', 'sample',
                      'sample_audvis_raw-trans.fif')

evoked = mne.read_evokeds(fname_ave, condition='Right Auditory',
                          baseline=(None, 0))
evoked.pick_types(meg=True, eeg=True)

info = evoked.info
cov = mne.read_cov(fname_cov)

# first time point has 2 (simultaneous) dipoles, second one is close to radial
# but second dipole is invalid (500 mm z-coord) and therefore omitted
times = np.array([0.0, 0.0, 0.001])  # first two, then one dipole
gof = 1.0
pos = np.array([[-0.060, -0.010, 0.060],
                [0.050, -0.010, 0.50],  # this dipole is INVALID (min_dist)
                [0.00, 0.050, 0.050]])  # (n_dipoles, 3) [m]
amp = np.array([20e-9, -20e-9, 40e-9])  # (n_dipoles,) [Am]
ori = np.array([[0., 1.0, 1.0],
                [0., 1.0, 1.0],
                [0., 1.0, 1.0]])  # (n_dipoles, 3)
ori_amp = np.diag(np.sqrt(np.dot(ori, ori.T)))
ori_norm = np.dot(np.diag(1./ori_amp), ori)

# create a Dipole, coordinates in HEAD frame
dipoles = mne.dipole.Dipole(times, pos, amp, ori_norm, gof)

# NB returning the stc and fwd means we don't have to recalculate when
# adding noise!!
stc, fwd = mne.simulation.simulate_stc_from_dipole(dipoles, fname_bem, info,
                                                   trans=fname_trans)

dipevo = mne.simulation.simulate_evoked(fwd, stc, info, None, snr=np.inf)
uniqts = dipevo.times
dipevo.plot_topomap(times=uniqts, ch_type='mag', outlines='skirt')
dipevo.plot_topomap(times=uniqts, ch_type='eeg', outlines='skirt')

# use covariance matrix to add a little noise (due to the way it's calculated
# the SNR here is not very realistic)
dipevo_noise = mne.simulation.simulate_evoked(fwd, stc, info,
                                              cov, snr=15.)
uniqts = dipevo_noise.times
dipevo_noise.plot_topomap(times=uniqts, ch_type='mag', outlines='skirt')
dipevo_noise.plot_topomap(times=uniqts, ch_type='eeg', outlines='skirt')