In [1]:
import sys
print(sys.version_info)
del sys
sys.version_info(major=2, minor=7, micro=13, releaselevel='final', serial=0)

dSPM

準備

In [2]:
import numpy as np
import matplotlib.pyplot as plt
import mne
from mne.datasets import sample
from mne.minimum_norm import (make_inverse_operator,apply_inverse,write_inverse_operator)
mne.set_log_level('INFO')
data_path=sample.data_path()
raw_fname=data_path+'\\MEG\\sample\\sample_audvis_filt-0-40_raw.fif'
In [3]:
raw=mne.io.read_raw_fif(raw_fname)
raw.set_eeg_reference()
events=mne.find_events(raw,stim_channel='STI 014')

event_id=dict(aud_r=1)
tmin=-0.2
tmax=0.5
raw.info['bads']=['MEG 2443','EEG 053']
picks=mne.pick_types(raw.info,meg=True,eeg=False,eog=True,exclude='bads')
baseline=(None,0)
reject=dict(grad=4000e-13,mag=4e-12,eog=150e-6)
epochs=mne.Epochs(raw,events,event_id,tmin,tmax,proj=True,picks=picks,baseline=baseline,reject=reject)
Opening raw data file C:\Users\akira\mne_data\MNE-sample-data\MEG\sample\sample_audvis_filt-0-40_raw.fif...
    Read a total of 4 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
        Average EEG reference (1 x 60)  idle
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Current compensation grade : 0
An average reference projection was already added. The data has been left untouched.
<ipython-input-3-95821d46ba12>:2: RuntimeWarning: An average reference projection was already added. The data has been left untouched.
  raw.set_eeg_reference()
319 events found
Events id: [ 1  2  3  4  5 32]
72 matching events found
Created an SSP operator (subspace dimension = 3)
4 projection items activated

分散行列の計算

In [4]:
noise_cov=mne.compute_covariance(epochs,tmax=0.,method=['shrunk','empirical'])
print(noise_cov.data.shape)
fig_cov,fig_spectra=mne.viz.plot_cov(noise_cov,raw.info);
Loading data for 72 events and 106 original time points ...
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on MAG : [u'MEG 1711']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
17 bad epochs dropped
Estimating covariance using SHRUNK
Done.
Estimating covariance using EMPIRICAL
Done.
Using cross-validation to select the best estimator.
Number of samples used : 1705
[done]
Number of samples used : 1705
[done]
log-likelihood on unseen data (descending order):
   shrunk: -1480.993
   empirical: -1628.225
selecting best estimator: shrunk
(305L, 305L)

加算波形の表示

In [5]:
evoked=epochs.average()
evoked.plot()
evoked.plot_topomap(times=np.linspace(0.05,0.15,5),ch_type='mag');
evoked.plot_white(noise_cov);
estimated rank (grad): 203
estimated rank (mag): 102
estimated rank (mag + grad): 305
    Created an SSP operator (subspace dimension = 3)
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.

逆問題を解く準備

導体モデルは予めFreeSurferで計算したメッシュを使います。

In [6]:
# 導体モデル読み込み
fname_fwd=data_path+'\\MEG\\sample\\sample_audvis-meg-oct-6-fwd.fif'

# 順問題を解く
fwd=mne.read_forward_solution(fname_fwd,surf_ori=True)
print(fwd)
fwd=mne.pick_types_forward(fwd,meg=True,eeg=False)
print(fwd)
print(fwd.keys())

# 逆問題を解く
info=evoked.info
inverse_operator=make_inverse_operator(info,fwd,noise_cov,loose=0.2,depth=0.8)
print(inverse_operator.keys())
write_inverse_operator('sample_audvis-meg-oct-6-inv.fif',inverse_operator)
Reading forward solution from C:\Users\akira\mne_data\MNE-sample-data\MEG\sample\sample_audvis-meg-oct-6-fwd.fif...
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    2 source spaces read
    Desired named matrix (kind = 3523) not available
    Read MEG forward solution (7498 sources, 306 channels, free orientations)
    Source spaces transformed to the forward solution coordinate frame
    Converting to surface-based source orientations...
    Average patch normals will be employed in the rotation to the local surface coordinates....
[done]
<Forward | MEG channels: 306 | EEG channels: 0 | Source space: Surface with 7498 vertices | Source orientation: Free>
    306 out of 306 channels remain after picking
<Forward | MEG channels: 306 | EEG channels: 0 | Source space: Surface with 7498 vertices | Source orientation: Free>
['info', 'sol_grad', 'nchan', 'source_nn', 'sol', 'source_rr', 'source_ori', 'surf_ori', 'coord_frame', '_orig_sol', 'mri_head_t', 'nsource', 'src', '_orig_source_ori']
Computing inverse operator with 305 channels.
    Created an SSP operator (subspace dimension = 3)
estimated rank (mag + grad): 302
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Total rank is 302
Creating the depth weighting matrix...
    203 planar channels
    limit = 7265/7498 = 10.037795
    scale = 2.52065e-08 exp = 0.8
Computing inverse operator with 305 channels.
Creating the source covariance matrix
Applying loose dipole orientations. Loose value of 0.2.
Whitening the forward solution.
Adjusting source covariance matrix.
Computing SVD of whitened and weighted lead field matrix.
    largest singular value = 4.65276
    scaling factor to adjust the trace = 1.03619e+19
['nave', 'methods', 'eigen_fields', 'eigen_leads_weighted', 'fmri_prior', 'orient_prior', 'mri_head_t', 'sing', 'noise_cov', 'nsource', 'info', 'src', 'source_cov', 'projs', 'source_nn', 'depth_prior', 'source_ori', 'coord_frame', 'units', 'eigen_leads']
Write inverse operator decomposition in sample_audvis-meg-oct-6-inv.fif...
    Write a source space...
    [done]
    Write a source space...
    [done]
    2 source spaces written
    Writing inverse operator info...
    Writing noise covariance matrix.
    Writing source covariance matrix.
    Writing orientation priors.
    [done]
In [7]:
# inverse_operator
print(inverse_operator['eigen_fields'])
print(inverse_operator['eigen_leads_weighted'])
print(inverse_operator['fmri_prior'])
print(inverse_operator['orient_prior'])
print(inverse_operator['mri_head_t'])
print(inverse_operator['sing'].shape)
print(inverse_operator['nsource'])
print(inverse_operator['info'])
print(inverse_operator['src'])
print(inverse_operator['source_cov'])
print(inverse_operator['eigen_leads'])
{'row_names': [], 'ncol': 305L, 'data': array([[ -1.04083409e-17,  -1.66533454e-16,   2.77555756e-17, ...,
         -2.17850920e-02,  -1.10465729e-01,  -1.12777425e-01],
       [  2.77555756e-17,   1.66533454e-16,  -1.11022302e-16, ...,
         -2.45658141e-02,  -2.57169245e-02,   2.27786326e-01],
       [ -5.55111512e-17,  -2.22044605e-16,  -1.11022302e-16, ...,
          1.99798856e-01,  -2.37947476e-01,  -8.74727465e-02],
       ..., 
       [  2.23704026e-01,   9.62248200e-01,  -1.55031965e-01, ...,
         -4.34097203e-14,   4.14113188e-14,   6.34665931e-14],
       [  3.15244699e-01,  -2.21947143e-01,  -9.22691848e-01, ...,
          8.50992454e-14,  -5.93830540e-14,  -1.56527569e-13],
       [ -9.22267471e-01,   1.57536876e-01,  -3.52994114e-01, ...,
          6.57460197e-15,  -3.86149446e-15,  -1.07067133e-14]]), 'nrow': 305L, 'col_names': [u'MEG 0113', u'MEG 0112', u'MEG 0111', u'MEG 0122', u'MEG 0123', u'MEG 0121', u'MEG 0132', u'MEG 0133', u'MEG 0131', u'MEG 0143', u'MEG 0142', u'MEG 0141', u'MEG 0213', u'MEG 0212', u'MEG 0211', u'MEG 0222', u'MEG 0223', u'MEG 0221', u'MEG 0232', u'MEG 0233', u'MEG 0231', u'MEG 0243', u'MEG 0242', u'MEG 0241', u'MEG 0313', u'MEG 0312', u'MEG 0311', u'MEG 0322', u'MEG 0323', u'MEG 0321', u'MEG 0333', u'MEG 0332', u'MEG 0331', u'MEG 0343', u'MEG 0342', u'MEG 0341', u'MEG 0413', u'MEG 0412', u'MEG 0411', u'MEG 0422', u'MEG 0423', u'MEG 0421', u'MEG 0432', u'MEG 0433', u'MEG 0431', u'MEG 0443', u'MEG 0442', u'MEG 0441', u'MEG 0513', u'MEG 0512', u'MEG 0511', u'MEG 0523', u'MEG 0522', u'MEG 0521', u'MEG 0532', u'MEG 0533', u'MEG 0531', u'MEG 0542', u'MEG 0543', u'MEG 0541', u'MEG 0613', u'MEG 0612', u'MEG 0611', u'MEG 0622', u'MEG 0623', u'MEG 0621', u'MEG 0633', u'MEG 0632', u'MEG 0631', u'MEG 0642', u'MEG 0643', u'MEG 0641', u'MEG 0713', u'MEG 0712', u'MEG 0711', u'MEG 0723', u'MEG 0722', u'MEG 0721', u'MEG 0733', u'MEG 0732', u'MEG 0731', u'MEG 0743', u'MEG 0742', u'MEG 0741', u'MEG 0813', u'MEG 0812', u'MEG 0811', u'MEG 0822', u'MEG 0823', u'MEG 0821', u'MEG 0913', u'MEG 0912', u'MEG 0911', u'MEG 0923', u'MEG 0922', u'MEG 0921', u'MEG 0932', u'MEG 0933', u'MEG 0931', u'MEG 0942', u'MEG 0943', u'MEG 0941', u'MEG 1013', u'MEG 1012', u'MEG 1011', u'MEG 1023', u'MEG 1022', u'MEG 1021', u'MEG 1032', u'MEG 1033', u'MEG 1031', u'MEG 1043', u'MEG 1042', u'MEG 1041', u'MEG 1112', u'MEG 1113', u'MEG 1111', u'MEG 1123', u'MEG 1122', u'MEG 1121', u'MEG 1133', u'MEG 1132', u'MEG 1131', u'MEG 1142', u'MEG 1143', u'MEG 1141', u'MEG 1213', u'MEG 1212', u'MEG 1211', u'MEG 1223', u'MEG 1222', u'MEG 1221', u'MEG 1232', u'MEG 1233', u'MEG 1231', u'MEG 1243', u'MEG 1242', u'MEG 1241', u'MEG 1312', u'MEG 1313', u'MEG 1311', u'MEG 1323', u'MEG 1322', u'MEG 1321', u'MEG 1333', u'MEG 1332', u'MEG 1331', u'MEG 1342', u'MEG 1343', u'MEG 1341', u'MEG 1412', u'MEG 1413', u'MEG 1411', u'MEG 1423', u'MEG 1422', u'MEG 1421', u'MEG 1433', u'MEG 1432', u'MEG 1431', u'MEG 1442', u'MEG 1443', u'MEG 1441', u'MEG 1512', u'MEG 1513', u'MEG 1511', u'MEG 1522', u'MEG 1523', u'MEG 1521', u'MEG 1533', u'MEG 1532', u'MEG 1531', u'MEG 1543', u'MEG 1542', u'MEG 1541', u'MEG 1613', u'MEG 1612', u'MEG 1611', u'MEG 1622', u'MEG 1623', u'MEG 1621', u'MEG 1632', u'MEG 1633', u'MEG 1631', u'MEG 1643', u'MEG 1642', u'MEG 1641', u'MEG 1713', u'MEG 1712', u'MEG 1711', u'MEG 1722', u'MEG 1723', u'MEG 1721', u'MEG 1732', u'MEG 1733', u'MEG 1731', u'MEG 1743', u'MEG 1742', u'MEG 1741', u'MEG 1813', u'MEG 1812', u'MEG 1811', u'MEG 1822', u'MEG 1823', u'MEG 1821', u'MEG 1832', u'MEG 1833', u'MEG 1831', u'MEG 1843', u'MEG 1842', u'MEG 1841', u'MEG 1912', u'MEG 1913', u'MEG 1911', u'MEG 1923', u'MEG 1922', u'MEG 1921', u'MEG 1932', u'MEG 1933', u'MEG 1931', u'MEG 1943', u'MEG 1942', u'MEG 1941', u'MEG 2013', u'MEG 2012', u'MEG 2011', u'MEG 2023', u'MEG 2022', u'MEG 2021', u'MEG 2032', u'MEG 2033', u'MEG 2031', u'MEG 2042', u'MEG 2043', u'MEG 2041', u'MEG 2113', u'MEG 2112', u'MEG 2111', u'MEG 2122', u'MEG 2123', u'MEG 2121', u'MEG 2133', u'MEG 2132', u'MEG 2131', u'MEG 2143', u'MEG 2142', u'MEG 2141', u'MEG 2212', u'MEG 2213', u'MEG 2211', u'MEG 2223', u'MEG 2222', u'MEG 2221', u'MEG 2233', u'MEG 2232', u'MEG 2231', u'MEG 2242', u'MEG 2243', u'MEG 2241', u'MEG 2312', u'MEG 2313', u'MEG 2311', u'MEG 2323', u'MEG 2322', u'MEG 2321', u'MEG 2332', u'MEG 2333', u'MEG 2331', u'MEG 2343', u'MEG 2342', u'MEG 2341', u'MEG 2412', u'MEG 2413', u'MEG 2411', u'MEG 2423', u'MEG 2422', u'MEG 2421', u'MEG 2433', u'MEG 2432', u'MEG 2431', u'MEG 2442', u'MEG 2441', u'MEG 2512', u'MEG 2513', u'MEG 2511', u'MEG 2522', u'MEG 2523', u'MEG 2521', u'MEG 2533', u'MEG 2532', u'MEG 2531', u'MEG 2543', u'MEG 2542', u'MEG 2541', u'MEG 2612', u'MEG 2613', u'MEG 2611', u'MEG 2623', u'MEG 2622', u'MEG 2621', u'MEG 2633', u'MEG 2632', u'MEG 2631', u'MEG 2642', u'MEG 2643', u'MEG 2641']}
False
None
{'dim': 22494, 'kind': 6, 'bads': [], 'projs': [], 'diag': True, 'nfree': 1, 'eig': None, 'eigvec': None, 'data': array([ 0.2,  0.2,  1. , ...,  0.2,  0.2,  1. ]), 'names': []}
<Transform  |  MRI (surface RAS)->head>
[[ 0.99930954  0.00998479 -0.03578702 -0.00316745]
 [ 0.01275934  0.81240475  0.58295429  0.00685511]
 [ 0.0348942  -0.58300853  0.81171638  0.02888404]
 [ 0.          0.          0.          1.        ]]
(305L,)
7498
<Info | 12 non-empty fields
    bads : list | 0 items
    ch_names : list | MEG 0113, MEG 0112, MEG 0111, MEG 0122, MEG 0123, ...
    chs : list | 306 items (GRAD: 204, MAG: 102)
    command_line : unicode | 369 items
    custom_ref_applied : bool | False
    dev_head_t : 'mne.transforms.Transform | 3 items
    meas_file : unicode | 21 items
    mri_file : unicode | 29 items
    mri_head_t : 'mne.transforms.Transform | 3 items
    mri_id : dict | 4 items
    nchan : int | 306
    working_dir : unicode | 79 items
    meas_id : NoneType
>
<SourceSpaces: [<surface (lh), n_vertices=155407, n_used=3732, coordinate_frame=head>, <surface (rh), n_vertices=156866, n_used=3766, coordinate_frame=head>]>
{'dim': 22494, 'kind': 2, 'bads': [], 'projs': [], 'diag': True, 'nfree': 1, 'names': [], 'eig': None, 'data': array([  2.57891014e-19,   2.57891014e-19,   1.28945507e-18, ...,
         4.75741249e-19,   4.75741249e-19,   2.37870624e-18]), 'eigvec': None}
{'row_names': [], 'ncol': 305L, 'data': array([[ -2.36380081e-04,   1.57956037e-03,  -3.68674003e-03, ...,
          0.00000000e+00,   9.79178253e-01,   0.00000000e+00],
       [  1.29835858e-03,   1.06021254e-03,   1.01300425e-02, ...,
          5.32719348e-02,   2.33236412e-02,  -9.63053092e-01],
       [ -1.53428250e-03,  -5.71567382e-03,  -6.24716934e-03, ...,
          8.26868452e-01,   4.73989620e-02,  -9.43302318e-03],
       ..., 
       [ -1.15396694e-03,  -1.11017075e-03,   7.32175468e-03, ...,
          1.13411702e-04,  -1.01822729e-05,  -1.65034633e-05],
       [  2.78521200e-03,   4.24980440e-03,  -2.99328070e-03, ...,
          1.66715067e-05,  -2.80029332e-06,   6.43944371e-06],
       [  9.25545923e-03,   1.71409166e-02,   6.32775633e-03, ...,
         -1.02517033e-04,  -1.33913434e-04,  -4.63036493e-05]]), 'col_names': [], 'nrow': 22494L}

逆問題を解く

In [8]:
method='dSPM'
snr=3.
lambda2=1./snr**2
stc=apply_inverse(evoked,inverse_operator,lambda2,method=method,pick_ori=None)
# del fwd,inverse_operator,epochs # メモリ節約
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 55
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 3)
    Created the whitener using a full noise covariance matrix (3 small eigenvalues omitted)
    Computing noise-normalization factors (dSPM)...
[done]
Picked 305 channels from the data
Computing inverse...
(eigenleads need to be weighted)...
combining the current components...
(dSPM)...
[done]
In [9]:
plt.plot(1e3*stc.times,stc.data[::100,:].T)
plt.xlabel('time (ms)')
plt.ylabel('%s value' % method)
plt.show()