加算波形データを読み込む¶
In [1]:
import mne
filename='c:\\170207\\hashizume_akira\\041029\\SEF_median_ave.fif'
meg1=mne.read_evokeds(filename,proj=True) # SSPはかけておきます。
読み込んだデータの構造を見る¶
In [2]:
print(meg1)# 作成したmeg1が何かを調べる
type(meg1)
Out[2]:
[・・・]なのでリストです。<クラス名: 中身・・・・>となっています。
In [3]:
type(meg1[0])
Out[3]:
mne.evokedモジュール(?)のEvokedというクラスのリストです。
In [4]:
meg1[0] # meg1[0]って何？
Out[4]:
[・・・]がなくなりました。
In [5]:
print(dir(meg1[0])) # 属性attribute、関数methodを調べます。
_xxx_はクラス内部で使われる属性・関数です。外部からは利用できません。 どれが属性でどれが関数かわかりません。素直にMNE-pythonのホームページを参考にするとch_names,data,info,nave,kind,first,last,comment,times,verboseが属性でした。 printで表示される項目とは一致しないようです。
In [6]:
meg1[0].info #属栄infoを見る
Out[6]:
クラスInfoです。MNE-pythonのホームページでは属性を調べるとbads,ch_nams,chs,comps,custom_ref_applid,events,hpi_results,meas_date,nchan,projs,sfreq,acq_pars,acq_stim,buffer_size_sec,ch_head_t,description,dev_ctf_t,dev_head_t,dig,experimentor,file_id,highpass,hpi_meas,line_freq,lowpass,meas_id,proj_id,proj_name,subject_info,proc_historyとなっていました。
In [7]:
meg1[0].info["sfreq"] #標本化周波数
Out[7]:
In [8]:
sfreq=meg1[0].info["sfreq"]
print(meg1[0].first/sfreq,meg1[0].last/sfreq) # 時間幅を表示 printを使わないと2行表示できない？
print(meg1[0].times[0],meg1[0].times[-1]) # times[-1]はMATLABでtimes(end)
In [9]:
print(len(meg1[0].info["projs"])) # SSPの数
print(meg1[0].info["projs"][0])# 最初のSSPの中身
クラスProjectionの属性はMNE pythonのホームページででは記載されていませんでした。 関数にkeys()があります。これは辞書の項目を調べる関数でもあります。
In [10]:
ssp1=meg1[0].info["projs"][0]
ssp1.keys()
Out[10]:
In [11]:
ssp1["data"].keys()
Out[11]:
In [12]:
ch_names=ssp1["data"]["col_names"]
print(ch_names) # チャンネル名が表示されます。
In [13]:
x=ssp1["data"]["data"]
print(x) # SSPのベクトルが表示される
波形表示¶
In [14]:
chs=meg1[0].info["chs"] #チャンネル情報
type(chs)
Out[14]:
In [15]:
type(chs[0]) # 何のリスト
Out[15]:
In [16]:
chs[0] # 実際にchs[0]を見てみる
Out[16]:
チャンネル1個1個の細かい情報のようです。 チャンネルの種類別に分類します。
In [17]:
GRA,MAG,EEG,STIM=[],[],[],[] # MATLABだとGRA=[];MAG=[];EEG=[];STIM=[]
for num in range(0,len(chs)): # lenはMATLABのlength関数
str=chs[num]["ch_name"]
if str[0:3]=='MEG':
if str[-1]=='1': # str[-1]はMATLABだとstr(end)
MAG.append(num) # MATLABだとMAG=[MAG;num]
else:
GRA.append(num)
elif str[0:3]=='EEG':
EEG.append(num)
elif str[0:3]=='STI':
STIM.append(num)
In [18]:
import numpy as np
print(np.array(GRA).T) # リストを配列に変換して転置
チャンネル情報を一個一個みるのが面倒な時は、以下のようにします。
In [19]:
MEG2=np.array(range(0,306));
MAG2=np.array(range(3,306,3));
GRA2=np.setdiff1d(MEG2,MAG2);
print(GRA2);
In [20]:
%matplotlib inline
# %matplotlib inlineがないとjupyterでは表示されない
import matplotlib.pyplot as plt
# plt.close('all') # ウインドウは全部閉じる 対話式のjupyterでは不要
meg1[0].times.shape,meg1[0].data.shape # 配列の大きさをチェック
Out[20]:
553×1と315×553の配列なので転置が必要です。
In [21]:
plt.plot(meg1[0].times,meg1[0].data[GRA,:].T); # セミコロンで<matplotlib.lines.Line 2D at ...>を隠す
plt.xlim([meg1[0].times[0],meg1[0].times[-1]]); # x軸をトリミングしてます
plt.title('planar gradiometer');
plt.xlabel('[ms]');
plt.ylabel('[T/m]');
In [22]:
plt.plot(meg1[0].times,meg1[0].data[MAG,:].T,'b'); # セミコロンで<matplotlib.lines.Line 2D at ...>を隠す
plt.xlim([meg1[0].times[0],meg1[0].times[-1]]) ;# x軸をトリミングしてます
plt.title('magnetometer');
plt.xlabel('[ms]');
plt.ylabel('[T]');
In [23]:
meg1[0].plot(picks=GRA); # クラスEvokedはplot()という関数を持っています。これを使う方が簡単です。
meg1[0].plot(picks=MAG);
In [24]:
meg1[0].plot_topomap(times=[0.023,0.03,0.08],ch_type='grad'); # クラスEvokedの関数を使えばトポグラフィーも簡単に表示できます
meg1[0].plot_topomap(times=[0.023,0.03,0.08],ch_type='mag');
平面型グラジオメータの場合のトポグラフィーは等傾斜地場のトポグラフィーとなります。
基線の補正¶
In [25]:
times=meg1[0].times;
data=meg1[0].data.copy();
%matplotlib inline
fig=plt.figure();
plt.plot(times,data[GRA,:].T,'b');
plt.title('original')
plt.xlim([times[0],times[-1]]);
plt.xlabel('[ms]');
plt.ylabel('[T/m]');
t=np.nonzero((times>-0.05)*(times<0.0));# 時間幅は-0.05～0.0 基線補正用
sz=data.shape;
B=data[:,t];
print(B.shape)# 三次元配列になってしまう
data=data-np.dot(np.mean(B,axis=2),np.ones((1,sz[1])));
print(data.shape)
fig=plt.figure();
plt.plot(times,data[GRA,:].T,'b');
plt.title('after baseline correction')
plt.xlim([times[0],times[-1]]);
plt.xlabel('[ms]');
plt.ylabel('[T/m]');
周波数フィルタの適応¶
In [26]:
Fs=meg1[0].info["sfreq"];
times=meg1[0].times;
data2=mne.filter.band_pass_filter(x=data,Fs=Fs,Fp1=2,Fp2=100); # 2～100Hzです
%matplotlib inline
for num in range(0,2):
fig=plt.figure();
if num==0:
plt.plot(times,data[GRA,:].T,'b');
plt.title('after baseline correction');
else:
plt.plot(times,data2[GRA,:].T,'b');
plt.title('2-100 Hz');
plt.xlim([times[0],times[-1]]);
plt.xlabel('[ms]');
plt.ylabel('[T/m]');
時間幅の指定¶
In [27]:
%matplotlib inline
t=np.nonzero((times>-0.05)*(times<0.1));# 時間幅は-0.05～0.1
T=times[t];
for num in range(0,2):
fig=plt.figure();
if num==0:
X=data[GRA,:].copy(); # どうも2回にわける必要あり
X=X[:,t].squeeze(); # squeezeがないと三次元配列になってしまう
plt.plot(T,X.T,'b');
plt.title('after baseline correction');
else:
X=data2[GRA,:].copy();
X=X[:,t].squeeze();
plt.plot(T,X.T,'b');
plt.title('2-100 Hz');
plt.xlim([T[0],T[-1]]);
plt.xlabel('[ms]');
plt.ylabel('[T/m]');
プロットした頭皮の点表示¶
In [28]:
p=meg1[0].info["dig"] # Polhemusでプロットした点
type(p) # pって何？
Out[28]:
In [29]:
type(p[0]) # 何のリスト？
Out[29]:
In [30]:
p[0] # どんな辞書
Out[30]:
p[0]["r"]が座標のようです
In [31]:
import numpy as np
p=meg1[0].info["dig"] # Polhemusでプロットした点
pp=np.zeros((3,len(p))) # ゼロ行列の作成 MATLABだとzeros(3,length(p))
for num in range(0,len(p)):
pp[:,num]=p[num]["r"]
%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
# plt.close('all) # jupyterでは不要
fig=plt.figure()
ax=fig.add_subplot(111,projection='3d')
ax.scatter3D(pp[0,:],pp[1,:],pp[2,:])
ax.axis('equal') # 本来これでMATLABのaxis equalと同じになるはずだけど・・・bugです
x=[-0.15,0.15]
y=[-0.15,0.15]
z=[-0.05,0.15]
ax.auto_scale_xyz(x,y,z) # MATLABのdaspect([1,1,1])の代わり zは適当に調整してください
ax.set_axis_off()
ax.view_init(elev=0,azim=0) # 視点を変える
頭座標です x,y,z軸の長さは適当です。axis('equal')はx軸、y軸にしか使えず、z軸には使えません。bugです。 正確にやるにはグリッドのボックスをプロットした点の最小値～最大値にピッタリ合わせればよいらしいです。
In [32]:
print(dir(ax))
センサの三次元表示¶
In [33]:
%matplotlib inline
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
fig=plt.figure();
# ax=Axes3D(fig);
ax1=fig.add_subplot(1,2,1,projection='3d')
ax2=fig.add_subplot(1,2,2,projection='3d')
# ax=fig.gca(projection='3d')
wd,wz=0.014,0.0003
# p=np.array(zeros(3,len(MAG)))
one5=np.ones((5,1))
triangles=np.array([[0,1,2],[2,3,0]]) #,[0,3,2],[2,1,0]])
for num in range(0,len(MAG)):
num2=MAG[num]
loc=np.matrix([meg1[0].info["chs"][num2]["loc"]]) # リストを行列にする
p=one5*(loc[0,0:3]+loc[0,9:12]*wz) # 行列だと*が使える！ @でもよい
ex=loc[0,3:6].copy()*wd # copy()がなくてもこの場合はOK 基本「参照渡し」なので注意！
ey=loc[0,6:9].copy()*wd
p[0,:]+=ex+ey
p[1,:]=p[0,:]-2*ex
p[2,:]=p[1,:]-2*ey
p[3,:]=p[2,:]+2*ex
p[4,:]=p[0,:]
p=np.array(p) # 行列を配列に変換
x,y,z=p[:,0],p[:,1],p[:,2]
x,y,z=np.hstack(x),np.hstack(y),np.hstack(z) # 二次元配列を一次元配列にする
ax1.plot(x,y,z,color='blue')
ax2.plot_trisurf(x,y,z,triangles=triangles,color='aqua',alpha=1,edgecolor='none')
x=[-0.15,0.15]
y=[-0.15,0.15]
z=[-0.1,0.13]
ax1.auto_scale_xyz(x,y,z) # MATLABのdaspect([1,1,1])の代わり zは適当に調整してください
ax1.set_axis_off()
ax1.view_init(elev=30,azim=30) # 視点を変える
ax2.auto_scale_xyz(x,y,z) # MATLABのdaspect([1,1,1])の代わり zは適当に調整してください
ax2.set_axis_off()
ax2.view_init(elev=30,azim=30) # 視点を変える
線表示plotだけだとちゃんと表示されるのですが、三角メッシュ表示 plot_trisrfはうまく表示されませんでした。たぶんbugです。 頭座標でなく、センサ座標です。
In [34]:
%matplotlib inline
fig=meg1[0].plot_sensors(kind='3d',ch_type='mag');
ax=fig.axes[0]
ax.view_init(elev=30,azim=30); # 残念ながら回転しません！
クラスEvokedのセンサ表示関数plot_sensorsを使ってみました。jupyterでは回転しませんでした。
MATLABのファイル形式で保存¶
In [35]:
import scipy.io
out_data={} # 辞書として定義
type(out_data)
out_data["data"]=meg1[0].data
out_data["time"]=meg1[0].times
scipy.io.savemat('C:\\170207\\hashizume_akira\\041029\\test.mat',out_data) # MATLAB 7.3以降の形式で書き出し
MATLABでは以下のように読み込まれます。
load('C:\170207\hashizume_akira\041029\test.mat') whos Name Size Bytes Class Attributes
data 315x553 1393560 double
time 1x553 4424 double
ou_data["time"]はpythonでは553×1ですが、MATLABで1×553となっています。
In [36]:
del(out_data) # out_dataを消去
out_data=scipy.io.loadmat('C:\\170207\\hashizume_akira\\041029\\test.mat'); # MATファイルの読み込み
print(out_data.keys()) # 辞書の項目名を表示
print(out_data["data"].shape)
print(out_data["time"].shape)