vb_job_bad_trial

 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
    2: 2D-Plot bad trial & channel for set of thresholds
    3-6: Interactive mode
    7,8: Stat mode to view distribution of trials with set of
      10: just save trial statics without channel/trial rejection

 [savefile   '.info.mat']
  : file name to save bad channel/trial info
    or save trial statics without channel/trial rejection (mode 7,8,10)
 [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
   For stat mode only, you can select multi values as list like
   thred_val = [8:11; 7:10];
               [vals_for_max(y); vals_for_diff(y)];

 - prob_val  = 0.023;
   threshold probability to select bad channel
   if more than (prob_val*Ntrial) trials are bad, the channel is rejected
   For stat mode only, you can select multi values as list like
   prob_val  = [0.05   0.01   0.001];

 - 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)
 
 2022-4-27 k_suzuki
 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,fig_handler] = vb_job_bad_trial(datafile,job_mode, ...
0002     thred_val,prob_val,thred_ch,savefile,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 %    2: 2D-Plot bad trial & channel for set of thresholds
0012 %    3-6: Interactive mode
0013 %    7,8: Stat mode to view distribution of trials with set of
0014 %      10: just save trial statics without channel/trial rejection
0015 %
0016 % [savefile   '.info.mat']
0017 %  : file name to save bad channel/trial info
0018 %    or save trial statics without channel/trial rejection (mode 7,8,10)
0019 % [stat_file '.info.mat']
0020 %  : file name having trial statics made at the previous job
0021 %
0022 % ---- Threshold value for bad channel selection
0023 % - thred_val = [9 8];
0024 %   threshold value for max(y) & diff(y) in one trial
0025 %   if (ymax/ystd) > thred_val(1) | difmax/difstd > thred_val(2),
0026 %      it is marked as bad
0027 %   For stat mode only, you can select multi values as list like
0028 %   thred_val = [8:11; 7:10];
0029 %               [vals_for_max(y); vals_for_diff(y)];
0030 %
0031 % - prob_val  = 0.023;
0032 %   threshold probability to select bad channel
0033 %   if more than (prob_val*Ntrial) trials are bad, the channel is rejected
0034 %   For stat mode only, you can select multi values as list like
0035 %   prob_val  = [0.05   0.01   0.001];
0036 %
0037 % - thred_ch  = 1/4;
0038 %   if more than (Nch*thred_ch) channels are bad, the trial is rejected,
0039 %   before checking bad channel.
0040 %
0041 % - std_mode  = 1;
0042 %   = 0: std is averaged over all trials
0043 %   = 1: std is averaged within one session (one BDF file)
0044 %
0045 % 2022-4-27 k_suzuki
0046 % 2009-2-1 Masa-aki Sato
0047 %
0048 % Copyright (C) 2011, ATR All Rights Reserved.
0049 % License : New BSD License(see VBMEG_LICENSE.txt)
0050 
0051 if ~exist('std_mode','var') , std_mode = 1; end;
0052 
0053 if ~isempty(datafile) && ~iscell(datafile),
0054     tmp = datafile;
0055     datafile = {tmp};
0056 end;
0057 
0058 if exist('stat_file','var') && exist([stat_file '.info.mat'],'file')
0059     % Load statics info (ratio. ymax)
0060     load([stat_file '.info.mat'],'fileinfo')
0061 else
0062     % Get file info for multiple session MEG/EEG files
0063     % All trials are loaded
0064     if  length(datafile)==1&&strcmp(datafile{1}(end-8:end), '.info.mat')
0065         % If datafile is info.mat, then unpack all session files
0066         info = load(datafile{1});
0067         datafile = info.fileinfo.filename;
0068     end
0069     fileinfo = vb_get_multi_fileinfo(datafile);
0070 end
0071 
0072 % for compatibility
0073 if isfield(fileinfo, 'channel_id')
0074     data_format = 'old';
0075 elseif isfield(fileinfo, 'ActiveChannel')
0076     data_format = 'new';
0077 end
0078 if strcmp(data_format, 'new')
0079     fileinfo = vb_fileinfo_active_field_convert_to('old', fileinfo, datafile{1});
0080 end
0081 
0082 % ystd: amplitude mean calculated by histrgam for all trials in one channel
0083 % ymax : Nch x Ntrial
0084 % ystd : Nch x 1
0085 
0086 if isfield(fileinfo,'ymax1')
0087     ymax1 = fileinfo.ymax1;
0088     ymax2 = fileinfo.ymax2;
0089     ystd1 = fileinfo.ystd1;
0090     ystd2 = fileinfo.ystd2;
0091     if ~isempty(datafile) 
0092         fileinfo = vb_get_multi_fileinfo(datafile);
0093     end
0094 else
0095     ymax1 = [];
0096     ymax2 = [];
0097     ystd1 = [];
0098     ystd2 = [];
0099     fprintf('Load data ')
0100     for n=1:length(datafile)
0101         % Load MEG/EEG data
0102         fprintf('-')
0103         data = vb_load_meg_data(datafile{n});
0104         [max1,max2,std1,std2] = vb_channel_statics(data);
0105         ymax1 = [ymax1 , max1];
0106         ymax2 = [ymax2 , max2];
0107         ystd1 = [ystd1 , std1];
0108         ystd2 = [ystd2 , std2];
0109     end
0110     fprintf('\n')
0111 end
0112 
0113 % Compile data among multi-files for interactive mode
0114 if job_mode==3||job_mode==4||job_mode==5||job_mode==6
0115     Nchannel = fileinfo.Nchannel;  % # of total channels
0116     Nsample  = fileinfo.Nsample ;  % # of samples
0117     Ntrial   = fileinfo.Ntrial  ;  % # of trials for each session
0118     Ntotal   = fileinfo.Ntotal  ;  % # of all trials
0119     % Subsampling step (sampled data are used for visualization)
0120     tstep = 1;
0121     t = 1:tstep:Nsample;
0122     data_plt = zeros(Nchannel,length(t),Ntotal);
0123     jtry = 1;
0124     fprintf('Load data ')
0125     for n=1:length(datafile)
0126         % Load EEG data
0127         fprintf('-')
0128         data_tmp = vb_load_meg_data(datafile{n});
0129         data_plt(:,:,jtry:jtry+Ntrial(n)-1) = data_tmp(:,t,:);
0130         jtry = jtry+Ntrial(n);
0131     end
0132     [Nch,T,Ntry] = size(data_plt);
0133     fprintf('[Nch,T,Ntry] = [%d, %d, %d] \n',Nch,T,Ntry)
0134 end
0135 
0136 if std_mode == 0 
0137     ystd1 = mean(ystd1,2);
0138     ystd2 = mean(ystd2,2);
0139 end
0140 
0141 ratio1 = ymax1;
0142 ratio2 = ymax2;
0143 
0144 ix1  = find(ystd1 > 0);
0145 ix2  = find(ystd2 > 0);
0146 
0147 if size(ystd1,2) == 1
0148     ratio1(ix1,:) = vb_repmultiply(ymax1(ix1,:), 1./ystd1(ix1));
0149     ratio2(ix2,:) = vb_repmultiply(ymax2(ix2,:), 1./ystd2(ix2));
0150 else
0151     ratio1(ix1) = ymax1(ix1)./ystd1(ix1);
0152     ratio2(ix2) = ymax2(ix2)./ystd2(ix2);
0153 end
0154 
0155 [Nch,Ntry] = size(ymax1);
0156 thred_num  = Nch * thred_ch;
0157 
0158 switch job_mode
0159 case    1
0160     % Remove error trials first (batch mode)
0161 
0162     % 1. Find noisy trials with multiple bad channels
0163     ch_act  = [1:Nch];
0164     [try_err ,flg]= vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val,thred_num);
0165     try_act = vb_setdiff2([1:Ntry],try_err);
0166 
0167     % 2. Find Bad Channel after removing noisy trials
0168     [ch_bad,flg] = vb_find_bad_channel(ratio1(:,try_act),ratio2(:,try_act), ...
0169                                     thred_val,prob_val);
0170     ch_act = vb_setdiff2([1:Nch],ch_bad);
0171     
0172     % 3. Find bad trials after removing bad channel
0173     try_bad  = vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val);
0174     try_act = vb_setdiff2([1:Ntry],try_bad);
0175 case    2
0176     %  2D-Plot bad trial & channel for set of thresholds
0177     vb_plot_max_ratio_image(ratio1,thred_val(1))
0178     vb_plot_max_ratio_image(ratio2,thred_val(2))
0179     
0180     ch_act  = 1:Nch;
0181     try_act = 1:Ntry;
0182     
0183 %---- 3-6 are interactive mode
0184 case    3
0185     % Check Bad channel interactively
0186     
0187     % 1. Find noisy trials with multiple bad channels
0188     ch_act  = [1:Nch];
0189     [try_err ,flg]= vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val,thred_num);
0190     try_act = vb_setdiff2([1:Ntry],try_err);
0191 
0192     % 2, Find Bad Channel with large max_ratio
0193     ch_bad = vb_find_bad_channel(ratio1(:,try_act),ratio2(:,try_act),...
0194                 thred_val,prob_val);
0195     
0196     % Check Bad channel interactively
0197     ch_bad = vb_plot_bad_channel(data_plt,ch_bad,flg);
0198     
0199     % Find bad trials with large max_ratio after removing bad channel
0200     ch_act  = vb_setdiff2([1:Nch],ch_bad);
0201     try_bad = vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val);
0202     try_bad = vb_plot_bad_trial(data_plt(ch_act,:,:),try_bad,flg(ch_act,:));
0203     try_act = vb_setdiff2([1:Ntry],try_bad);
0204 case    4
0205     % Check Bad trials interactively
0206     % 1. Find noisy trials with multiple bad channels
0207     ch_act  = [1:Nch];
0208     [try_err ,flg]= vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val,thred_num);
0209     
0210     % Check Bad trials interactively
0211     try_bad = vb_plot_bad_trial(data_plt(ch_act,:,:),try_err,flg(ch_act,:));
0212     try_act = vb_setdiff2([1:Ntry],try_bad);
0213     
0214     % Find Bad Channel with large max_ratio
0215     ch_bad = vb_find_bad_channel(ratio1(:,try_act),ratio2(:,try_act),...
0216                 thred_val,prob_val);
0217     
0218     ch_act = vb_setdiff2([1:Nch],ch_bad);
0219 case    5
0220     % Check Bad trials interactively
0221     % 1. Find noisy trials with multiple bad channels
0222     ch_act  = [1:Nch];
0223     [try_err ,flg]= vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val,thred_num);
0224     try_act = vb_setdiff2([1:Ntry],try_err);
0225     
0226     % Find Bad Channel with large max_ratio
0227     ch_bad = vb_find_bad_channel(ratio1(:,try_act),ratio2(:,try_act),...
0228                 thred_val,prob_val);
0229     ch_act = vb_setdiff2([1:Nch],ch_bad);
0230     
0231     % Check Bad trials interactively
0232     try_bad = vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val);
0233     try_bad = vb_plot_bad_trial(data_plt(ch_act,:,:),try_bad,flg(ch_act,:));
0234     try_act = vb_setdiff2([1:Ntry],try_bad);
0235 case    6
0236     % Check trials interactively
0237     % Find Bad Channel with large max_ratio
0238     [ch_bad,flg] = vb_find_bad_channel(ratio1,ratio2,thred_val,prob_val);
0239     
0240     % Find bad trials with large max_ratio after removing bad channel
0241     ch_act  = vb_setdiff2([1:Nch],ch_bad);
0242     try_bad = vb_find_bad_trial(ratio1,ratio2,ch_act,thred_val);
0243 
0244     % Check Bad trials interactively
0245     try_bad = vb_plot_bad_trial(data_plt(ch_act,:,:),try_bad,flg(ch_act,:));
0246     try_act = vb_setdiff2([1:Ntry],try_bad);
0247     try_err = [];
0248     
0249 %---- 7 and 8 are stat mode
0250 case    7
0251     % Plot ymax histgram
0252     fig_handler(1)=figure;
0253     subplot 121
0254     [ylist1,h1,y1] = inner_plot_stat_hist(ratio1 ,[]);
0255     hold on
0256     title('Max/std of amplitude')
0257     subplot 122
0258     [ylist2,h2,y2] = inner_plot_stat_hist(ratio2 ,[]);
0259     hold on
0260     title('Max/std of derivative')
0261     % set output variable (No change)
0262     %flg  = [ylist1; ylist2];
0263     %data = [h1(:), y1(:), h2(:), y2(:)];
0264     ch_act  = 1:Nch;
0265     try_act = 1:Ntry;
0266     try_err = [];
0267     flg = [];
0268 case    8
0269     if size(thred_val,1)==2
0270         thred_val_lst1 = thred_val(1,:);
0271         thred_val_lst2 = thred_val(2,:);
0272     else
0273         thred_val_lst1 = thred_val(1);
0274         thred_val_lst2 = thred_val(2);
0275     end
0276     % Plot ymax histgram
0277     fig_handler(1)=figure;
0278     subplot 121
0279     [ylist1,h1,y1] = inner_plot_stat_hist(ratio1 ,[]);
0280     hold on
0281     title('Max/std of amplitude')
0282     subplot 122
0283     [ylist2,h2,y2] = inner_plot_stat_hist(ratio2 ,[]);
0284     hold on
0285     title('Max/std of derivative')
0286     % set output variable (No change)
0287     %flg  = [ylist1; ylist2];
0288     %data = [h1(:), y1(:), h2(:), y2(:)];
0289     ch_act  = 1:Nch;
0290     try_act = 1:Ntry;
0291     try_err = [];
0292     flg = [];
0293     %  Plot bad trial & channel for set of thresholds
0294     vb_plot_max_ratio_image(ratio1,thred_val_lst1)
0295     fig_handler(2)=gcf;
0296     vb_plot_max_ratio_image(ratio2,thred_val_lst2)
0297     fig_handler(3)=gcf;
0298     vb_plot_channel_trial(ratio1,ratio2, ...
0299             thred_val_lst1,thred_val_lst2,prob_val,thred_num)
0300     fig_handler(4)=gcf;
0301 %----- Just save trial statics
0302 case    10
0303     ch_act  = 1:Nch;
0304     try_act = 1:Ntry;
0305     try_err = [];
0306     flg = [];
0307 end
0308 
0309 Ntry = length(try_act);
0310 Nact = length(ch_act);
0311 
0312 if job_mode==1 || job_mode==2 || job_mode==3 || job_mode==4 || job_mode==5 || job_mode==6
0313     fprintf('Active %d channel & %d epoch\n',Nact, Ntry)
0314     fprintf('Err trial= %d\n',length(try_err));
0315 end
0316 
0317 % channel & trial statistics
0318 fileinfo.ymax1 = ymax1;
0319 fileinfo.ymax2 = ymax2;
0320 fileinfo.ystd1 = ystd1;
0321 fileinfo.ystd2 = ystd2;
0322 
0323 fileinfo.channel_id = ch_act;
0324 fileinfo.act_trial  = try_act;
0325 fileinfo = vb_fileinfo_active_field_convert_to(data_format, fileinfo, datafile{1});
0326 
0327 if ~exist('savefile','var')||isempty(savefile), return;end
0328 
0329 vb_save_active_info(fileinfo, [savefile '.info.mat']);
0330 
0331 return
0332 end
0333 %% ---- END ---- %%
0334 
0335 % Inner function for stat mode
0336 function [ylist,h,yy] = inner_plot_stat_hist(y,plist)
0337 %
0338 % Copyright (C) 2011, ATR All Rights Reserved.
0339 % License : New BSD License(see VBMEG_LICENSE.txt)
0340 
0341 [Nch, Ntry] = size(y);
0342 Nbin = 500;
0343 ymax = max(y(:));
0344 
0345 [h, yy] =hist(y(:),Nbin);
0346 
0347 [ylist, p] = vb_find_hist_threshold(h,yy,plist);
0348 
0349 % sum( pd * dy) = 1
0350 dy = yy(2) - yy(1);
0351 pd = p / dy;
0352 
0353 hold on
0354 type  = {'r--', 'b--', 'm--', 'c--'}; % Max num of plist=4
0355 
0356 N = length(ylist);
0357 lgn = cell(N-1,1);
0358 % Plot thresholds
0359 for n=1:N
0360     plot([ylist(n) ylist(n)],[0 max(pd)], type{n})
0361     lgn{n} = ['p=' num2str(plist(n))];
0362 end
0363 % legend(lgn)
0364 xlim([0 ymax])
0365 
0366 % plot histgram
0367 plot(yy,pd);
0368 end