脳磁図データの表示（evoked data）
加算波形データを表示する方法です。 いろんな方法があると思いますがaxes１個のものにしてみました。 データは広島大学大学院医歯薬保健学研究院生体環境適応科学研究室（弓削研究室）の 大鶴先生から頂いた体性感覚刺激の２つのFIFFファイルをMATファイルに変換したものです。 被験者名と検査日は削除しています。 データはhttp://meg.aalip/matlab/data/SEF1.matとhttp://meg.aalip.matlab/data/SEF2.matから入手できます。 vv_lout.matはレイアウト表示用のファイルです。
Contents
脳磁図データとレイアウトデータの読み込み
ファイルはそれぞれ
http://meg.aalip.jp/matlab/data/SEF1.mat
http://meg.aalip.jp/matlab/data/SEF2.mat
http://meg.aalip.jp/matlab/data/vv_lout.mat
から入手してください。 データファイルが隣のフォルダ、../dataフォルダにあると仮定しています。 データですが、頭とDewarの位置関係は変わらず、標本化周波数は同一と仮定しています。
load('../data/SEF1'); SEF1=ave_0001; MTX1=MEG2HEAD; load('../data/SEF2'); SEF2=ave_0001; MTX2=MEG2HEAD; load('../data/vv_lout');
SSP適応
Signal Space Projectionを適応します。 PCA_xxxxのチャンネル順とデータのチャンネル順が同じと仮定しています。 またSEF1とSEF2のSSPが同一のものであると仮定しています。 SSPの効果が一目瞭然と思います。
if exist('PCA_0001','var') SSP=eye(length(PCA_0001)); for n=1:999 nPCA=sprintf('PCA_%04.0f',n); if exist(nPCA,'var') nPCA=eval(nPCA);% 文字列を実際の固有ベクトルに SSP=SSP-nPCA*nPCA'; else break; end; end; ch=1:length(PCA_0001); figure('color',[1,1,1]); subplot(221);plot(SEF1(:,end),SEF1(:,ch),'color',[0,0,1]); title('SSP前'); subplot(222);plot(SEF2(:,end),SEF2(:,ch),'color',[0,0.5,0]); title('SSP前'); SEF1(:,ch)=SEF1(:,ch)*SSP; SEF2(:,ch)=SEF2(:,ch)*SSP; subplot(223);plot(SEF1(:,end),SEF1(:,ch),'color',[0,0,1]); title('SSP後'); subplot(224);plot(SEF2(:,end),SEF2(:,ch),'color',[0,0.5,0]); title('SSP後'); haxes=findobj(gcf,'type','axes');% 単位記入 for n=1:length(haxes); set(gcf,'CurrentAxes',haxes(n)); axis tight; xlabel('[ms]');ylabel('[fT or fT/cm]'); end; end;
2～200Hzの周波数帯域通過フィルタの作成
最小二乗線形位相 FIR フィルターであるfirls関数を使って周波数フィルタを設計し、 周波数応答を見てみます。
sfreq=1000/(SEF1(2,end)-(SEF1(1,end))); bandpass=[2,200]; % 通過周波数帯域 f=[0,1,3,199,201,sfreq/2]/(sfreq/2); a=1;% FIRなので分母は1 n=256;% 次数は128 bp=[0,0,1,1,0,0]; b=firls(n,f,bp); [h,w]=freqz(b,a,1000,sfreq); figure('color',[1,1,1]); subplot(211);plot(w,abs(h));grid on; set(gca,'xlim',[0,sfreq/2]); xlabel('[Hz]'); subplot(212);plot(w,abs(h));grid on; set(gca,'xlim',[0,10]); xlabel('[Hz]');
周波数フィルタの適応
低周波帯域のフィルタ特性はかなり甘甘です。 実際にフィルタを適応してみます。位相がズレないようfiltfilt関数を使います。
figure('color',[1,1,1]); for n=1:2 switch n case 1;data=SEF1;color=[0,0,1]; case 2;data=SEF2;color=[0,0.5,0]; end; time=data(:,end);ch=1:306; subplot(2,2,1+(n-1));plot(time,data(:,ch),'color',color); axis tight;grid on;title('フィルタ前');xlabel('[ms]'); data(:,ch)=filtfilt(b,a,data(:,ch)); subplot(2,2,3+(n-1));plot(time,data(:,ch),'color',color); axis tight;grid on;title('フィルタ後');xlabel('[ms]'); switch n case 1;SEF1=data; case 2;SEF2=data; end; end;
全波形表示
gradiometerのMEGxxx2とMEGxxx3は対応してないかもしれません。 MEGxxx1はmagnetometerです。magnetomterとgradiometerで単位が異なるので補正してます。
time=union(SEF1(:,end),SEF2(:,end));% ２つの時間帯を統合 lout=vv_lout; wd=lout(1,4);% 波形の幅 ht=lout(1,5);% 波形の高さ ch=1:306; chlout=1:306;% layout用 amp=0.01*ht; amp102=0.2;% magnetometerは単位が違うので5分の1で表示 X0=linspace(0,wd,length(time))'; for n=1:2; switch n case 1;data=SEF1; hf=figure('color',[1,1,1]); color=[0,0,1]; case 2;data=SEF2; color=[0,0.5,0]; end; [~,t]=intersect(time,data(:,end));% 統合した時間のなかでdataの時間部分を表示 data=amp*data(:,ch); data(:,3:3:306)=amp102*data(:,3:3:306);% magnetometerだけ補正 smp=length(t); X=ones(smp,1)*lout(chlout,2)'+X0(t,:)*ones(1,size(data,2)); Y=ones(smp,1)*lout(chlout,3)'+data; plot(X,Y,'color',color);hold on; end; set(gca,'Ydir','normal'); axis tight;axis off; title('全波形表示');
gradiometer波形の表示
MEGxxx1の成分、magnetometerの成分を除去してます。
time=union(SEF1(:,end),SEF2(:,end));% ２つの時間帯を統合 lout=vv_gra; wd=lout(1,4);% 波形の幅 ht=lout(1,5);% 波形の高さ ch=1:306; ch(3:3:306)=[];% magnetometerは除去 chlout=1:204;% layout用 amp=0.01*ht; X0=linspace(0,wd,length(time))'; for n=1:2; switch n case 1;data=SEF1; hf=figure('color',[1,1,1]); color=[0,0,1]; case 2;data=SEF2; color=[0,0.5,0]; end; [~,t]=intersect(time,data(:,end));% 統合した時間のなかでdataの時間部分を表示 data=amp*data(:,ch); smp=length(t); X=ones(smp,1)*lout(chlout,2)'+X0(t,:)*ones(1,size(data,2)); Y=ones(smp,1)*lout(chlout,3)'+data; plot(X,Y,'color',color);hold on; end; set(gca,'Ydir','normal'); axis tight;axis off; title('gradiometer波形表示');
gradiometer波形のroot mean square波形表示
MEGxxx2とMEGxxx3のroot mean square波形を表示します。magnetometerのlayoutを使います。
time=union(SEF1(:,end),SEF2(:,end));% ２つの時間帯を統合 lout=vv_mag; wd=lout(1,4);% 波形の幅 ht=lout(1,5);% 波形の高さ ch=1:306; ch(3:3:306)=[];% magnetometerは除去 chlout=1:102;% layout用 amp=0.01*ht; X0=linspace(0,wd,length(time))'; for n=1:2; switch n case 1;data=SEF1; hf=figure('color',[1,1,1]); color=[0,0,1]; case 2;data=SEF2; color=[0,0.5,0]; end; [~,t]=intersect(time,data(:,end));% 統合した時間のなかでdataの時間部分を表示 data=amp*data(:,ch); smp=length(t); data=sqrt(data(:,1:2:end).^2+data(:,2:2:end).^2); X=ones(smp,1)*lout(chlout,2)'+X0(t,:)*ones(1,size(data,2)); Y=ones(smp,1)*lout(chlout,3)'+data; plot(X,Y,'color',color);hold on; end; set(gca,'Ydir','normal'); axis tight;axis off; title('gradiometer波形のroot mean square波形表示');
