内容

はじめに... 1

Spyderで使ってみる... 1

脳磁図データの読み込み... 1

データの切り出し... 1

波形表示... 1

FIFファイルで保存... 1

エポックの切り出し... 1

MatlabのMATファイルへ書き出し... 1

加算平均波形の表示... 1

その他もろもろ... 1

各エポックの最大値を求め方... 1

加算平均されたFIFファイルの読み方... 1

差分波形の作成... 1

加算波形の加算... 1

エポック数を同じにした差分波形 trial-count equalization. 1

時間周波数解析... 1

空間フィルタの作成... 1

 

 

はじめに

MNEtutorialを使ってみました。

TutorialsIntroduction to MNE and PythonBasic MEG and EEG processingというところを試してみました。

これを試してみました。

MNE Pytonでできることとできないことが列挙してあります。意訳してみます。

MNE Pythonでできること

l  生データの可視化 C言語で書かれたmne_browse_rawの機能拡張

l  Epochの切り出し epochの切り出し、基線などの調整

l  加算平均

l  主成分分析(SSP) 心磁図や眼磁図の影響除去

l  独立成分分析 いろんな雑音除去

l  Maxwelフィルタ 環境雑音除去

l  境界要素モデル作成 一層・三層の境界要素モデルの作成

l  順問題を解く 境界要素などの導体モデルの作成

l  線形逆問題を解く dSPM sLORETA MNE LCMV DICSなど

l  疎行列の逆問題を解く L1/L2ノルムMxNE Gamma Map 時間周波数MxNE

l  連結性推定 センサ単位 電流源単位

l  センサ信号や電流源信号の可視化

l  時間周波数解析 Morletwavelet変換

l  皮質・皮質下で電流源推定

l  電流双極子推定

l  復号 脳波・脳磁図トポグラフィーの多変数パターン解析

l  差の強調 条件間 センサ間 被検者間の

l  非正規分布統計処理 時間 空間 周波数

l  プログラミング バッチ処理 高速分散処理

l  各種データのFIFCTF BTI/4D KIT 各種脳波データのFIF

MNE-Pythonでできないこと

l  脳の切り出し FreeSurferを使いましょう

l  記録時の動き補正 ElektaMaxfilterTMを使いましょう

これだけできればとりあえず十分かなとは思いますが、脳の切り出しがFreeSurfer依存なのは・・・・ここがbottle neckですかね。臨床では使いにくいです。

 

Spyderで使ってみる

脳磁図データの読み込み

Spyderを起動します。IPythonコンソールではなくPython consoleを使うことにします。

#はコメント文です。

>>> >>> import mne # MNE-pythonの組み込み

>>> mne.get_config_path() # mneの場所表示

'C:\\Users\\akira\\AppData\\Roaming\\.mne\\mne-python.json'

LinuxOSXではフォルダとフォルダの仕切りは/ですがWindowsだと\\(バックスラッシュ・バックスラッシュ)になります。

>>> from mne.datasets import sample # mne.datasetからsampleの読み込み

次に以下の文を実行すると、2GBのデータをダウンロードしに行きます。相当時間がかかります。

>>> data_path=sample.data_path() #2GBのデータのダウンロードに相当時間がかかります。

>>> raw_fname=data_path+’\\Users\\akira\\mne_data\\MNE-sample-data\\MEG\\sample\\sample_audvis_filt-0-40_raw.fif'

一旦データを読み込んだなら、2回目以降は以下のコマンドでデータを読みに行く時間が省けます。

>>> raw_fname='C:\\Users\\akira\\mne_data\\MNE-sample-data\\MEG\\sample\\sample_audvis_filt-0-40_raw.fif'

>>> print(raw)

C:\Users\akira\mne_data\MNE-sample-data\\MEG\sample\sample_audvis_filt-0-40_raw.fif

では実際にファイルの情報を読みにいきます。

>>> raw=mne.io.read_raw_fif(raw_fname,add_eeg_ref=False)

読み込んだ内容の表示です。

>>> print(raw)

<Raw  |  sample_audvis_filt-0-40_raw.fif, n_channels x n_times : 376 x 41700 (277.7 sec), ~3.6 MB, data not loaded>

376 ch×41700サンプル(277.7) 3.6MBと表示されます。

raw.infoにはもっと詳しい情報があります。

>>> print(raw.info)

<Info | 20 non-empty fields

    bads : 'list | MEG 2443, EEG 053

    buffer_size_sec : 'numpy.float64 | 13.3196808772

    ch_names : 'list | MEG 0113, MEG 0112, MEG 0111, MEG 0122, MEG 0123, ...

    chs : 'list | 376 items (EOG: 1, EEG: 60, STIM: 9, MAG: 102, GRAD: 204)

    comps : 'list | 0 items

    custom_ref_applied : 'bool | False

    dev_head_t : 'mne.transforms.Transform | 3 items

    dig : 'list | 146 items

    events : 'list | 0 items

    file_id : 'dict | 4 items

    filename : 'str | C:\Users\a.../sample_audvis_filt-0-40_raw.fif

    highpass : 'float | 0.10000000149011612 Hz

    hpi_meas : 'list | 1 items

    hpi_results : 'list | 1 items

    lowpass : 'float | 40.0 Hz

    meas_date : 'numpy.ndarray | 2002-12-04 04:01:10

    meas_id : 'dict | 4 items

    nchan : 'int | 376

    projs : 'list | PCA-v1: off, PCA-v2: off, PCA-v3: off, ...

    sfreq : 'float | 150.15374755859375 Hz

    acq_pars : 'NoneType

    acq_stim : 'NoneType

    ctf_head_t : 'NoneType

    description : 'NoneType

    dev_ctf_t : 'NoneType

    experimenter : 'NoneType

    hpi_subsystem : 'NoneType

    kit_system_id : 'NoneType

    line_freq : 'NoneType

    proj_id : 'NoneType

    proj_name : 'NoneType

    subject_info : 'NoneType

    xplotter_layout : 'NoneType

チャンネル名はraw.ch.namesに格納されています。

>>> print(raw.ch_names)

['MEG 0113', 'MEG 0112', 'MEG 0111', 'MEG 0122', 'MEG 0123', 'MEG 0121', 'MEG 0132', 'MEG 0133', 'MEG 0131', 'MEG 0143', 'MEG 0142', 'MEG 0141', 'MEG 0213', 'MEG 0212', 'MEG 0211', 'MEG 0222', 'MEG 0223', 'MEG 0221', 'MEG 0232', 'MEG 0233', 'MEG 0231', 'MEG 0243', 'MEG 0242', 'MEG 0241', 'MEG 0313', 'MEG 0312', 'MEG 0311', 'MEG 0322', 'MEG 0323', 'MEG 0321', 'MEG 0333', 'MEG 0332', 'MEG 0331', 'MEG 0343', 'MEG 0342', 'MEG 0341', 'MEG 0413', 'MEG 0412', 'MEG 0411', 'MEG 0422', 'MEG 0423', 'MEG 0421', 'MEG 0432', 'MEG 0433', 'MEG 0431', 'MEG 0443', 'MEG 0442', 'MEG 0441', 'MEG 0513', 'MEG 0512', 'MEG 0511', 'MEG 0523', 'MEG 0522', 'MEG 0521', 'MEG 0532', 'MEG 0533', 'MEG 0531', 'MEG 0542', 'MEG 0543', 'MEG 0541', 'MEG 0613', 'MEG 0612', 'MEG 0611', 'MEG 0622', 'MEG 0623', 'MEG 0621', 'MEG 0633', 'MEG 0632', 'MEG 0631', 'MEG 0642', 'MEG 0643', 'MEG 0641', 'MEG 0713', 'MEG 0712', 'MEG 0711', 'MEG 0723', 'MEG 0722', 'MEG 0721', 'MEG 0733', 'MEG 0732', 'MEG 0731', 'MEG 0743', 'MEG 0742', 'MEG 0741', 'MEG 0813', 'MEG 0812', 'MEG 0811', 'MEG 0822', 'MEG 0823', 'MEG 0821', 'MEG 0913', 'MEG 0912', 'MEG 0911', 'MEG 0923', 'MEG 0922', 'MEG 0921', 'MEG 0932', 'MEG 0933', 'MEG 0931', 'MEG 0942', 'MEG 0943', 'MEG 0941', 'MEG 1013', 'MEG 1012', 'MEG 1011', 'MEG 1023', 'MEG 1022', 'MEG 1021', 'MEG 1032', 'MEG 1033', 'MEG 1031', 'MEG 1043', 'MEG 1042', 'MEG 1041', 'MEG 1112', 'MEG 1113', 'MEG 1111', 'MEG 1123', 'MEG 1122', 'MEG 1121', 'MEG 1133', 'MEG 1132', 'MEG 1131', 'MEG 1142', 'MEG 1143', 'MEG 1141', 'MEG 1213', 'MEG 1212', 'MEG 1211', 'MEG 1223', 'MEG 1222', 'MEG 1221', 'MEG 1232', 'MEG 1233', 'MEG 1231', 'MEG 1243', 'MEG 1242', 'MEG 1241', 'MEG 1312', 'MEG 1313', 'MEG 1311', 'MEG 1323', 'MEG 1322', 'MEG 1321', 'MEG 1333', 'MEG 1332', 'MEG 1331', 'MEG 1342', 'MEG 1343', 'MEG 1341', 'MEG 1412', 'MEG 1413', 'MEG 1411', 'MEG 1423', 'MEG 1422', 'MEG 1421', 'MEG 1433', 'MEG 1432', 'MEG 1431', 'MEG 1442', 'MEG 1443', 'MEG 1441', 'MEG 1512', 'MEG 1513', 'MEG 1511', 'MEG 1522', 'MEG 1523', 'MEG 1521', 'MEG 1533', 'MEG 1532', 'MEG 1531', 'MEG 1543', 'MEG 1542', 'MEG 1541', 'MEG 1613', 'MEG 1612', 'MEG 1611', 'MEG 1622', 'MEG 1623', 'MEG 1621', 'MEG 1632', 'MEG 1633', 'MEG 1631', 'MEG 1643', 'MEG 1642', 'MEG 1641', 'MEG 1713', 'MEG 1712', 'MEG 1711', 'MEG 1722', 'MEG 1723', 'MEG 1721', 'MEG 1732', 'MEG 1733', 'MEG 1731', 'MEG 1743', 'MEG 1742', 'MEG 1741', 'MEG 1813', 'MEG 1812', 'MEG 1811', 'MEG 1822', 'MEG 1823', 'MEG 1821', 'MEG 1832', 'MEG 1833', 'MEG 1831', 'MEG 1843', 'MEG 1842', 'MEG 1841', 'MEG 1912', 'MEG 1913', 'MEG 1911', 'MEG 1923', 'MEG 1922', 'MEG 1921', 'MEG 1932', 'MEG 1933', 'MEG 1931', 'MEG 1943', 'MEG 1942', 'MEG 1941', 'MEG 2013', 'MEG 2012', 'MEG 2011', 'MEG 2023', 'MEG 2022', 'MEG 2021', 'MEG 2032', 'MEG 2033', 'MEG 2031', 'MEG 2042', 'MEG 2043', 'MEG 2041', 'MEG 2113', 'MEG 2112', 'MEG 2111', 'MEG 2122', 'MEG 2123', 'MEG 2121', 'MEG 2133', 'MEG 2132', 'MEG 2131', 'MEG 2143', 'MEG 2142', 'MEG 2141', 'MEG 2212', 'MEG 2213', 'MEG 2211', 'MEG 2223', 'MEG 2222', 'MEG 2221', 'MEG 2233', 'MEG 2232', 'MEG 2231', 'MEG 2242', 'MEG 2243', 'MEG 2241', 'MEG 2312', 'MEG 2313', 'MEG 2311', 'MEG 2323', 'MEG 2322', 'MEG 2321', 'MEG 2332', 'MEG 2333', 'MEG 2331', 'MEG 2343', 'MEG 2342', 'MEG 2341', 'MEG 2412', 'MEG 2413', 'MEG 2411', 'MEG 2423', 'MEG 2422', 'MEG 2421', 'MEG 2433', 'MEG 2432', 'MEG 2431', 'MEG 2442', 'MEG 2443', 'MEG 2441', 'MEG 2512', 'MEG 2513', 'MEG 2511', 'MEG 2522', 'MEG 2523', 'MEG 2521', 'MEG 2533', 'MEG 2532', 'MEG 2531', 'MEG 2543', 'MEG 2542', 'MEG 2541', 'MEG 2612', 'MEG 2613', 'MEG 2611', 'MEG 2623', 'MEG 2622', 'MEG 2621', 'MEG 2633', 'MEG 2632', 'MEG 2631', 'MEG 2642', 'MEG 2643', 'MEG 2641', 'STI 001', 'STI 002', 'STI 003', 'STI 004', 'STI 005', 'STI 006', 'STI 014', 'STI 015', 'STI 016', 'EEG 001', 'EEG 002', 'EEG 003', 'EEG 004', 'EEG 005', 'EEG 006', 'EEG 007', 'EEG 008', 'EEG 009', 'EEG 010', 'EEG 011', 'EEG 012', 'EEG 013', 'EEG 014', 'EEG 015', 'EEG 016', 'EEG 017', 'EEG 018', 'EEG 019', 'EEG 020', 'EEG 021', 'EEG 022', 'EEG 023', 'EEG 024', 'EEG 025', 'EEG 026', 'EEG 027', 'EEG 028', 'EEG 029', 'EEG >> >030', 'EEG 031', 'EEG 032', 'EEG 033', 'EEG 034', 'EEG 035', 'EEG 036', 'EEG 037', 'EEG 038', 'EEG 039', 'EEG 040', 'EEG 041', 'EEG 042', 'EEG 043', 'EEG 044', 'EEG 045', 'EEG 046', 'EEG 047', 'EEG 048', 'EEG 049', 'EEG 050', 'EEG 051', 'EEG 052', 'EEG 053', 'EEG 054', 'EEG 055', 'EEG 056', 'EEG 057', 'EEG 058', 'EEG 059', 'EEG 060', 'EOG 061']

 

データの切り出し

100115秒のサンプル番号を取得します。

>>> start,stop=raw.time_as_index([100,115]) #100115

>>> start

15015

>>> stop

17267

データを切り出し、大きさを調べてみます。

>>> data,times=raw[:,start:stop]

>>> print(data.shape)

(376, 2252)

>>> print(times.shape)

(2252,)

データ数は376×2252、時系列は2252サンプルと表示されます。

チャンネル数を制限する時は以下のようにします。

>>> data,times=raw[2:20:3,start:stop]

>>> print(data.shape)

(6, 2252)

Matlabの書き方だと2:3:20です。実際に選択されたチャンネルは以下になります。

>>> print(raw.ch_names[2:20:3])

['MEG 0111', 'MEG 0121', 'MEG 0131', 'MEG 0141', 'MEG 0211', 'MEG 0221']

 

波形表示

以下のコマンドでウインドウが開いて波形が表示されます。尚IPythonコンソールだとウインドウで表示された時の絵が出るだけです。

>>> raw.plot()

<matplotlib.figure.Figure object at 0x0000014910019860>

赤のカーソルのあるバーをクリックすると選択された時間・チャンネルが変わります。カーソルバーのドラッグには対応していません。縦バーの色はチャンネルを表し、青色が平面型グラジオメータ、藍色がマグネトメータ、黒が脳波です。

Helpを押すと使い方が表示されます。

スケールは変更できますが、チャンネル別のスケールの変更はできそうもありません。

Projを押すとSSPなどが表示されます。

 

FIFファイルで保存

チャンネル番号を選択し、0150秒を脳波チャンネルなしで保存します。

>>> picks=mne.pick_types(raw.info,meg=True,eeg=False,stim=True,exclude='bads')

>>> picks

array([  0,   1,   2, ..., 312, 313, 314])

>>> raw.save('sample_audvis_meg_raw.fif',tmin=0,tmax=150,picks=picks,overwrite=True)

>>> mne.get_config_path()

'C:\\Users\\akira\\AppData\\Roaming\\.mne\\mne-python.json'

Matlabでは番号は1,2,3・・・となりますが、Pythonでは0番から始まり、0,1,2,・・・となります。

パスが指定されてないので保存先はmne.get_config_path()で表示されるフォルダですが、ない!

C:\\Users\akiraに保存されてました。

 

エポックの切り出し

トリガーチャンネルを指定します。

>>> events=mne.find_events(raw,stim_channel='STI 014')

>>> events.shape

(319, 3)

>>> print(events[:6])

[[6994    0    2]

 [7086    0    3]

 [7192    0    1]

 [7304    0    4]

 [7413    0    2]

 [7506    0    3]]

MATLABだとevent(1:6,:)です。319×3の配列です。1列目はサンプル(時間)、3列目はトリガー信号の値です。

初期設定でSTI101をトリガー信号にする場合は

>>> mne.set_config('MNE_STIM_CHANNEL','STI101',set_env=True)

となります。初期設定があればmne.find_eventsstim_channelを省略できます。

>>> mne.set_config('MNE_STIM_CHANNEL','STI 014',set_env=True)

>>> events=mne.find_events(raw)

>>> events

array([[ 6994,     0,     2],

       [ 7086,     0,     3],

       [ 7192,     0,     1],

       ...,

       [38520,     0,     1],

       [38621,     0,     4],

       [42168,     0,    32]], dtype=int64)

エポック条件を定義します。

>>> event_id=dict(aud_l=1, aud_r=2) #左が1、右が2

>>> tmin=-0.2 #開始時間

>>> tmax=0.5 #時間幅

>>> raw.info['bads']

['MEG 2443', 'EEG 053']

>>> raw.info['bads']+=['MEG 2443', 'EEG 053'] # MNEtutorialbad channelを追加するようにしてあったので追加しました

>>> raw.info['bads']

['MEG 2443', 'EEG 053', 'MEG 2443', 'EEG 053'] # ダブってます。

>>> picks=mne.pick_types(raw.info,meg=True,eeg=True,eog=True,stim=False,exclude='bads')

>>> baseline=(None,0)  # 基線の設定 なし 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,preload=False,reject=reject,add_eeg_ref=False)

>>> epochs # tutorialだとprint(epochs)

<Epochs  |  n_events : 145 (good & bad), tmin : -0.199795213158 (s), tmax : 0.499488032896 (s), baseline : (None, 0), ~3.6 MB, data not loaded,

 'aud_l': 72, 'aud_r': 73

左だけのデータを見てみます。

>>> epochs_data=epochs['aud_l'].get_data()

>>> print(epochs_data.shape) # 前項が計算が終わるまで表示が反映されません。

(55, 365, 106)

epochs_dataは三次元配列で55エポック、365チャンネル、106サンプル時間です。

 

MatlabMATファイルへ書き出し

scipyを使います。

>>> from scipy import io

>>> io.savemat('epochs_data.mat',dict(epochs_data=epochs_data),oned_as='row') #ちょっと時間がかかります。

これもC:\\Users\akiraに保存されていました。

Matlabで読んでみました。

ちゃんと読めました。

 

加算平均波形の表示

エポックを一旦sample-epo.fifというファイル名で保存し読み込んだのち、加算平均をして波形を表示させます。

>>> epochs.save('c:\\Users\\akira\\sample-epo.fif') # データの書き出し

>>> saved_epochs=mne.read_epochs('C:\\Users\\akira\\sample-epo.fif') #データの読み込み

>>> print(saved_epochs) #データの表示

<EpochsFIF  |  n_events : 116 (all good), tmin : -0.199795213158 (s), tmax : 0.499488032896 (s), baseline : (-0.1997952163219452, 0.0), ~37.8 MB, data loaded,

 'aud_l': 55, 'aud_r': 61>

>>> evoked=epochs['aud_l'].average() # これで加算平均

>>> print(evoked)

<Evoked  |  comment : 'aud_l', kind : average, time : [-0.199795, 0.499488], n_epochs : 55, n_channels x n_times : 364 x 106, ~3.8 MB>

>>> evoked.plot()

<matplotlib.figure.Figure object at 0x00000207565F82B0>

ウインドウが開いて波形が表示されます。チャンネルの種類別になっています。

 

その他もろもろ

各エポックの最大値を求め方

>>> max_in_each_epoch=[e.max() for e in epochs['aud_l']]

>>> print(max_in_each_epoch[:5]) # 0,1,2,3,4番目のepochの最大値を求める

[1.9375167201631341e-05, 1.6405516986429127e-05, 1.8545377810380145e-05, 2.0412807568093327e-05, 2.0922971488786929e-05]

 

加算平均されたFIFファイルの読み方

加算平均されたFIFファイルを読みに行きます。名前が紛らわしいですがsample_audvis_ave.fifでなく、sample_audvis-ave.fifです。

左右の波形を表示させてみました。

>>> evoked_fname='c:\\Users\\akira\\mne_data\\MEG\\sample\\sample_audvis-ave.fif'

>>> evoked1=mne.read_evokeds(evoked_fname,condition='Left Auditory',baseline=(None,0),proj=True)

>>> evoked1.plot()

<matplotlib.figure.Figure object at 0x000002075685C978>

>>> evoked2=mne.read_evokeds(evoked_fname,condition='Right Auditory',baseline=(None,0),proj=True)

>>> evoked2.plot()

<matplotlib.figure.Figure object at 0x00000207566921D0>

 

左です。

右です。

 

差分波形の作成

>>> contrast=mne.combine_evoked([evoked1,evoked2],weights=[0.5,-0.5])

>>> contrast=mne.combine_evoked([evoked1,-evoked2],weights='equal') #でも同じ処理

>>> print(contrast)

<Evoked  |  comment : '0.500 * Left Auditory + -0.500 * Right Auditory', kind : average, time : [-0.199795, 0.499488], n_epochs : 116, n_channels x n_times : 376 x 421, ~4.8 MB>

>>> contrast.plot()

<matplotlib.figure.Figure object at 0x0000020758D339E8>

 

加算波形の加算

>>> average=mne.combine_evoked([evoked1,evoked2],weights='nave')

>>> print(contrast)

<Evoked  |  comment : '0.500 * Left Auditory + 0.500 * -Right Auditory', kind : average, time : [-0.199795, 0.499488], n_epochs : 116, n_channels x n_times : 376 x 421, ~4.8 MB>

>>> average.plot()

<matplotlib.figure.Figure object at 0x0000020755D97668>

 

エポック数を同じにした差分波形 trial-count equalization

>>> epochs_eq=epochs.copy().equalize_event_counts(['aud_l','aud_r'])[0]

>>> evoked1,evoked2=epochs_eq['aud_l'].average(),epochs_eq['aud_r'].average()

>>> print(evoked1)

<Evoked  |  comment : 'aud_l', kind : average, time : [-0.199795, 0.499488], n_epochs : 55, n_channels x n_times : 364 x 106, ~3.8 MB>

>>> print(evoked2)

<Evoked  |  comment : 'aud_r', kind : average, time : [-0.199795, 0.499488], n_epochs : 55, n_channels x n_times : 364 x 106, ~3.8 MB>

>>> contrast=mne.combine_evoked([evoked1,-evoked2],weights='equal')

>>> print(contrast)

<Evoked  |  comment : '0.500 * aud_l + 0.500 * -aud_r', kind : average, time : [-0.199795, 0.499488], n_epochs : 110, n_channels x n_times : 364 x 106, ~3.8 MB>

>>> contrast.plot()

<matplotlib.figure.Figure object at 0x0000020756895278>

 

時間周波数解析

Morletwavelet変換を使って時間周波数解析をします。まずMorletが使えるように準備します。

>>> import numpy as np

>>> n_cycles=2 # Morlet waveletは2周期

>>> freqs=np.arange(7,30,3) # 目的とする周波数帯域

>>> freqs

array([ 7, 10, 13, 16, 19, 22, 25, 28])

induced powerphase-locking valueを計算しMEG 1332の時間周波数マップを表示させます。

>>> from mne.time_frequency import tfr_morlet

>>> power,itc=tfr_morlet(epochs,freqs=freqs,n_cycles=n_cycles,return_itc=True,decim=3,n_jobs=1)

>>> power.plot([power.ch_names.index('MEG 1332')])

<matplotlib.figure.Figure object at 0x0000020762E36470>

 

空間フィルタの作成

FreeSurferで脳の切り出しが終わり、導体モデルが既にできていることを前提としています。

空間フィルタ作成のため、mne.minimum_normを導入し、ファイルを読み込みます。

>>> from mne.minimum_norm import apply_inverse,read_inverse_operator

>>> fname_inv='c:\\Users\\akira\\mne_data\\MEG\\sample\\sample_audvis-meg-oct-6-meg-inv.fif'

逆問題のための条件を設定していきます。

>>> inverse_operator=read_inverse_operator(fname_inv)

>>> print(inverse_operator)

<InverseOperator | MEG channels: 305 | EEG channels: 0 | Source space: Surface with 7498 vertices | Source orientation: Free>

>>> snr=3.0 # 信号雑音比

>>> lambda2=1.0 / snr**2

>>> method="dSPM"

電流源推定をし、ファイルに保存します。

>>> stc=apply_inverse(evoked,inverse_operator,lambda2,method)

>>> stc

<SourceEstimate  |  7498 vertices, subject : sample, tmin : -199.79521315838787 (ms), tmax : 499.488032896 (ms), tstep : 6.659840438612929 (ms), data size : 7498 x 106>

>>> stc.save('c:\\Users\\akira\\mne_dSPM_inverse')

mne_dSPM_inverse-lh.stcmne_dSPM_inverse-rh.stcの2つのファイルができました。

次に局所(label)での電流源の時系列変化を見ていきます。

まずは準備とlabelの読み込みです。

>>> fname_label='c:\\Users\\akira\\mne_data\\MEG\\sample\\labels\\Aud-lh.label'

>>> label=mne.read_label(fname_label)

>>> from mne.minimum_norm import apply_inverse_raw

015秒のdSPMを見てみます。

>>> start,stop=raw.time_as_index([0,15]) # 015秒の raw data

>>> stc=apply_inverse_raw(raw,inverse_operator,lambda2,method,label,start,stop)

>>> stc

<SourceEstimate  |  33 vertices, subject : sample, tmin : 0.0 (ms), tmax : 14991.3008273 (ms), tstep : 6.659840438612929 (ms), data size : 33 x 2252>

STCファイルに保存します。

>>> stc.save('c:\\Users\\akira\\mne_dSPM_raw_inverse_Aud')

mne_dSPM_raw_inverse_Aud-lh.stcmne_dSPM_raw_inverse_Aud-rh.stcの2つのファイルができました。