*fifファイルを読み込んでdata.gradの情報が取り出せないことがある

ft_preprocessing, ft_read_headerとかでとれてくるデータ

例. hdr=ft_read_header('hogehoge.fif')

では通常

hdr.gradの情報が含まれてますが、ないことがあります。gradがないとforward solution (lead field行列の計算=ft_prepare_leadfield)とかで処理ができずエラーが返ってきます。

原因

脳波のチャンネルが含まれるときgradが返ってこないことがあります (おそらく脳波が60ch以上あるとき)。

MEGのチャンネルだけだと問題が生じません。

脳波が含まれていてもgradが返ってきたり、来なかったりするのは、脳波のチャンネルが60を越えているかどうかによると推測します。(下記スクリプトの内容から)

大まかな説明

gradを取り出すメインの関数はmne2grad.mのようです(fieldtrip-fileio-privateの中)

下の方に130行くらいから最後までスクリプトの一部を張り付けてます。

numel(dig_eeg)<numel(ch_eeg)の状態になり、forループがエラーになるのが原因でした。

mne2grad.mはft_read_headerから呼び出すとき、try~catchでくくられてるのでエラーが出なかったようです。

もう少し説明

orig=fiff_read_meas_info('hogehoge.fif')で取り出された情報で、

dig_eegはデジタイザから取り出したEEGチャンネルの情報

chn_eegはチャンネル情報(orig.chsにはch_nameとかいろいろ含まれる)から取り出したEEGチャンネルの情報

を示しています。

デジタイザからの情報であるdig_eegにはreferenceチャンネルが含まれる(1chがreference)のに対し、

チャンネル情報のch_eegにはreferenceが含まれていない、

という前提でスクリプト中のdig_eeg(1)=[];で消されているようです。

neuromagのデータは元からdig_eegにreferenceの情報は含まれておらず(たぶん)、

dig_eeg(1)=[];(下記スクリプト参照)でnumel(dig_eeg)<numel(ch_eeg)になるため、forループがエラーになり、結果としてgradが返ってこなくなるようです。

エラーが出力されればむしろ原因がつかみやすかったですが、mne2grad.mがエラーになっても、それを呼び出すft_read_headerが

try

mne2grad(hoge)

catch

云々

end

でtry~catch~endでくくっているためエラーが出力されず、単にhdr.gradが含まれていないという結果になっていたようです。

対策

1. 下のスクリプト中の黄色の部分を追記します。chn_eeg云々の行は、元々は黄色のif~endブロックより下にあるので、エラーにならないよう前に持ってきます。

2. 単に下のスクリプト中のオレンジの行をコメントアウトしてもよいと思います。

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ mne2grad.m の一部 (だいたい130行目くらいから最後まで ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

if nEEG>0

elec.pnt = zeros(nEEG,3);

elec.unit = 'cm';

elec.label = cell(nEEG,1);

% Amendments to overcome problem with fiff_read_meas_info.m when >60

% channels (thanks to Rik Henson)

dig_eeg = find([orig.dig.kind]==3); % Find EEG digitisations

chn_eeg = find([orig.chs.kind]==2); % Find EEG channels

if ~isempty(dig_eeg)

if numel(dig_eeg)>numel(chn_eeg) %<-追加

dig_eeg(1) = []; % Remove reference

end %<-追加

end

%disp(dig_eeg);%チェック

k=0;

for n = chn_eeg

k=k+1;

if nEEG<=60

elec.pnt(k,1:3)=100*orig.chs(n).eeg_loc(1:3); % multiply by 100 to get cm

else

elec.pnt(k,1:3)=100*orig.dig(dig_eeg(k)).r; % multiply by 100 to get cm

end

elec.label{k}=deblank(orig.ch_names{n});

end

%check we've got all the EEG channels:

if k ~= (nEEG)

error('Number of EEG channels identified does not match number of channels in elec structure!!!!!');

end

end