epochでicaをかけた後、各componentのマップと加算平均波形を表示するfunction: plot_comp_map_and_evoked のスクリプト

epoch_ica_plot_component.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# function設定(componentの任意の条件(neutral/centerとか)とcomponentの番号から
#  mapとevoked(componentから作成)のプロット
def plot_comp_map_and_evoked(comps, cnd, idx_comp, is_peak_need):
	"""
	component (icaの)を各条件、component番号で加算平均波形をプロットする。

	comps: icaで出てくるcomponentのデータ. e.g.comps=ica.get_sources(epochs)
	cnd: comps.event_idで出てくる条件の名前のどれか
	idx_comp: ica.plot_components()で出てくるcomponentの番号のどれか
	"""
	comp_data = comps[cnd].get_data()
	evoked_comp_data = np.mean(comp_data[:,idx_comp,:],axis=0)
	times = comps.times*1000 # ms
	fig_ic = ica.plot_components(picks=[idx_comp,0])
	ax_evo=fig_ic.axes[1] #右側のaxes
	fig_ic.delaxes(ax_evo) #右側消す
	fig_ic.add_subplot(122) #右側再度axes設定
	plot(times, evoked_comp_data)
	title('evoked component')
	xlim((times[0],times[-1]))
	suptitle(cnd + ', component no.:'+str(idx_comp))
	# peakの潜時が欲しいときは以下
	if is_peak_need:
		coord = plt.ginput(1)
		while len(coord)>0:
			lat = coord[0][0]
			seg_evo = evoked_comp_data[(times > lat-10) & (times < lat+10)] # evoked segment
			seg_idx = np.where((times > lat-10) & (times < lat+10))[0] # index segment

			# is convex or concave 指した場所が上に凸のピークか判断
			is_convex = np.mean(seg_evo - np.linspace(seg_evo[0],seg_evo[-1],num=len(seg_evo))) > 0
			if is_convex:
				idx_peak_seg = np.min(np.where(seg_evo == max(seg_evo)))
			else:
				idx_peak_seg = np.min(np.where(seg_evo == min(seg_evo)))
			idx_peak = idx_peak_seg + seg_idx[0]
			# plot
			plt.cla()
			plot(times, evoked_comp_data)
			plot(times[idx_peak],evoked_comp_data[idx_peak],'r*')
			xlim((times[0],times[-1]))
			title('%s %0.1f %s' %('latency:',times[idx_peak],'ms'))
			coord = plt.ginput(1)
	return fig_ic