Last active
December 20, 2015 00:40
-
-
Save cosacog/1513852da9d37680693b to your computer and use it in GitHub Desktop.
emg(筋電図)の立ち上が(onset)りを適当なチャンネル(e.g. EEG063)から決めていくfunctionいろいろ. MRCFを想定しています。要 mne python , numpy, matplotlib. 下の方のコメント欄にも説明書いてます
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python | |
# -*- coding: utf-8 -*- | |
# この中で実際に使うのは | |
# 1. マウスクリックなどで決めたトリガーを基に取り出す | |
# detect_emg_onset_1ch_by_stim_ch(raw, ch_emg, stim_channel, event_id) | |
# 2. 視察で筋電図から大まかにピークを取り出して、次に再度各区間で立ち上がり時刻を正確に決める | |
# detect_emg_onset_1ch_vis_inspect(raw, ch_emg) | |
# 3. 閾値で筋電図からピークを検出して、次に再度各区間で立ち上がり時刻を正確に決める | |
# detect_emg_onset_1ch_by_thr(raw, ch_emg, win_mv) # e.g. win_mv=0.3 | |
# のどれかです. | |
# eventsのチェック用のfunctionを追加しました。 | |
# 4. eventsでちゃんと筋電図が立ち上がりが合っているか確認します。 | |
# check_plot_emg_events(raw, ch_emg, events) | |
# 5. eventsでちゃんとepochのトリガーに筋電図立ち上がりが合っているか確認します。 | |
# check_plot_epochs_events(raw, ch_emg, events) | |
# その他いくつかfunctionがありますが、上記1-3のfunctionから呼び出すためのものです | |
def raw2emg(raw, ch_emg, freq_hp): | |
""" | |
rawインスタンスから筋電図のチャンネルを取り出してhigh passフィルタをかけ、 | |
整流したものを出力 | |
:param raw: mne.io.RawFIFでできるrawインスタンス | |
:param ch_emg: 筋電図の入ったチャンネル. e.g.'EEG063' | |
:param freq_hp: 基線の揺れを抑えるためハイパスをかける周波数 (Hz), | |
e.g. 3とか1とかくらいがよろしいかと。5以上は歪みが大きくなる懸念を感じます | |
:return: emg_rect_data: 筋電図のデータを取り出したnumpy array | |
times: 時間を示すnumpy array | |
emg: rawからch_emgのチャンネルに限定したrawインスタンス | |
""" | |
import numpy as np | |
# EMG取り出し | |
emg = raw.copy() | |
# ch_emgがlistでなかったらlistにする | |
if type(ch_emg) is not list: | |
ch_emg = [ch_emg] | |
emg.pick_channels(ch_emg) | |
emg.filter(freq_hp,None) # 5 hz high pass | |
emg_data, times = emg[:,:] | |
# rectification | |
emg_rect_data = np.abs(emg_data) | |
# データが1行だったら1列に変更(転置) | |
if emg_rect_data.shape[0]==1: | |
emg_rect_data = emg_rect_data.T | |
return(emg_rect_data, times, emg) | |
def plot_ginput(emg, t_peaks, dur_seg, event_id): | |
""" | |
emgのデータとおよその立ち上がり時刻で4秒程度区間を切り出して筋電図の正確な立ち上がり時刻を目視で取り出す | |
他のfunctionから呼び出すのに使う | |
:param emg: rawからemgのチャンネルを取り出したインスタンス | |
:param t_peaks: おおよその筋電図立ち上がり時刻(index) | |
:param dur_seg: t_peaksを基準に前後取り出す時間幅(s) | |
:param event_id: 出力するeventsで使うトリガー番号 | |
:return:events_jerk: 出力するevents (mne.find_eventsと同じ形式のnumpy array) | |
""" | |
import numpy as np | |
import matplotlib.pyplot as plt | |
emg_data, times = emg[:,:] #データをrawインスタンスから取り出し | |
# データが1行だったら1列に変更(転置) | |
if emg_data.shape[0]==1: | |
emg_data = emg_data.T | |
sf = emg.info['sfreq'] | |
# 視察区間のサンプル数 | |
dur_idx_pre = int(dur_seg[0]*sf) # 前 | |
dur_idx_pst = int(dur_seg[1]*sf) # 後 | |
vec_t_emg_onset =[] #これに立ち上がり時刻(配列)を入れていく | |
for ii in np.arange(len(t_peaks)): | |
## 解析区間のindex (idx)取り出し | |
# idxの最初が0より小くならないようにする | |
idx_init_tmp_vis = np.max([0,t_peaks[ii]-dur_idx_pre]) | |
# idxの最後が区間のサイズより大きくならないようにする | |
idx_end_tmp_vis = np.min([len(emg_data),t_peaks[ii]+dur_idx_pst]) | |
# 視察の区間のindexを配列にする | |
idx_tmp_vis = np.arange(idx_init_tmp_vis,idx_end_tmp_vis).astype(int) | |
plt.clf() | |
# emgを区間で切り取って整流 | |
emg_seg = np.abs(emg_data[idx_tmp_vis]-np.median(emg_data[idx_tmp_vis])) | |
times_seg = times[idx_tmp_vis] | |
plt.plot(times_seg, emg_seg) | |
plt.xlim((times_seg[0],times_seg[-1])) | |
plt.title(str(ii+1)+'/'+ str(len(t_peaks))) | |
# マウスで立ち上がり時刻を決める | |
coord = plt.ginput(1,timeout=300) | |
t_onset_vis = np.array([]) # クリックしなかった場合の判定のために設定 | |
# 画面に縦線(潜時)を表示: Enterを押すとloopを抜ける | |
while coord: | |
t_onset_vis = coord[0][0] | |
plt.clf() | |
plot(times_seg, emg_seg) | |
xlim((times_seg[0],times_seg[-1])) | |
title(str(ii+1)+'/'+ str(len(t_peaks))) | |
axvline(t_onset_vis, color='red') | |
coord = ginput(1,timeout=300) | |
# マウスクリックしないでescを押すとt_onset_vis.size => 0になる | |
if t_onset_vis.size > 0: | |
vec_t_emg_onset = np.r_[vec_t_emg_onset, t_onset_vis] | |
plt.close() | |
# 時刻をindexに変換-> eventsに変換 | |
events_jerk = zeros((len(vec_t_emg_onset),3),dtype=int) # events | |
for ii in np.arange(len(vec_t_emg_onset)): | |
events_jerk[ii,0]=min(where(times > vec_t_emg_onset[ii])[0]) | |
events_jerk[:,0]=events_jerk[:,0]+raw.first_samp # 補正しないとうまく動かない | |
# eventsの3列目: トリガーの番号 | |
events_jerk[:,2]= event_id | |
return events_jerk | |
def detect_emg_onset_1ch_by_stim_ch(raw, ch_emg, stim_channel, event_id): | |
""" | |
トリガーのチャンネルから筋電図の時間を大まかに切り出しし、各区間で筋電図立ち上がり時刻をマウスで合わせる | |
:param raw: mne.io.RawFIFのインスタンス | |
:param ch_emg: raw.ch_namesで出てくるチャンネルのうち筋電図のチャンネル1個 (e.g. 'EEG063') | |
:param stim_channel: raw.ch_namesで出てくるチャンネルのうちトリガーのチャンネル(通常'STI101') | |
:param event_id: マウスリリースなど手の動きを示すトリガーの番号 | |
:return:events:mne.find_eventsで出てくる形のnumpy array | |
""" | |
import mne | |
import numpy as np | |
## EMG取り出し: emg_rect_data:, emg | |
freq_highpass =1 # Hz | |
emg_rect_data, times, emg = raw2emg(raw, ch_emg, freq_highpass) | |
## stim_channelからevents取り出し | |
events_all = mne.find_events(raw, stim_channel=stim_channel) | |
# マウスのボタン押しなどに対応するeventsだけ取り出す | |
events = events_all[events_all[:,2]==event_id,:] | |
t_peaks = events[:,0] - raw.first_samp # onset index | |
## 筋電図の立ち上がりを視察で同定 | |
# 設定 | |
dur_seg = [2,2] # sec: 筋電図を観察するための区間(前後dur_seg秒の幅だけ表示する) | |
event_id_out =1 # 出力するevents_jerkで使うトリガー番号 | |
# emgデータとおおよその立ち上がり時刻から正確なemg立ち上がり時刻(のevents)を取り出す | |
events_jerk = plot_ginput(emg, t_peaks, dur_seg,event_id_out) | |
return(events_jerk) | |
def detect_emg_onset_1ch_vis_inspect(raw, ch_emg): | |
""" | |
筋電図のチャンネルから立ち上がり時刻=運動開始時刻を決めてeventsで出力するfunction | |
最初にEMG全体の画面が出てきて、視察-クリックで筋電図の各頂点を同定する | |
Enterを押すと各頂点を含む区間の筋電図が表示されるのでクリックで立ち上がり時刻を決めてEnter, | |
無視するときはクリックする前にEscを押す | |
ファイルにsaveするときはnp.save('events.npy')あたりで | |
c.f. detect_emg_onset_1ch_by_thr | |
:param raw: mne.io.RawFIFで出てくるrawインスタンス | |
:param ch_emg: 筋電図の入ったチャンネル1個 (e.g. 'EEG063') | |
:return:events_jerk: events= mne.find_eventsで出てくる形式のnumpy array | |
""" | |
import mne | |
import numpy as np | |
import matplotlib.pyplot as plt | |
## EMG取り出し | |
freq_highpass =1 # Hz | |
emg_rect_data, times, emg = raw2emg(raw, ch_emg, freq_highpass) | |
## visual inspection for detecting peaks | |
plt.figure() | |
plt.plot(emg_rect_data) | |
plt.title('Click to set the time for onset detection') | |
t_peaks = [] # time of peaks to detect onset | |
t_tmp =ginput(1) | |
while t_tmp: | |
plt.plot(t_tmp[0][0],t_tmp[0][1],'r*') | |
t_peaks = np.r_[t_peaks,t_tmp[0][0]] | |
t_tmp = ginput(1,timeout=300) | |
t_peaks = np.round(np.sort(t_peaks)) | |
## 筋電図の立ち上がりを視察で同定 | |
# 設定 | |
dur_seg = [2,2] # sec: 筋電図を観察するための区間(前後dur_seg秒の幅だけ表示する) | |
event_id_out =1 # 出力するevents_jerkで使うトリガー番号 | |
# emgデータとおおよその立ち上がり時刻から正確なemg立ち上がり時刻(のevents)を取り出す | |
events_jerk = plot_ginput(emg, t_peaks, dur_seg,event_id_out) | |
return(events_jerk) | |
def detect_emg_onset_1ch_by_thr(raw, ch_emg, win_mv): | |
""" | |
raw: mne.ioRawFIFのオブジェクト | |
ch_emg: emgのチャンネル e.g.['EEG063'] | |
筋電図のチャンネルから立ち上がり時刻=運動開始時刻を決めてeventで出力するfunction | |
最初にEMG全体の画面が出てきて、クリックで閾値を決める. | |
検出された頂点の数はタイトルに表示 | |
Enterを押すと各頂点が表示されるのでクリックで立ち上がり時刻を決めてEnter, | |
無視するときはクリックする前にEscを押す | |
ファイルにsaveするときはnp.save('events.npy')あたりで | |
c.f detect_emg_onset_1ch_vis_inspect | |
:param raw:raw instance | |
:param ch_emg: EMG channel. e.g. ['EEG063'] | |
:param:win_mv: 包絡線(envelope)を作るときの時間窓(sec), 0.1-0.3くらいがいいぽい | |
:return:events_jerk: mne.find_eventsで出てくるeventsと同じ形式のnumpy array | |
""" | |
import mne | |
import numpy as np | |
import matplotlib.pyplot as plt | |
# EMG取り出し | |
freq_highpass =1 # Hz | |
emg_rect_data, times, emg = raw2emg(raw, ch_emg, freq_highpass) | |
## visual inspection for detecting peaks -> threshold to detect peaks | |
sf = emg.info['sfreq'] | |
# win_mv = 0.1 # sec | |
n_pnt_win = np.int(sf*win_mv/2) # points for moving window | |
emg_env = np.zeros(emg_rect_data.shape) # envelope | |
print('estimating envelope. it takes some time ...') | |
for ii in np.arange(len(emg_env)): | |
idx_seg = np.arange(ii-n_pnt_win, ii+n_pnt_win) | |
idx_seg_valid = idx_seg[(idx_seg >=0) & (idx_seg <len(emg_env))] | |
emg_env[ii] = np.max(emg_rect_data[idx_seg_valid]) | |
plt.figure() | |
plt.plot(times,emg_env) | |
plt.xlim((times[0],times[-1])) | |
plt.title('Click to set the threshold for onset detection') | |
#t_peaks = [] # time of peaks to detect onset | |
thr_tmp =plt.ginput(1) | |
while thr_tmp: | |
plt.clf() | |
plt.plot(times,emg_env) | |
thr = thr_tmp[0][1] | |
plt.axhline(thr,color='red') | |
plt.xlim((times[0],times[-1])) | |
# 区間決め | |
emg_bin = emg_env > thr | |
idx_emg_bin_onset = np.array(where((emg_bin[1:] ==True) & (emg_bin[:-1]==False)))[0] +1 # onset | |
idx_emg_bin_ofset = np.array(where((emg_bin[1:] ==False)& (emg_bin[:-1]==True)))[0] # offset | |
# discard if the last onset is not back to baseline | |
if idx_emg_bin_onset[-1] > idx_emg_bin_ofset[-1]: | |
idx_emg_bin_onset = idx_emg_bin_onset[:-1] | |
times_onset = times[idx_emg_bin_onset] | |
plt.plot(times_onset,ones(len(times_onset))*thr,'r*') | |
plt.title('cross points '+str(len(times_onset))) | |
thr_tmp = ginput(1,timeout=300) | |
#t_peaks = np.round(np.sort(t_peaks)) | |
t_peaks = idx_emg_bin_onset # onset index | |
## 筋電図の立ち上がりを視察で同定 | |
# 設定 | |
dur_seg = [2,2] # sec: 筋電図を観察するための区間(前後dur_seg秒の幅だけ表示する) | |
event_id_out =1 # 出力するevents_jerkで使うトリガー番号 | |
# emgデータとおおよその立ち上がり時刻から正確なemg立ち上がり時刻(のevents)を取り出す | |
events_jerk = plot_ginput(emg, t_peaks, dur_seg,event_id_out) | |
return(events_jerk) | |
def check_plot_emg_events(raw, ch_emg, events): | |
""" | |
筋電図から取り出した立ち上がり時刻が正しいか確認するplot | |
numpy arrayを重ね書きして表示 | |
0の時刻で筋電図が立ち上がっていればok.もし違っていれば処理を確認. | |
:param raw: mne.io.RawFIFでできるrawインスタンス | |
:param ch_emg: 筋電図の入ったチャンネル. e.g.'EEG063' | |
:param events: mne.find_eventsでできるeventsと同じ形式のnumpy array | |
:return: なし | |
""" | |
import mne | |
import matplotlib.pyplot as plt | |
event_id = int(events[0,2]) | |
epochs = mne.Epochs(raw, events, event_id=event_id,tmin=-2.0, tmax=2.0, | |
preload=True) | |
# ch_emgがlistでなければlistにする | |
if type(ch_emg) is not list: | |
ch_emg = [ch_emg] | |
epochs.pick_channels(ch_emg) | |
epochs_data = epochs.get_data() | |
times = epochs.times | |
plt.figure() | |
plt.plot(times, epochs_data[:,0,:].T) | |
title('EMG plot') | |
def check_plot_epochs_events(raw, ch_emg, events): | |
""" | |
筋電図から取り出した立ち上がり時刻が正しいか確認するplot | |
mne.Epochsのインスタンスのプロット | |
0の時刻で筋電図が立ち上がっていればok.もし違っていれば処理を確認. | |
:param raw: mne.io.RawFIFでできるrawインスタンス | |
:param ch_emg: 筋電図の入ったチャンネル. e.g.'EEG063' | |
:param events: mne.find_eventsでできるeventsと同じ形式のnumpy array | |
:return: なし | |
""" | |
import mne | |
import matplotlib.pyplot as plt | |
event_id = int(events[0,2]) | |
epochs = mne.Epochs(raw, events, event_id=event_id,tmin=-2.0, tmax=2.0, | |
preload=True) | |
# ch_emgがlistでなければlistにする | |
if type(ch_emg) is not list: | |
ch_emg = [ch_emg] | |
epochs.pick_channels(ch_emg) | |
epochs.plot() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
#1. マウスクリックなどで決めたトリガーを基に筋電図立ち上がり時刻を取り出す
筋電図のチャンネルと別にマウスクリックとか、外部刺激のトリガーとか(トリガーチャンネル)がデータに入っていて、その前後で筋電図が入るのがわかっている時に、トリガーを基準に筋電図を切り出して筋電図の立ち上がり時刻を視察で求めます。
式
引数(入力するデータ)
返り値(戻り値、出力するデータ)