emg(筋電図)の立ち上が(onset)りを適当なチャンネル(e.g. EEG063)から決めていくfunctionいろいろ. MRCFを想定しています。要 mne python , numpy, matplotlib. 下の方のコメント欄にも説明書いてます

detect_emg_onset.py

#!/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()

Readmeの代わり

筋電図の立ち上がり時刻を視察で決めてeventsの形式で出力する方法(function)3つ

1. マウスクリックなどの入った(通常STI101)トリガーチャンネルのデータを基に筋電図データの区間を取り出し、視察で立ち上がり時刻を取り出す

   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)

eventsのチェック用のfunction2つ

1. eventsでちゃんと筋電図の立ち上がりに合っているか確認します。

   check_plot_emg_events(raw, ch_emg, events)

2. eventsでちゃんとepochのトリガーに筋電図立ち上がりが合っているか確認します。

   check_plot_epochs_events(raw, ch_emg, events)

順次各項目の説明を追加します。

#1. マウスクリックなどで決めたトリガーを基に筋電図立ち上がり時刻を取り出す

筋電図のチャンネルと別にマウスクリックとか、外部刺激のトリガーとか(トリガーチャンネル)がデータに入っていて、その前後で筋電図が入るのがわかっている時に、トリガーを基準に筋電図を切り出して筋電図の立ち上がり時刻を視察で求めます。

式

events = detect_emg_onset_1ch_by_stim_ch(raw, ch_emg, stim_channel, event_id)

引数(入力するデータ)

• raw: mne.io.RawFIFのインスタンス。
• ch_emg: 筋電図のチャンネル1個 (e.g. 'EEG063')。 raw.ch_namesで出てくるチャンネルでチャンネルは一覧できます。
• stim_channel: トリガーチャンネル(通常'STI101')
• event_id: マウスリリースなど手の動きを示すトリガーの番号

返り値(戻り値、出力するデータ)

• events：mne.find_eventsで出てくる形のnumpy array