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

SSP

準備です。

In [2]:
import numpy as np
import mne
from mne.datasets import sample
from mne.preprocessing import compute_proj_ecg,compute_proj_eog
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,preload=True)
raw.set_eeg_reference()
raw.pick_types(meg=True,ecg=True,eog=True,stim=True)
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
Reading 0 ... 41699  =      0.000 ...   277.709 secs...
An average reference projection was already added. The data has been left untouched.
<ipython-input-3-60410e94f0a2>:2: RuntimeWarning: An average reference projection was already added. The data has been left untouched.
  raw.set_eeg_reference()
Out[3]:
<Raw  |  sample_audvis_filt-0-40_raw.fif, n_channels x n_times : 315 x 41700 (277.7 sec), ~103.8 MB, data loaded>

心電図から固有ベクトル抽出

心電図で加算平均して、主成分分析して、平面型グラジオメータ1個、マグネトメータ1個の固有値ベクトルを求めます。

In [4]:
projs,events=compute_proj_ecg(raw,n_grad=1,n_mag=1,average=True)
projs
Including 4 SSP projectors from raw file
Running ECG SSP computation
Reconstructing ECG signal from Magnetometers
Setting up band-pass filter from 5 - 35 Hz
Number of ECG events detected : 285 (average pulse 61 / min.)
Computing projector
Setting up band-pass filter from 1 - 35 Hz
285 matching events found
Created an SSP operator (subspace dimension = 3)
4 projection items activated
Loading data for 285 events and 91 original time points ...
    Rejecting  epoch based on MAG : [u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1411', u'MEG 1421']
    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 MAG : [u'MEG 1711']
    Rejecting  epoch based on MAG : [u'MEG 1411', u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1411']
    Rejecting  epoch based on MAG : [u'MEG 1421']
    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 1411']
    Rejecting  epoch based on MAG : [u'MEG 1441']
    Rejecting  epoch based on EOG : [u'EOG 061']
15 bad epochs dropped
No EEG channels found. Forcing n_eeg to 0
Adding projection: planar--0.200-0.400-PCA-01
Adding projection: axial--0.200-0.400-PCA-01
Done.
Out[4]:
[<Projection  |  PCA-v1, active : False, n_channels : 102>,
 <Projection  |  PCA-v2, active : False, n_channels : 102>,
 <Projection  |  PCA-v3, active : False, n_channels : 102>,
 <Projection  |  Average EEG reference, active : False, n_channels : 60>,
 <Projection  |  ECG-planar--0.200-0.400-PCA-01, active : False, n_channels : 203>,
 <Projection  |  ECG-axial--0.200-0.400-PCA-01, active : False, n_channels : 102>]

新たに求めた固有ベクトルがECG-planar--0.200・・・とECG-axial--0.200・・・です。

In [5]:
events.shape
Out[5]:
(285L, 3L)
In [6]:
ecg_projs=projs[-2:] # MATLABだとprojs[end-1:end]
print(ecg_projs[0])
print(ecg_projs[1])
mne.viz.plot_projs_topomap(ecg_projs);
<Projection  |  ECG-planar--0.200-0.400-PCA-01, active : False, n_channels : 203>
<Projection  |  ECG-axial--0.200-0.400-PCA-01, active : False, n_channels : 102>

眼電図から固有ベクトル抽出

In [7]:
projs,events=compute_proj_eog(raw,n_grad=1,n_mag=1,average=True)
projs
Including 4 SSP projectors from raw file
Running EOG SSP computation
EOG channel index for this subject is: [314]
Filtering the data to remove DC offset to help distinguish blinks from saccades
Setting up band-pass filter from 2 - 45 Hz
Setting up band-pass filter from 1 - 10 Hz
Now detecting blinks and generating corresponding events
Number of EOG events detected : 46
Computing projector
Setting up band-pass filter from 1 - 35 Hz
46 matching events found
Created an SSP operator (subspace dimension = 3)
4 projection items activated
Loading data for 46 events and 61 original time points ...
    Rejecting  epoch based on MAG : [u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1411', u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1411', u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1411']
    Rejecting  epoch based on MAG : [u'MEG 1411', u'MEG 1421']
5 bad epochs dropped
No EEG channels found. Forcing n_eeg to 0
Adding projection: planar--0.200-0.200-PCA-01
Adding projection: axial--0.200-0.200-PCA-01
Done.
Out[7]:
[<Projection  |  PCA-v1, active : False, n_channels : 102>,
 <Projection  |  PCA-v2, active : False, n_channels : 102>,
 <Projection  |  PCA-v3, active : False, n_channels : 102>,
 <Projection  |  Average EEG reference, active : False, n_channels : 60>,
 <Projection  |  EOG-planar--0.200-0.200-PCA-01, active : False, n_channels : 203>,
 <Projection  |  EOG-axial--0.200-0.200-PCA-01, active : False, n_channels : 102>]
In [8]:
eog_projs=projs[-2:] # MATLABだとprojs[end-1:end]
mne.viz.plot_projs_topomap(eog_projs);

SSPの適応

projsにSSPを追加します。

In [9]:
raw.info["projs"]+=eog_projs+ecg_projs
raw.info["projs"]
Out[9]:
[<Projection  |  PCA-v1, active : False, n_channels : 102>,
 <Projection  |  PCA-v2, active : False, n_channels : 102>,
 <Projection  |  PCA-v3, active : False, n_channels : 102>,
 <Projection  |  Average EEG reference, active : False, n_channels : 60>,
 <Projection  |  EOG-planar--0.200-0.200-PCA-01, active : False, n_channels : 203>,
 <Projection  |  EOG-axial--0.200-0.200-PCA-01, active : False, n_channels : 102>,
 <Projection  |  ECG-planar--0.200-0.400-PCA-01, active : False, n_channels : 203>,
 <Projection  |  ECG-axial--0.200-0.400-PCA-01, active : False, n_channels : 102>]
In [10]:
events=mne.find_events(raw,stim_channel='STI 014')
reject=dict(grad=4000e-13,mag=4e-12,eog=150e-6)

event_id={'auditory/left':1}
epochs_no_proj=mne.Epochs(raw,events,event_id,tmin=-0.2,tmax=0.5,proj=False,baseline=(None,0),reject=reject)
epochs_no_proj.average().plot(spatial_colors=True);
319 events found
Events id: [ 1  2  3  4  5 32]
72 matching events found
Created an SSP operator (subspace dimension = 7)
    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']
C:\Users\akira\Anaconda2\lib\site-packages\matplotlib\figure.py:1744: UserWarning: This figure includes Axes that are not compatible with tight_layout, so its results might be incorrect.
  warnings.warn("This figure includes Axes that are not "
In [11]:
epochs_proj=mne.Epochs(raw,events,event_id,tmin=-0.2,tmax=0.5,proj=True,baseline=(None,0),reject=reject)
epochs_proj.average().plot(spatial_colors=True);
72 matching events found
Created an SSP operator (subspace dimension = 7)
8 projection items activated
    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']

delayed SSPモード

In [12]:
epochs_delayed=mne.Epochs(raw,events,event_id,tmin=-0.2,tmax=0.5,proj='delayed',baseline=(None,0),reject=reject)
evoked=epochs_delayed.average();
evoked.plot(spatial_colors=True);
72 matching events found
Entering delayed SSP mode.
Created an SSP operator (subspace dimension = 7)
    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']