vb_job_preprocess_meg3

---
 function vb_job_preprocess_meg3(proj_root, parm)

 This is exclusive use of SBI system.
 --
 Reject noisy sensors and trials, and divide trials according to
 experimental conditions in division file. If division file
 (divfile) was not specified, MEG data will not be divided. 
 If parm was not specified, recommended values will be used. 

 Input parameters (recommended value)
 - parm.rm_trial  : It can be used to directly specify trials to
                    be rejected. 
 - parm.verify    : If it is ON (constant), the script shows
                    trials that should be rejected and the user
                    must verify each trial by clicking 'OK' button
                    on a dialog (ON). 
 - parm.jobmode   : Select jobs and order of that in this
                    preprocessing. 
                    NOTICE that the result of preprocessing would
                    be changed if you specify the order of jobs.
                    The following (default) order is recomended.
                    1 -- Reject dead sensors 
                    2 -- Remove drift and correct baseline 
                    3 -- Remove trials in which EOG signal exceeds
                         a threshold.
                    4 -- Reject sensors on which MEG time courses
                         of some trials are noisy.  
                    5 -- Reject noisy trials
 - Ohter parameters for each job
 -- 1 -- Reject dead sensors 

 -- 2 -- Remove drift and correct baseline 
 - parm.freq      : MEG data frequency.
 - parm.driftremovalmode
                  : 'linear' -- remove drift with linear fitting.
                  : 'quadra' -- remove drift with quadra fitting.
                  : 'sin' -- remove drift with sin-cos fitting.
 - parm.basiocycle: when parm.driftrimovalmode = 'sin', you need
                    specify this parameter as vector of
                    basio-cycles, for example parm.basiocycle =
                    [4,8,16,32] (second). It is recomended that
                    basiocycle > trial time length.
 - parm.baseline_pretrig
                  : If this parameter is ON, baseline correction
                    by using median of MEG data at pretrigger
                    period for each sensor.

 -- 3 -- Remove trials in which EOG signal exceeds a threshold.
 - parm.num_eog_channels : the number of eog channels
 - parm.eog_channel_id : channel ID in which EOG data were stored
                         (for YOKOGAWA system)
 - parm.th_eog    : Threshold for EOG signal. This parameter must
                    be vector at the same length as the number of
                    'parm.num_eog_channels'. Each element of the
                    vector corresponds to the element of
                    'parm.eog_channel_id'.
                    ([0.75 0.75])

 -- 4,5 -- Reject sensors on which MEG time courses of some trials are noisy.
           Reject noisy trials.
 - parm.th_meg    : Threshold for noise in MEG data (8)
 - parm.meg_coef  : Threshold for power of MEG data considered as
                    a candidate of rejected sensor (5.0)
 - parm.Ntrial_th : The number of trials over which MEG time
                    courses exhibit noise on the same sensor (5)

 2004-12-28 Taku Yoshioka
 2005-08-16 Taku Yoshioka
 2006-06-14 Dai Kawawaki
 2006-12-19 Dai Kawawaki (bug fix)
 2007-04-17 Dai Kawawaki (bug fix)
---

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

0001 function vb_job_preprocess_meg3(proj_root, parm)
0002 %---
0003 % function vb_job_preprocess_meg3(proj_root, parm)
0004 %
0005 % This is exclusive use of SBI system.
0006 % --
0007 % Reject noisy sensors and trials, and divide trials according to
0008 % experimental conditions in division file. If division file
0009 % (divfile) was not specified, MEG data will not be divided.
0010 % If parm was not specified, recommended values will be used.
0011 %
0012 % Input parameters (recommended value)
0013 % - parm.rm_trial  : It can be used to directly specify trials to
0014 %                    be rejected.
0015 % - parm.verify    : If it is ON (constant), the script shows
0016 %                    trials that should be rejected and the user
0017 %                    must verify each trial by clicking 'OK' button
0018 %                    on a dialog (ON).
0019 % - parm.jobmode   : Select jobs and order of that in this
0020 %                    preprocessing.
0021 %                    NOTICE that the result of preprocessing would
0022 %                    be changed if you specify the order of jobs.
0023 %                    The following (default) order is recomended.
0024 %                    1 -- Reject dead sensors
0025 %                    2 -- Remove drift and correct baseline
0026 %                    3 -- Remove trials in which EOG signal exceeds
0027 %                         a threshold.
0028 %                    4 -- Reject sensors on which MEG time courses
0029 %                         of some trials are noisy.
0030 %                    5 -- Reject noisy trials
0031 % - Ohter parameters for each job
0032 % -- 1 -- Reject dead sensors
0033 %
0034 % -- 2 -- Remove drift and correct baseline
0035 % - parm.freq      : MEG data frequency.
0036 % - parm.driftremovalmode
0037 %                  : 'linear' -- remove drift with linear fitting.
0038 %                  : 'quadra' -- remove drift with quadra fitting.
0039 %                  : 'sin' -- remove drift with sin-cos fitting.
0040 % - parm.basiocycle: when parm.driftrimovalmode = 'sin', you need
0041 %                    specify this parameter as vector of
0042 %                    basio-cycles, for example parm.basiocycle =
0043 %                    [4,8,16,32] (second). It is recomended that
0044 %                    basiocycle > trial time length.
0045 % - parm.baseline_pretrig
0046 %                  : If this parameter is ON, baseline correction
0047 %                    by using median of MEG data at pretrigger
0048 %                    period for each sensor.
0049 %
0050 % -- 3 -- Remove trials in which EOG signal exceeds a threshold.
0051 % - parm.num_eog_channels : the number of eog channels
0052 % - parm.eog_channel_id : channel ID in which EOG data were stored
0053 %                         (for YOKOGAWA system)
0054 % - parm.th_eog    : Threshold for EOG signal. This parameter must
0055 %                    be vector at the same length as the number of
0056 %                    'parm.num_eog_channels'. Each element of the
0057 %                    vector corresponds to the element of
0058 %                    'parm.eog_channel_id'.
0059 %                    ([0.75 0.75])
0060 %
0061 % -- 4,5 -- Reject sensors on which MEG time courses of some trials are noisy.
0062 %           Reject noisy trials.
0063 % - parm.th_meg    : Threshold for noise in MEG data (8)
0064 % - parm.meg_coef  : Threshold for power of MEG data considered as
0065 %                    a candidate of rejected sensor (5.0)
0066 % - parm.Ntrial_th : The number of trials over which MEG time
0067 %                    courses exhibit noise on the same sensor (5)
0068 %
0069 % 2004-12-28 Taku Yoshioka
0070 % 2005-08-16 Taku Yoshioka
0071 % 2006-06-14 Dai Kawawaki
0072 % 2006-12-19 Dai Kawawaki (bug fix)
0073 % 2007-04-17 Dai Kawawaki (bug fix)
0074 %---
0075 %
0076 % Copyright (C) 2011, ATR All Rights Reserved.
0077 % License : New BSD License(see VBMEG_LICENSE.txt)
0078 
0079 proj_root = vb_rm_trailing_slash(proj_root);
0080 
0081 megfile = [proj_root filesep parm.megfile];
0082 divfile = parm.divfile;
0083 
0084 [MEGinfo, MEGdata, divNo] = init_preprocessing_meg(megfile,divfile);
0085 
0086 for job = 1:length(parm.jobmode)
0087   switch parm.jobmode(job)
0088    case 1 
0089     % Reject dead sensors
0090     [MEGinfo,MEGdata] = reject_dead_sensors(parm,MEGinfo,MEGdata);
0091     
0092    case 2
0093     % Remove drift and correct baseline
0094     [MEGdata] = remove_drift_correct_baseline(parm,MEGinfo,MEGdata);
0095     
0096    case 3
0097     % Remove trials in which EOG signal exceeds a threshold.
0098     [MEGinfo,MEGdata] = remove_trials_when_eye_moved(parm,MEGinfo,MEGdata);
0099    case 4
0100     % Reject sensors on which MEG time courses of some trials are noisy.
0101     [MEGinfo,MEGdata] = reject_noisy_sensors(parm,MEGinfo,MEGdata);
0102    case 5
0103     % Reject noisy trials
0104     [MEGinfo,MEGdata] = reject_noisy_trials(parm,MEGinfo,MEGdata);
0105    otherwise
0106     error('You cannot specify jobmode any other than [1:5].');
0107   end %switch(parm.jobmode)
0108 end %for job = 1:length(parm.jobmode)
0109 
0110 %save division
0111 save_division(parm, megfile, MEGinfo, MEGdata, divNo);
0112 return
0113 
0114 
0115 
0116 
0117 %internal functions
0118 
0119 %
0120 % Load MEG data
0121 %
0122 function [MEGinfo,MEGdata,divNo] = init_preprocessing_meg(megfile,divfile)
0123 fprintf('--- initializing... \n');
0124 
0125 MEGinfo = vb_load_meg_info(megfile);
0126 bexp = vb_load_meg_data(megfile); 
0127 % eog = load_eog_data(megfile);
0128 eog = vb_load_megeeg_data(megfile);
0129 [pick, Qpick, CoilWeight] = vb_load_sensor(megfile);
0130 Ntrial = MEGinfo.Nrepeat_org;
0131 
0132 % Load division file
0133 divNo = load_division_file(bexp,divfile,Ntrial);
0134 
0135 % set meg data
0136 MEGdata = set_meg_data(bexp,eog,pick,Qpick,CoilWeight);
0137 
0138 fprintf('--- done.\n');
0139 return
0140 
0141 %
0142 % Load division file
0143 %
0144 function divisionNumber = load_division_file(bexp,divfile,Ntrial)
0145 if isempty(divfile),
0146   divisionNumber = zeros(Ntrial,1);
0147 else
0148   divisionNumber = load('-ascii',divfile);
0149   if length(divisionNumber)~=Ntrial;
0150     error(['The number of trials does not match with the number of ' ...
0151        'experimental conditions.']);
0152   end
0153 end
0154 return
0155 
0156 %
0157 % Reject dead sensors
0158 %
0159 function [MEGinfo,MEGdata] = reject_dead_sensors(parm,MEGinfo,MEGdata)
0160 fprintf('--- searching dead sensosrs...\n');
0161 
0162 [bexp,eog,pick,Qpick,CoilWeight]= get_meg_data(MEGdata);
0163 chNo = MEGinfo.MEGch_name;%fixed by D.K. at 06-12-19
0164 trNo = MEGinfo.TrialNumber;
0165 sampleTime = MEGinfo.Nsample;
0166 Ntrial = MEGinfo.Nrepeat;
0167 Nchannel = MEGinfo.Nchannel;
0168 rejectChannel = [];
0169 
0170 h       = waitbar(0,'Check dead sensor');
0171 prg     = 0;
0172 prg_all = Nchannel;
0173 vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0174 refresh(h);
0175 if parm.verify==ON, hh = figure; hold on; end;
0176 
0177 for ixch = 1:Nchannel
0178   waitbar(ixch/Nchannel);
0179   
0180   for ixtr=1:Ntrial
0181     max_value = max(bexp(ixch,:,ixtr));
0182     min_value = min(bexp(ixch,:,ixtr));
0183     if max_value==min_value
0184       if parm.verify==ON, 
0185     figure(hh);clf(hh);hold on; 
0186     plot(bexp(:,:,ixtr)',':');
0187     plot(bexp(ixch,:,ixtr)','linewidth',4);
0188     title(sprintf('trial %d - sensor %d',trNo(ixtr),chNo(ixch)));
0189     xlim([1 sampleTime]);
0190     button = ...
0191         questdlg(['Reject sensor ' num2str(chNo(ixch)) '?'],...
0192              'Sensor rejection','Yes','No','Yes');
0193       else
0194     button = 'Yes';
0195       end
0196       if strcmp(button,'Yes')
0197     disp(sprintf('---- sensor %d was rejected at trial %d.',...
0198              chNo(ixch), trNo(ixtr)));
0199     rejectChannel = [rejectChannel; chNo(ixch)];
0200     break;
0201       end
0202     end
0203   end
0204   
0205   for ii=1:15; vb_disp_nonl(sprintf('\b')); end
0206   vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0207   prg = prg+1;
0208 end
0209 
0210 close(h); 
0211 vb_disp_nonl(sprintf('\n'));
0212 
0213 if parm.verify==ON, close(hh); end;
0214 % modified by DK at 070417
0215 if isempty(vb_setdiff2(MEGinfo.MEGch_id, rejectChannel));
0216   error('All sensors were rejected. Please check parameters.');
0217 end
0218 [MEGinfo,bexp,pick,Qpick,CoilWeight] = ...
0219     vb_channel_rejection(MEGinfo,bexp,pick,Qpick,rejectChannel,CoilWeight);
0220 MEGdata = set_meg_data(bexp,eog,pick,Qpick,CoilWeight);
0221 fprintf('--- done.\n');
0222 return
0223 
0224 
0225 %
0226 % remove drift & baseline correction
0227 %
0228 function [MEGdata] = remove_drift_correct_baseline(parm,MEGinfo,MEGdata)
0229 fprintf('--- removing drift component... \n');
0230 
0231 [bexp,eog,pick,Qpick,CoilWeight]= get_meg_data(MEGdata);
0232 chNo = MEGinfo.MEGch_name;%fixed by D.K. at 06-12-19
0233 trNo = MEGinfo.TrialNumber;
0234 Ntrial = MEGinfo.Nrepeat;
0235 Nchannel = MEGinfo.Nchannel;
0236 sampleTime = MEGinfo.Nsample;
0237 time = [1:1000/parm.freq:sampleTime];
0238 pre_trigger = MEGinfo.Pretriger;
0239 
0240 h       = waitbar(0,'remove drift & correct baseline');
0241 prg     = 0;
0242 prg_all = Ntrial;
0243 vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0244 refresh(h);
0245 
0246 for ixtr=1:Ntrial
0247   waitbar(ixtr/Ntrial);
0248 
0249   %remove drift
0250   switch lower(parm.driftremovalmode)
0251    case 'sin'
0252     bexp(:,:,ixtr) = vb_polysincosfilt(bexp(:,:,ixtr),...
0253                     parm.basiocycle,...
0254                     parm.freq);
0255 
0256    case 'linear'
0257     for ixch=1:Nchannel
0258       P = polyfit(time,bexp(ixch,:,ixtr),1);
0259       dr = P(1)*time + P(2);
0260       bexp(ixch,:,ixtr) = bexp(ixch,:,ixtr) - dr;
0261     end
0262     
0263    case 'quadra'
0264     for ixch=1:Nchannel
0265       P = polyfit(time,bexp(ixch,:,ixtr),2);
0266       dr = P(1)*time.^2 + P(2)*time + P(3);
0267       bexp(ixch,:,ixtr) = bexp(ixch,:,ixtr) - dr;
0268     end
0269 
0270    otherwise %D.K. 06-12-19
0271     error('invalid usage of parameter for drift rimoval.');
0272     
0273   end
0274   
0275   %baselien correction with pretrigger
0276   if parm.baseline_pretrig == ON
0277     base_line = median(bexp(:,1:pre_trigger,ixtr),2);
0278     bexp(:,:,ixtr) = bexp(:,:,ixtr) - repmat(base_line, 1, sampleTime);
0279   end
0280   
0281   for ii=1:15; vb_disp_nonl(sprintf('\b')); end
0282   vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0283   prg = prg+1;
0284 end
0285 close(h);
0286 vb_disp_nonl(sprintf('\n'));
0287 
0288 MEGdata = set_meg_data(bexp,eog,pick,Qpick,CoilWeight);
0289 
0290 fprintf('--- done.\n');
0291 return
0292 
0293 
0294 %
0295 % Remove trials in which EOG signal exceeds a threshold.
0296 %
0297 function [MEGinfo,MEGdata] = remove_trials_when_eye_moved(parm,MEGinfo,MEGdata);
0298 fprintf('--- rejecting trials in which subject''s eyes were moved...\n');
0299 
0300 [bexp,eog,pick,Qpick,CoilWeight]= get_meg_data(MEGdata);
0301 chNo = MEGinfo.MEGch_name;%fixed by D.K. at 06-12-19
0302 trNo = MEGinfo.TrialNumber;
0303 Ntrial = MEGinfo.Nrepeat;
0304 sampleTime = MEGinfo.Nsample;
0305 
0306 Neog = parm.num_eog_channels;
0307 eog_id = parm.eog_channel_id;
0308 th_eog = parm.th_eog;
0309 ix_eog = chname2index_eog(MEGinfo,Neog,eog_id);
0310 
0311 rejectTrial = [];
0312 
0313 h       = waitbar(0,'remove drift & correct baseline');
0314 prg     = 0;
0315 prg_all = Ntrial;
0316 vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0317 refresh(h);
0318 if parm.verify==ON, hh = figure; hold on; end;
0319 
0320 if ~isempty(eog) 
0321   for ixtr=1:Ntrial
0322     waitbar(ixtr/Ntrial); 
0323     
0324     for n = 1:Neog
0325       diff_eog = max(eog(ix_eog(n),:,ixtr))-min(eog(ix_eog(n),:,ixtr));
0326       if diff_eog > th_eog(n) %fixed 06-12-19 D.K.
0327     if parm.verify==ON, 
0328       figure(hh);clf(hh);hold on;
0329       subplot(2,1,1); 
0330       plot(eog(ix_eog,:,ixtr)');
0331       legend('eog channel 1','eog channel 2','...'); 
0332       xlim([1 sampleTime]);
0333       subplot(2,1,2); 
0334       plot(bexp(:,:,ixtr)'); 
0335       xlim([1 sampleTime]);
0336       title(['Trial ' num2str(ixtr)]);
0337       button = questdlg(['Reject trial ' num2str(trNo(ixtr)) '?'],...
0338                 '   EOG rejection','Yes','No','Yes');
0339     else
0340       button = 'Yes';
0341     end 
0342     
0343     if strcmp(button, 'Yes')
0344       rejectTrial = [rejectTrial; trNo(ixtr)];
0345       fprintf('---- trial %d was rejected.\n',trNo(ixtr));
0346       break;
0347     end
0348       end 
0349     end %for n = 1:Neog
0350     
0351     for ii=1:15; vb_disp_nonl(sprintf('\b')); end
0352     vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0353     prg = prg+1;
0354   end %for ixtr=1:Ntrial
0355 end 
0356 
0357 close(h);
0358 vb_disp_nonl(sprintf('\n'));
0359 
0360 if parm.verify==ON, close(hh); end; 
0361 % modified by DK at 070417
0362 if isempty(vb_setdiff2(MEGinfo.TrialNumber, rejectTrial)) 
0363   error('All trials were rejected. Please check parameters.');
0364 end
0365 [MEGinfo,bexp,eog] = vb_trial_rejection(MEGinfo,bexp,eog,rejectTrial);
0366 MEGdata = set_meg_data(bexp,eog,pick,Qpick,CoilWeight);
0367 fprintf('--- done.\n');
0368 return
0369 
0370 function ix_eog = chname2index_eog(MEGinfo,Neog,eog_id)
0371 switch (lower(MEGinfo.device))
0372  case 'sbi'
0373   ix_eog = 1:Neog;
0374  case 'yokogawa'
0375   chname = MEGinfo.EEGinfo.all_ch_name;
0376   [tmp,ix_eog] = intersect(chname,eog_id);
0377   if isempty(tmp)
0378     error('There are no such channel name.');
0379   end
0380   if length(tmp) ~= length(eog_id)
0381     warning('Some channel names you specified did not match.');
0382   end
0383  otherwise
0384   error('MEG system is unknown. Please check the MEG system you used.');
0385 end
0386 return
0387 
0388 
0389 %
0390 % Reject sensors on which MEG time courses of some trials are
0391 % noisy.
0392 %
0393 function [MEGinfo,MEGdata] = reject_noisy_sensors(parm,MEGinfo,MEGdata)
0394 fprintf('--- rejecting noisy sensors...\n');
0395 [bexp,eog,pick,Qpick,CoilWeight]= get_meg_data(MEGdata);
0396 %fixed by D.K. at 06-12-19
0397 %fixed by D.K. at 07-04-17
0398 chNo = MEGinfo.MEGch_id;
0399 trNo = MEGinfo.TrialNumber;
0400 Nchannel = MEGinfo.Nchannel;
0401 Ntrial = MEGinfo.Nrepeat;
0402 sampleTime = MEGinfo.Nsample;
0403 
0404 rejectChannel = [];
0405 
0406 ixcount = cell(Nchannel,1);
0407 ratio_all = [];
0408 
0409 h       = waitbar(0,'Check noisy sensor');
0410 prg     = 0;
0411 prg_all = Ntrial;
0412 vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0413 refresh(h);
0414 
0415 % count the number of noisy trials in each of sensors
0416 for ixtr=1:Ntrial
0417   waitbar(ixtr/Ntrial); 
0418   
0419   tmp = abs(bexp(:,:,ixtr));  
0420   Z0 = median(tmp(:));
0421   for ixch=1:Nchannel
0422     ratio = max(tmp(ixch,:))/median(tmp(ixch,:));
0423     ratio_all = [ratio_all; ratio];
0424     if max(tmp(ixch,:)) > parm.meg_coef*Z0 & ratio > parm.th_meg
0425       ixcount{ixch} = [ixcount{ixch}; ixtr];
0426     end
0427   end
0428   
0429   for ii=1:15; vb_disp_nonl(sprintf('\b')); end
0430   vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0431   prg = prg+1;
0432 end
0433 
0434 clear tmp
0435 close(h); 
0436 vb_disp_nonl(sprintf('\n'));
0437 
0438 % confirm sensor rejection
0439 if parm.verify==ON, hh = figure; hold on; end
0440 ntrial_th = parm.Ntrial_th;
0441 for n=1:Nchannel
0442   if length(ixcount{n})>=ntrial_th % fixed 06-12-19 D.K.
0443     % Calculate ratio of max to median for all trials
0444     r = [];
0445     for c=1:length(ixcount{n})
0446       tmp = abs(bexp(n,:,ixcount{n}(c)));
0447       r = [r; max(tmp)/median(tmp)];
0448     end
0449     [r,ii] = sort(r); r = flipud(r); ii = flipud(ii); %fixed at 06-12-19 by D.K.
0450     
0451     if parm.verify==ON, 
0452       figure(hh);clf(hh);hold on;
0453       for trcount=1:ntrial_th
0454     subplot(ntrial_th,1,trcount); hold on;
0455     plot(bexp(:,:,ixcount{n}(ii(end-ntrial_th+trcount)))',':');
0456     xlim([1 sampleTime]);
0457     plot(bexp(n,:,ixcount{n}(ii(end-ntrial_th+trcount)))','LineWidth',4);
0458     hold off;
0459     title(['Trial ' num2str(trNo(ixcount{n}(ii(end-ntrial_th+trcount)))) ...
0460            ' max/med:' sprintf('%2.1f',r(end-ntrial_th+trcount))]);
0461       end
0462       button = questdlg(['Reject sensor ' num2str(chNo(n)) '?'],...
0463             'Sensor rejection','Yes','No','Yes');
0464     else
0465       button = 'Yes'; 
0466     end
0467     
0468     if strcmp(button,'Yes')
0469       rejectChannel = [rejectChannel; chNo(n)];
0470       fprintf('---- channel %d was rejected.\n',chNo(n));
0471     end
0472   end
0473 end
0474 rejectChannel = sort(rejectChannel);
0475 
0476 if parm.verify==ON, close(hh); end; 
0477 clear tmp;
0478 
0479 % modified by DK at 070417
0480 if isempty(vb_setdiff2(MEGinfo.MEGch_id, rejectChannel));
0481   error('All sensors were rejected. Please check parameters.');
0482 end
0483 [MEGinfo,bexp,pick,Qpick,CoilWeight] = ...
0484     vb_channel_rejection(MEGinfo,bexp,pick,Qpick,rejectChannel,CoilWeight);
0485 MEGdata = set_meg_data(bexp,eog,pick,Qpick,CoilWeight);
0486 fprintf('--- done.\n');
0487 return
0488 
0489 
0490 %
0491 % Reject noisy trials
0492 %
0493 function [MEGinfo,MEGdata] = reject_noisy_trials(parm,MEGinfo,MEGdata)
0494 fprintf('--- rejecting noisy trials...\n');
0495 [bexp,eog,pick,Qpick,CoilWeight]= get_meg_data(MEGdata);
0496 %fixed by D.K. at 06-12-19
0497 %fixed by D.K. at 07-04-17
0498 chNo = MEGinfo.MEGch_id;
0499 trNo = MEGinfo.TrialNumber;
0500 Nchannel = MEGinfo.Nchannel;
0501 Ntrial = MEGinfo.Nrepeat;
0502 sampleTime = MEGinfo.Nsample;
0503 
0504 rejectTrial = [];
0505 
0506 ratio_max = zeros(Ntrial,1);
0507 %ix = setdiff(1:Nsensor,ix_sensor_rej);
0508 h       = waitbar(0,'Check noisy trial'); 
0509 prg     = 0;
0510 prg_all = Ntrial;
0511 vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0512 
0513 for ixtr = 1:Ntrial
0514   waitbar(ixtr/Ntrial); 
0515   
0516   tmp = abs(bexp(:,:,ixtr));
0517   Z0 = median(tmp(:));
0518   for ixch=1:Nchannel
0519     ratio = max(tmp(ixch,:))/median(tmp(ixch,:));
0520     if max(tmp(ixch,:)) > parm.meg_coef*Z0 
0521       ratio_max(ixtr) = max([ratio_max(ixtr) ratio]);
0522     end
0523   end
0524   
0525   for ii=1:15; vb_disp_nonl(sprintf('\b')); end
0526   vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0527   prg = prg+1;
0528 end
0529 
0530 clear tmp;
0531 close(h); 
0532 vb_disp_nonl(sprintf('\n'));
0533 
0534 if parm.verify==ON, hh=figure; end;
0535 for ixtr = 1:Ntrial
0536   if ratio_max(ixtr)>parm.th_meg | ~isempty(find(parm.rm_trial==trNo(ixtr))), 
0537     if parm.verify==ON,
0538       figure(hh);clf(hh);hold on;
0539       plot(bexp(:,:,ixtr)'); 
0540       xlim([1 sampleTime]);
0541       title(['Trial ' num2str(trNo(ixtr)) ...
0542          sprintf(' max/med:%2.1f',ratio_max(ixtr))]);
0543       button = questdlg(['Reject trial ' num2str(trNo(ixtr)) '?'],...
0544             'Trial rejection','Yes','No','Yes');
0545     else
0546       button = 'Yes';
0547     end
0548     
0549     if strcmp(button,'Yes')
0550       rejectTrial = [rejectTrial; trNo(ixtr)];
0551       fprintf('---- trial %d was rejected.\n',trNo(ixtr));
0552     end
0553   end
0554 end
0555 if parm.verify == ON, close(hh); end;
0556 
0557 % modified by DK at 070417
0558 if isempty(vb_setdiff2(MEGinfo.TrialNumber, rejectTrial)) 
0559   error('All trials were rejected. Please check parameters.');
0560 end
0561 [MEGinfo,bexp,eog] = vb_trial_rejection(MEGinfo,bexp,eog,rejectTrial);
0562 MEGdata = set_meg_data(bexp,eog,pick,Qpick,CoilWeight);
0563 fprintf('--- done.\n');
0564 return
0565 
0566 
0567 %
0568 % Trial division
0569 %
0570 function save_division(parm, megfile, MEGinfo, MEGdata, divNo)
0571 fprintf('--- now saving...\n');
0572 [bexp_org,eog_org,pick_org,Qpick,CoilWeight]=get_meg_data(MEGdata);
0573 chNo = MEGinfo.MEGch_name;%fixed by D.K. at 06-12-19
0574 trNo = MEGinfo.TrialNumber;
0575 Nchannel = MEGinfo.Nchannel;
0576 Ntrial = MEGinfo.Nrepeat;
0577 MEGinfo.TrialNumber_org = trNo;
0578 
0579 %reject rejected trials in divNo
0580 divNo = divNo(trNo);
0581 
0582 division = unique(divNo); 
0583 Ncond = length(division);
0584 disp(['---- ' num2str(length(division)) ' conditions in division file']);
0585 
0586 for div=1:length(division)
0587   ixtr = find(divNo==division(div));
0588 
0589   %divide trials
0590   bexp = bexp_org(:,:,ixtr);
0591   if ~isempty(eog_org),eog = eog_org(:,:,ixtr);else eog=[]; end;
0592   
0593   MEGinfo.Nrepeat = length(ixtr);
0594   MEGinfo.TrialNumber = trNo(ixtr);
0595 
0596   disp(['---- ' num2str(length(ixtr)) ' trials were remained ' ...
0597     'in condition ' num2str(division(div))]);
0598 
0599   if parm.verify==ON, 
0600     answer = inputdlg(['Name for condition ' num2str(division(div)) ...
0601                '?'],['Condition ' num2str(division(div))]);
0602     if isempty(answer) % if "cancel" button was pushed.
0603       error('Process was aborted.');
0604     end
0605   else
0606     answer{1} = ['cond' num2str(division(div))];
0607   end
0608   
0609   filename = strrep(megfile,'.meg.mat',['_' answer{1} '.meg.mat']);
0610   disp(['    saved as ' filename]);
0611   
0612   % Save MEG data
0613   switch lower(MEGinfo.device)
0614    case 'sbi'
0615     pick = cat( 2, pick_org(1:Nchannel,:), pick_org([Nchannel+1:Nchannel*2],:)  );
0616     pick = cat( 2, pick,                   Qpick(1:Nchannel,:) );
0617     vb_fsave(filename,'pick','bexp','MEGinfo','eog'); 
0618    case 'yokogawa'
0619     eeg = eog;
0620     pick = pick_org;
0621     vb_fsave(filename,'pick','Qpick','bexp','MEGinfo','eeg'); 
0622    otherwise
0623     pick = cat( 2, pick_org(1:Nchannel,:), pick_org([Nchannel+1:Nchannel*2],:)  );
0624     pick = cat( 2, pick,                   Qpick(1:Nchannel,:) );
0625     vb_fsave(filename,'pick','bexp','MEGinfo','eog'); 
0626   end
0627 fprintf('--- done.\n');
0628 end
0629 
0630 function MEGdata = set_meg_data(bexp,eog,pick,Qpick,CoilWeight)
0631 MEGdata.bexp = bexp;
0632 MEGdata.eog = eog;
0633 MEGdata.pick = pick;
0634 MEGdata.Qpick = Qpick;
0635 MEGdata.CoilWeight = CoilWeight;
0636 return
0637 
0638 function [bexp,eog,pick,Qpick,CoilWeight]=get_meg_data(MEGdata)
0639 bexp = MEGdata.bexp;
0640 eog = MEGdata.eog;
0641 pick = MEGdata.pick;
0642 Qpick = MEGdata.Qpick;
0643 CoilWeight = MEGdata.CoilWeight;
0644 return