vb_job_bad_trial_interact

 Find bad channel/trials by checking max value in each channel/trial
 --- Usage
 vb_job_bad_trial(datafile,job_mode, ...
    thred_val,prob_val,thred_ch,actfile,stat_file,std_mode)
 --- Input
 datafile: cell string for multiple MEG/EEG epoch data files
 job_mode: 
      1: channel/trial rejection in batch mode
     10: save trial statics without channel/trial rejection
 [actfile   '.info.mat']
  : file name to save bad channel/trial info for batch mode
    or save trial statics without channel/trial rejection
 [stat_file '.info.mat']
  : file name having trial statics made at the previous job
 
 ---- Threshold value for bad channel selection
 - thred_val = [9 8]; 
   threshold value for max(y) & diff(y) in one trial
   if (ymax/ystd) > thred_val(1) | difmax/difstd > thred_val(2), 
      it is marked as bad
 - prob_val  = 0.023;
   threshold probability to select bad channel
   if more than (prob_val*Ntrial) trials are bad, the channel is rejected
 - thred_ch  = 1/4;
   if more than (Nch*thred_ch) channels are bad, the trial is rejected,
   before checking bad channel.
 - std_mode  = 1;
   = 0: std is averaged over all trials
   = 1: std is averaged within one session (one BDF file)
 
 2009-2-1 Masa-aki Sato

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

0001 function  [fileinfo,flg] = vb_job_bad_trial_interact(datafile,job_mode, ...
0002     thred_val,prob_val,thred_ch,actfile,stat_file,std_mode)
0003 % Find bad channel/trials by checking max value in each channel/trial
0004 % --- Usage
0005 % vb_job_bad_trial(datafile,job_mode, ...
0006 %    thred_val,prob_val,thred_ch,actfile,stat_file,std_mode)
0007 % --- Input
0008 % datafile: cell string for multiple MEG/EEG epoch data files
0009 % job_mode:
0010 %      1: channel/trial rejection in batch mode
0011 %     10: save trial statics without channel/trial rejection
0012 % [actfile   '.info.mat']
0013 %  : file name to save bad channel/trial info for batch mode
0014 %    or save trial statics without channel/trial rejection
0015 % [stat_file '.info.mat']
0016 %  : file name having trial statics made at the previous job
0017 %
0018 % ---- Threshold value for bad channel selection
0019 % - thred_val = [9 8];
0020 %   threshold value for max(y) & diff(y) in one trial
0021 %   if (ymax/ystd) > thred_val(1) | difmax/difstd > thred_val(2),
0022 %      it is marked as bad
0023 % - prob_val  = 0.023;
0024 %   threshold probability to select bad channel
0025 %   if more than (prob_val*Ntrial) trials are bad, the channel is rejected
0026 % - thred_ch  = 1/4;
0027 %   if more than (Nch*thred_ch) channels are bad, the trial is rejected,
0028 %   before checking bad channel.
0029 % - std_mode  = 1;
0030 %   = 0: std is averaged over all trials
0031 %   = 1: std is averaged within one session (one BDF file)
0032 %
0033 % 2009-2-1 Masa-aki Sato
0034 %
0035 % Copyright (C) 2011, ATR All Rights Reserved.
0036 % License : New BSD License(see VBMEG_LICENSE.txt)
0037 
0038 if ~exist('std_mode','var') , std_mode = 1; end;
0039 
0040 if ~isempty(datafile) && ~iscell(datafile),
0041     tmp = datafile;
0042     datafile = {tmp};
0043 end;
0044 
0045 if exist('stat_file','var') && exist([stat_file '.info.mat'],'file')
0046     % Load statics info (ratio. ymax)
0047     load([stat_file '.info.mat'],'fileinfo')
0048 else
0049     % Get file info for multiple session MEG/EEG files
0050     % All trials are loaded
0051     fileinfo = vb_get_multi_fileinfo(datafile);
0052 end
0053 
0054 
0055 % ystd: amplitude mean calculated by histrgam for all trials in one channel
0056 % ymax : Nch x Ntrial
0057 % ystd : Nch x 1
0058 
0059 if isfield(fileinfo,'ymax1')
0060     ymax1 = fileinfo.ymax1;
0061     ymax2 = fileinfo.ymax2;
0062     ystd1 = fileinfo.ystd1;
0063     ystd2 = fileinfo.ystd2;
0064     if ~isempty(datafile) 
0065         fileinfo = vb_get_multi_fileinfo(datafile);
0066     end
0067 else
0068     ymax1 = [];
0069     ymax2 = [];
0070     ystd1 = [];
0071     ystd2 = [];
0072     fprintf('Load data ')
0073     for n=1:length(datafile)
0074         % Load MEG/EEG data
0075         fprintf('-')
0076         data_tmp = vb_load_meg_data(datafile{n});
0077         [max1,max2,std1,std2] = vb_channel_statics(data_tmp);
0078         ymax1 = [ymax1 , max1];
0079         ymax2 = [ymax2 , max2];
0080         ystd1 = [ystd1 , std1];
0081         ystd2 = [ystd2 , std2];
0082     end
0083     fprintf('\n')
0084 end
0085 
0086 if std_mode == 0 
0087     ystd1 = mean(ystd1,2);
0088     ystd2 = mean(ystd2,2);
0089 end
0090 
0091 ratio1 = ymax1;
0092 ratio2 = ymax2;
0093 
0094 ix1  = find(ystd1 > 0);
0095 ix2  = find(ystd2 > 0);
0096 
0097 if size(ystd1,2) == 1
0098     ratio1(ix1,:) = vb_repmultiply(ymax1(ix1,:), 1./ystd1(ix1));
0099     ratio2(ix2,:) = vb_repmultiply(ymax2(ix2,:), 1./ystd2(ix2));
0100 else
0101     ratio1(ix1) = ymax1(ix1)./ystd1(ix1);
0102     ratio2(ix2) = ymax2(ix2)./ystd2(ix2);
0103 end
0104 
0105 if ~exist('data','var')
0106     Nchannel = fileinfo.Nchannel;  % # of total channels
0107     Nsample  = fileinfo.Nsample ;  % # of samples
0108     Ntrial   = fileinfo.Ntrial  ;  % # of trials for each session
0109     Ntotal   = fileinfo.Ntotal  ;  % # of all trials
0110     % Subsampling step for EEG data storage
0111     tstep = 5;
0112     t = 1:tstep:Nsample;
0113     data = zeros(Nchannel,length(t),Ntotal);
0114     jtry = 1;
0115     fprintf('Load data ')
0116     for n=1:length(datafile)
0117         % Load EEG data
0118         fprintf('-')
0119         data_tmp = vb_load_meg_data(datafile{n});
0120         data(:,:,jtry:jtry+Ntrial(n)-1) = data_tmp(:,t,:);
0121         jtry = jtry+Ntrial(n);
0122     end
0123     [Nch,T,Ntry] = size(data);
0124     fprintf('[Nch,T,Ntry] = [%d, %d, %d] \n',Nch,T,Ntry)
0125 else
0126     data = [];
0127 end
0128 
0129 [Nch,Ntry] = size(ymax1);
0130 thred_num  = Nch * thred_ch;
0131 
0132 switch job_mode
0133 case    5
0134     % Check Bad channel interactively
0135     
0136     % 1. Find noisy trials with multiple bad channels
0137     ch_act  = [1:Nch];
0138     [try_err ,flg]= vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val,thred_num);
0139     try_act = vb_setdiff2([1:Ntry],try_err);
0140 
0141     % 2, Find Bad Channel with large max_ratio
0142     ch_bad = vb_find_bad_channel(ratio1(:,try_act),ratio2(:,try_act),...
0143                 thred_val,prob_val);
0144     
0145     % Check Bad channel interactively
0146     ch_bad = vb_plot_bad_channel(data,ch_bad,flg);
0147     
0148     % Find bad trials with large max_ratio after removing bad channel
0149     ch_act  = vb_setdiff2([1:Nch],ch_bad);
0150     try_bad = vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val);
0151     try_bad = vb_plot_bad_trial(data(ch_act,:,:),try_bad,flg(ch_act,:));
0152     try_act = vb_setdiff2([1:Ntry],try_bad);
0153 case    6
0154     % Check Bad trials interactively
0155     % 1. Find noisy trials with multiple bad channels
0156     ch_act  = [1:Nch];
0157     [try_err ,flg]= vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val,thred_num);
0158     
0159     % Check Bad trials interactively
0160     try_bad = vb_plot_bad_trial(data(ch_act,:,:),try_err,flg(ch_act,:));
0161     try_act = vb_setdiff2([1:Ntry],try_bad);
0162     
0163     % Find Bad Channel with large max_ratio
0164     ch_bad = vb_find_bad_channel(ratio1(:,try_act),ratio2(:,try_act),...
0165                 thred_val,prob_val);
0166     
0167     ch_act = vb_setdiff2([1:Nch],ch_bad);
0168 case    7
0169     % Check Bad trials interactively
0170     % 1. Find noisy trials with multiple bad channels
0171     ch_act  = [1:Nch];
0172     [try_err ,flg]= vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val,thred_num);
0173     try_act = vb_setdiff2([1:Ntry],try_err);
0174     
0175     % Find Bad Channel with large max_ratio
0176     ch_bad = vb_find_bad_channel(ratio1(:,try_act),ratio2(:,try_act),...
0177                 thred_val,prob_val);
0178     ch_act = vb_setdiff2([1:Nch],ch_bad);
0179     
0180     % Check Bad trials interactively
0181     try_bad = vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val);
0182     try_bad = vb_plot_bad_trial(data(ch_act,:,:),try_bad,flg(ch_act,:));
0183     try_act = vb_setdiff2([1:Ntry],try_bad);
0184 case    8
0185     % Check trials interactively
0186     
0187     % Find Bad Channel with large max_ratio
0188     [ch_bad,flg] = vb_find_bad_channel(ratio1,ratio2,thred_val,prob_val);
0189     
0190     % Find bad trials with large max_ratio after removing bad channel
0191     ch_act  = vb_setdiff2([1:Nch],ch_bad);
0192     try_bad = vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val);
0193 
0194     % Check Bad trials interactively
0195     try_bad = vb_plot_bad_trial(data(ch_act,:,:),try_bad,flg(ch_act,:));
0196     try_act = vb_setdiff2([1:Ntry],try_bad);
0197     try_err = [];
0198 end
0199 
0200 Ntry = length(try_act);
0201 Nact = length(ch_act);
0202 
0203 if job_mode==1 || job_mode==2
0204     fprintf('Active %d channel & %d epoch\n',Nact, Ntry)
0205     fprintf('Err trial= %d\n',length(try_err));
0206 end
0207 
0208 % channel & trial statistics
0209 fileinfo.ymax1 = ymax1;
0210 fileinfo.ymax2 = ymax2;
0211 fileinfo.ystd1 = ystd1;
0212 fileinfo.ystd2 = ystd2;
0213 
0214 fileinfo.channel_id = ch_act;
0215 fileinfo.act_trial  = try_act;
0216 fileinfo.err_trial  = try_err;
0217 
0218 if ~exist('actfile','var')||isempty(actfile), return;end
0219 
0220 vb_save_active_info(fileinfo, [actfile '.info.mat']);
0221 
0222 return
0223 %% ---- END ---- %%