vb_gen_megdata

 シミュレーション用MEGデータを作成する
 2003-11-14 Taku Yoshioka

 parm.[dir/file].brain : 全脳モデルファイル名
 parm.[dir/file].basis : リードフィールド行列ファイル名
 MEGinfo.nprsamp   : ノイズサンプルステップ数
 MEGinfo.Noise     : ノイズモデル(VB推定における同パラメータ参照)
 MEGinfo.Noise_file: ノイズ共分散ファイル名
 MEGinfo.T         : MEGデータステップ数(時間)
 MEGinfo.SN        : 観測ノイズSN比
 MEGinfo.ix0       : 電流源中心インデックス(電流源数x1)
 MEGinfo.ix        : 電流源インデックス(電流源数x1、セル配列)
 MEGinfo.pattern   : 電流パターン情報(電流源数x1、セル配列)
 MEGinfo.pattern{i}: (発火回数x3、発火開始・終了時刻、強度)
 MEGinfo.Rfilt     : 平滑化フィルタ半径(Rmax=2*Rfilt)

 注)TはMEG信号の全時間長(ステップ数)であり、
 その中の1:nprsampは観測ノイズとみなされる。
 nprsampの値は作成されるMEG信号データに何の影響も与えない。

 bexp0   : ノイズ無しMEGデータ
 bexp    : 作成されたMEGデータ(bexp0+ノイズ)
 J       : MEGデータに対応する電流分布
 Iextract: Jの各行に対応する全脳モデル頂点インデックス


 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [bexp0,bexp,J,Iextract] = vb_gen_megdata(parm,MEGinfo)
0002 %
0003 % シミュレーション用MEGデータを作成する
0004 % 2003-11-14 Taku Yoshioka
0005 %
0006 % parm.[dir/file].brain : 全脳モデルファイル名
0007 % parm.[dir/file].basis : リードフィールド行列ファイル名
0008 % MEGinfo.nprsamp   : ノイズサンプルステップ数
0009 % MEGinfo.Noise     : ノイズモデル(VB推定における同パラメータ参照)
0010 % MEGinfo.Noise_file: ノイズ共分散ファイル名
0011 % MEGinfo.T         : MEGデータステップ数(時間)
0012 % MEGinfo.SN        : 観測ノイズSN比
0013 % MEGinfo.ix0       : 電流源中心インデックス(電流源数x1)
0014 % MEGinfo.ix        : 電流源インデックス(電流源数x1、セル配列)
0015 % MEGinfo.pattern   : 電流パターン情報(電流源数x1、セル配列)
0016 % MEGinfo.pattern{i}: (発火回数x3、発火開始・終了時刻、強度)
0017 % MEGinfo.Rfilt     : 平滑化フィルタ半径(Rmax=2*Rfilt)
0018 %
0019 % 注)TはMEG信号の全時間長(ステップ数)であり、
0020 % その中の1:nprsampは観測ノイズとみなされる。
0021 % nprsampの値は作成されるMEG信号データに何の影響も与えない。
0022 %
0023 % bexp0   : ノイズ無しMEGデータ
0024 % bexp    : 作成されたMEGデータ(bexp0+ノイズ)
0025 % J       : MEGデータに対応する電流分布
0026 % Iextract: Jの各行に対応する全脳モデル頂点インデックス
0027 %
0028 %
0029 % Copyright (C) 2011, ATR All Rights Reserved.
0030 % License : New BSD License(see VBMEG_LICENSE.txt)
0031 
0032 MEGinfo.nprsamp = MEGinfo.nprsamp;
0033 Noise = MEGinfo.Noise;
0034 Noise_file = MEGinfo.Noise_file;
0035 ix = MEGinfo.ix;
0036 ix0 = MEGinfo.ix;
0037 Source = MEGinfo.pattern;
0038 Rfilt = MEGinfo.Rfilt;
0039 Rmax = Rfilt*2;
0040 
0041 % 隣接点情報
0042 load([parm.dir.brain parm.file.brain],'nextDD','nextIX');
0043 
0044 % 電流強度が0より大きくなり得る頂点集合
0045 Iextract = [];
0046 for i = 1:length(ix)
0047   Iextract = [Iextract; ix{i}];
0048 end
0049 Iextract = unique(Iextract);
0050 
0051 tmp = []; % 電源からRmax以内の点
0052 for i = 1:length(Iextract)
0053   tmp_ix = find(nextDD{Iextract(i)}<=Rmax);
0054   tmp = [tmp; nextIX{Iextract(i)}(tmp_ix)];
0055 end
0056 Iextract = unique(tmp);
0057 I = length(Iextract);
0058 
0059 % 全脳インデックスから部分脳(Iextract)インデックスへの変換テーブル
0060 Itrans = zeros(length(nextIX),1);
0061 Itrans(Iextract) = 1:length(Iextract);
0062 
0063 % 平滑化フィルタ
0064 CW = vb_smooth_filter_norm(parm,Rfilt,Rmax,Iextract);
0065 
0066 % リードフィールド
0067 load([parm.dir.basis parm.file.basis],'basis');
0068 basis = basis(Iextract,:);
0069 
0070 % 電流パターンの作成
0071 J = zeros(I,MEGinfo.T);
0072 
0073 for i = 1:length(ix)
0074   si = Source{i};
0075   for j = 1:size(si,1)
0076     t1 = si(j,1);
0077     t2 = si(j,2);
0078     dt = t2-t1;
0079     JJ = si(j,3) * sin(pi*((0:dt)./dt)).^2;
0080     J(Itrans(ix{i}),t1:t2) = J(Itrans(ix{i}),t1:t2) + ...
0081     repmat(JJ,[length(ix{i}) 1]);
0082   end
0083 end
0084 
0085 % ノイズ無しMEG観測データ
0086 J = CW * J;
0087 bexp0 = basis' * J;
0088 
0089 % 観測ノイズを加える
0090 tmp = bexp0(:,MEGinfo.nprsamp+1:size(bexp0,2));
0091 tmp = mean(abs(tmp'));
0092 nlevel = MEGinfo.SN * max(tmp);
0093 
0094 bexp = bexp0;
0095 switch MEGinfo.Noise
0096  case 1,
0097   bexp = bexp + nlevel * randn(size(bexp));
0098  case {2,3},
0099   load(MEGinfo.Noise_file,'cv0');
0100   N = size(cv0,1);
0101   cv0 = cv0 ./ (sum(diag(cv0)) / N);
0102   if MEGinfo.Noise == 2
0103     cv0 = diag(diag(cv0));
0104   end
0105   R = chol(cv0);
0106   bexp = bexp + nlevel * R * randn(size(bexp));
0107 end
0108