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

CTFのAEFサンプル

準備

In [2]:
# Authors: Mainak Jas <mainak.jas@telecom-paristech.fr>
#          Eric Larson <larson.eric.d@gmail.com>
#          Jaakko Leppakangas <jaeilepp@student.jyu.fi>
#
# License: BSD (3-clause)
In [3]:
import os.path as op
import pandas as pd
import numpy as np
import mne
mne.set_log_level('INFO')
from mne import combine_evoked
from mne.minimum_norm import apply_inverse
from mne.datasets.brainstorm import bst_auditory
from mne.io import read_raw_ctf
from mne.filter import notch_filter, filter_data
print(__doc__)
Automatically created module for IPython interactive environment

CTFデータの読み込み

In [4]:
use_precomputed = True # 予め計算したデータを利用する 計算しなおすときはFalse
In [5]:
data_path=bst_auditory.data_path() # データがないときには1.5GBダウンロードしに行きます。

CTF275 systemで2400Hzで計測し、600Hzにdownsampleするそうです。

In [6]:
subject='bst_auditory'
subjects_dir=data_path+'\\subjects'
raw_fname1=data_path+'\\MEG\\bst_auditory\\S01_AEF_20131218_01.ds'
raw_fname2=data_path+'\\MEG\\bst_auditory\\S01_AEF_20131218_02.ds'
erm_fname=data_path+'\\MEG\\bst_auditory\\S01_Noise_20131218_01.ds'
In [7]:
preload=not use_precomputed
raw=read_raw_ctf(raw_fname1,preload=preload)
n_times_run1=raw.n_times
mne.io.concatenate_raws([raw,read_raw_ctf(raw_fname2,preload=preload)])
raw_erm=read_raw_ctf(erm_fname,preload=preload)
ds directory : C:\Users\neuromag\mne_data\MNE-brainstorm-data\bst_auditory\MEG\bst_auditory\S01_AEF_20131218_01.ds
    res4 data read.
    hc data read.
    Separate EEG position data file read.
    Quaternion matching (desired vs. transformed):
       2.51   74.26    0.00 mm <->    2.51   74.26    0.00 mm (orig :  -56.69   50.20 -264.38 mm) diff =    0.000 mm
      -2.51  -74.26    0.00 mm <->   -2.51  -74.26    0.00 mm (orig :   50.89  -52.31 -265.88 mm) diff =    0.000 mm
     108.63    0.00    0.00 mm <->  108.63    0.00    0.00 mm (orig :   67.41   77.68 -239.53 mm) diff =    0.000 mm
    Coordinate transformations established.
    Reading digitizer points from [u'C:\\Users\\neuromag\\mne_data\\MNE-brainstorm-data\\bst_auditory\\MEG\\bst_auditory\\S01_AEF_20131218_01.ds\\S01_20131218_01.pos']...
    Polhemus data for 3 HPI coils added
    Device coordinate locations for 3 HPI coils added
    5 extra points added to Polhemus data.
    Measurement info composed.
Finding samples for C:\Users\neuromag\mne_data\MNE-brainstorm-data\bst_auditory\MEG\bst_auditory\S01_AEF_20131218_01.ds\S01_AEF_20131218_01.meg4: 
    System clock channel is available, checking which samples are valid.
    360 x 2400 = 864000 samples from 340 chs
Current compensation grade : 3
ds directory : C:\Users\neuromag\mne_data\MNE-brainstorm-data\bst_auditory\MEG\bst_auditory\S01_AEF_20131218_02.ds
    res4 data read.
    hc data read.
    Separate EEG position data file read.
    Quaternion matching (desired vs. transformed):
       2.64   74.60    0.00 mm <->    2.64   74.60    0.00 mm (orig :  -58.07   49.23 -263.11 mm) diff =    0.000 mm
      -2.64  -74.60    0.00 mm <->   -2.64  -74.60    0.00 mm (orig :   49.94  -53.82 -265.07 mm) diff =    0.000 mm
     108.24    0.00    0.00 mm <->  108.24   -0.00    0.00 mm (orig :   66.67   76.99 -243.39 mm) diff =    0.000 mm
    Coordinate transformations established.
    Reading digitizer points from [u'C:\\Users\\neuromag\\mne_data\\MNE-brainstorm-data\\bst_auditory\\MEG\\bst_auditory\\S01_AEF_20131218_02.ds\\S01_20131218_01.pos']...
    Polhemus data for 3 HPI coils added
    Device coordinate locations for 3 HPI coils added
    5 extra points added to Polhemus data.
    Measurement info composed.
Finding samples for C:\Users\neuromag\mne_data\MNE-brainstorm-data\bst_auditory\MEG\bst_auditory\S01_AEF_20131218_02.ds\S01_AEF_20131218_02.meg4: 
    System clock channel is available, checking which samples are valid.
    360 x 2400 = 864000 samples from 340 chs
Current compensation grade : 3
ds directory : C:\Users\neuromag\mne_data\MNE-brainstorm-data\bst_auditory\MEG\bst_auditory\S01_Noise_20131218_01.ds
    res4 data read.
    hc data read.
    Separate EEG position data file read.
    Quaternion matching (desired vs. transformed):
       0.00   80.00    0.00 mm <->    0.00   80.00    0.00 mm (orig :  -56.57   56.57 -270.00 mm) diff =    0.000 mm
       0.00  -80.00    0.00 mm <->    0.00  -80.00    0.00 mm (orig :   56.57  -56.57 -270.00 mm) diff =    0.000 mm
      80.00    0.00    0.00 mm <->   80.00   -0.00    0.00 mm (orig :   56.57   56.57 -270.00 mm) diff =    0.000 mm
    Coordinate transformations established.
    Polhemus data for 3 HPI coils added
    Device coordinate locations for 3 HPI coils added
    Measurement info composed.
Finding samples for C:\Users\neuromag\mne_data\MNE-brainstorm-data\bst_auditory\MEG\bst_auditory\S01_Noise_20131218_01.ds\S01_Noise_20131218_01.meg4: 
    System clock channel is available, checking which samples are valid.
    15 x 4800 = 72000 samples from 301 chs
Current compensation grade : 3
In [8]:
print(raw)
print(raw.info) 
# print(raw.info['ch_names'])
<RawCTF  |  S01_AEF_20131218_01.meg4, n_channels x n_times : 340 x 1728000 (720.0 sec), ~1.3 MB, data not loaded>
<Info | 21 non-empty fields
    bads : list | 0 items
    buffer_size_sec : numpy.float64 | 1.0
    ch_names : list | UDIO001, UPPT001, UTRG001, SCLK01-177, BG1-4408, ...
    chs : list | 340 items (REF_MEG: 26, EEG: 5, STIM: 3, MISC: 32, MAG: 274)
    comps : list | 5 items
    ctf_head_t : 'mne.transforms.Transform | 3 items
    custom_ref_applied : bool | False
    dev_ctf_t : 'mne.transforms.Transform | 3 items
    dev_head_t : 'mne.transforms.Transform | 3 items
    dig : list | 263 items
    events : list | 0 items
    experimenter : unicode | 3 items
    highpass : float | 0.0 Hz
    hpi_meas : list | 0 items
    hpi_results : list | 1 items
    lowpass : float | 1200.0 Hz
    meas_id : dict | 4 items
    nchan : int | 340
    projs : list | 0 items
    sfreq : float | 2400.0 Hz
    subject_info : dict | 1 items
    acq_pars : NoneType
    acq_stim : NoneType
    description : NoneType
    file_id : NoneType
    hpi_subsystem : NoneType
    kit_system_id : NoneType
    line_freq : NoneType
    meas_date : NoneType
    proj_id : NoneType
    proj_name : NoneType
    xplotter_layout : NoneType
>
In [9]:
raw.set_channel_types({'HEOG':'eog','VEOG':'eog','ECG':'ecg'}) # 水平と垂直方向の眼電図HEOGとVEOGをeogとする
In [10]:
if not use_precomputed:
    raw.pick_types(meg=True,eeg=False,stim=True,misc=True,eog=True,ecg=True)
In [11]:
annotations_df=pd.DataFrame()
offset=n_times_run1
print(annotations_df,offset)
(Empty DataFrame
Columns: []
Index: [], 864000)
In [12]:
for idx in [1,2]:
    csv_fname=data_path+'\\MEG\\bst_auditory\\events_bad_0%s.csv' % idx
    print(csv_fname)
    df=pd.read_csv(csv_fname,header=None,names=['onset','duration','id','label']) # pdはpandas
    print('Events from run {0}:'.format(idx))
    print(df)
    df['onset']+=offset*(idx-1)
    annotations_df=pd.concat([annotations_df,df],axis=0)
saccades_events=df[df['label']=='saccade'].values[:,:3].astype(int)
# サンプルから秒へ
onsets=annotations_df['onset'].values/raw.info['sfreq']
durations=annotations_df['duration'].values/raw.info['sfreq']
descriptions=annotations_df['label'].values

annotations=mne.Annotations(onsets,durations,descriptions)
raw.annotations=annotations
del onsets,durations,descriptions
C:\Users\neuromag\mne_data\MNE-brainstorm-data\bst_auditory\MEG\bst_auditory\events_bad_01.csv
Events from run 1:
     onset  duration  id label
0     7625      2776   1   BAD
1   142459       892   1   BAD
2   216954       460   1   BAD
3   345135      5816   1   BAD
4   357687      1053   1   BAD
5   409101      3736   1   BAD
6   461110       179   1   BAD
7   479866       426   1   BAD
8   764914     11500   1   BAD
9   798174      6589   1   BAD
10  846880      5383   1   BAD
11  858863      5136   1   BAD
C:\Users\neuromag\mne_data\MNE-brainstorm-data\bst_auditory\MEG\bst_auditory\events_bad_02.csv
Events from run 2:
     onset  duration  id    label
0        9      5583   1      BAD
1     9256      3114   1      BAD
2    14287      3456   1      BAD
3   116432       228   1      BAD
4   134489      1329   1      BAD
5   464527      4727   1      BAD
6   494136      4519   1      BAD
7   749288       189   1      BAD
8   788623      7937   1      BAD
9    21179         0   1  saccade
10   72993         0   1  saccade
11  134527         0   1  saccade
12  196555         0   1  saccade
13  249894         0   1  saccade
14  343357         0   1  saccade
15  400771         0   1  saccade
16  450256         0   1  saccade
17  593101         0   1  saccade
18  733942         0   1  saccade
19  765939         0   1  saccade
20  789476         0   1  saccade
21  792852         0   1  saccade
22  833208         0   1  saccade
23  859869         0   1  saccade
24  862888         0   1  saccade
In [13]:
saccade_epochs=mne.Epochs(raw,saccades_events,1,0.,0.5,preload=True,reject_by_annotation=False)
projs_saccade=mne.compute_proj_epochs(saccade_epochs,n_mag=1,n_eeg=0,desc_prefix='saccade')
if use_precomputed:
    proj_fname=data_path+'\\MEG\\bst_auditory\\bst_auditory-eog-proj.fif'
    projs_eog=mne.read_proj(proj_fname)[0]
else:
    projs_eog,_=mne.preprocessing.compute_proj_eog(raw.load.data(),n_mag=1,n_eeg=0)
raw.add_proj(projs_saccade)
raw.add_proj(projs_eog)
del saccade_epochs,saccades_events,projs_eog,projs_saccade
16 matching events found
0 projection items activated
Loading data for 16 events and 1201 original time points ...
1 bad epochs dropped
No gradiometers found. Forcing n_grad to 0
Adding projection: axial-saccade-PCA-01
    Read a total of 1 projection items:
        EOG-axial-998--0.200-0.200-PCA-01 (1 x 274)  idle
1 projection items deactivated
1 projection items deactivated
In [14]:
raw.plot(block=True);

フィルタ処理

In [15]:
if not use_precomputed:
    meg_picks=mne.pick_types(rwa.info,meg=True,eeg=False)
    raw.plot_psd(tmax=np.inf,picks=meg_picks)
    notches=np.arange(60,181,60)
    raw.notch_filter(notches)
    raw.plot_psd(tmax=np.inf,picks=meg_picks)
In [16]:
if not use_precomputed:
    raw.filter(None,100.,h_trans_bandwidth=0.5,filter_length='10s',phase='zero-double')

試行抽出

In [17]:
tmin,tmax=-0.1,0.5
event_id=dict(standard=1,deviant=2)
reject=dict(mag=4e-12,eog=250e-6)
events=mne.find_events(raw,stim_channel='UPPT001')
480 events found
Events id: [1 2]
In [18]:
raw.info['bads']=['MLO52-4408','MRT51-4408','MLO42-4408','MLO43-4408'] #0ゼロでなくOオー
picks=mne.pick_types(raw.info,meg=True,eeg=False,stim=False,eog=True,exclude='bads')
epochs=mne.Epochs(raw,events,event_id,tmin,tmax,picks=picks,baseline=(None,0),reject=reject,preload=False,proj=True)
480 matching events found
Created an SSP operator (subspace dimension = 2)
2 projection items activated
In [19]:
epochs.drop_bad()
epochs_standard=mne.concatenate_epochs([epochs['standard'][range(40)],epochs['standard'][182:222]])
epochs_standard.load_data()
epochs_standard.resample(600,npad='auto')
epochs_deviant=epochs['deviant'].load_data()
epochs_deviant.resample(600,npad='auto')
del epochs,picks
Loading data for 480 events and 1441 original time points ...
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on MAG : [u'MLP52-4408']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'HEOG']
    Rejecting  epoch based on EOG : [u'HEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on MAG : [u'MLP52-4408']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'HEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'HEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'HEOG']
    Rejecting  epoch based on EOG : [u'HEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'HEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
    Rejecting  epoch based on EOG : [u'VEOG']
39 bad epochs dropped
Loading data for 40 events and 1441 original time points ...
Loading data for 40 events and 1441 original time points ...
80 matching events found
Created an SSP operator (subspace dimension = 2)
0 bad epochs dropped
Loading data for 75 events and 1441 original time points ...
In [20]:
evoked_std=epochs_standard.average()
evoked_dev=epochs_deviant.average()
del epochs_standard,epochs_deviant
In [21]:
if use_precomputed:
    sfreq=evoked_std.info['sfreq']
    notches=[60,120,180]
    for evoked in (evoked_std,evoked_dev):
        evoked.data[:]=notch_filter(evoked.data,sfreq,notches)
        evoked.data[:]=filter_data(evoked.data,sfreq,l_freq=None,h_freq=100.)
Setting up band-stop filter
Filter length of 7920 samples (13.200 sec) selected
filter_length (7921) is longer than the signal (360), distortion is likely. Reduce filter length or filter a longer signal.
Setting up low-pass filter at 1e+02 Hz
h_trans_bandwidth chosen to be 25.0 Hz
Filter length of 158 samples (0.263 sec) selected
Setting up band-stop filter
Filter length of 7920 samples (13.200 sec) selected
<ipython-input-21-1b2e0a3fcaf0>:5: RuntimeWarning: filter_length (7921) is longer than the signal (360), distortion is likely. Reduce filter length or filter a longer signal.
  evoked.data[:]=notch_filter(evoked.data,sfreq,notches)
filter_length (7921) is longer than the signal (360), distortion is likely. Reduce filter length or filter a longer signal.
Setting up low-pass filter at 1e+02 Hz
h_trans_bandwidth chosen to be 25.0 Hz
Filter length of 158 samples (0.263 sec) selected
<ipython-input-21-1b2e0a3fcaf0>:5: RuntimeWarning: filter_length (7921) is longer than the signal (360), distortion is likely. Reduce filter length or filter a longer signal.
  evoked.data[:]=notch_filter(evoked.data,sfreq,notches)

波形・トポ描画

In [22]:
evoked_std.plot(window_title='Standard',gfp=True);
evoked_dev.plot(window_title='Deviant',gfp=True);
In [23]:
times=np.arange(0.05,0.301,0.025)
evoked_std.plot_topomap(times=times,title='Standard');
evoked_dev.plot_topomap(times=times,title='Deviant');