vb_job_bad_trial_stat

 Plot bad channel/trials for set of threshold
 --- Usage
  vb_job_bad_trial_stat(datafile,job_mode,actfile)
  vb_job_bad_trial_stat(datafile,job_mode,actfile,stat_file)
  vb_job_bad_trial_stat(datafile,job_mode,actfile,stat_file,  ...
    thred_val_lst,prob_val_lst,plist,thred_ch,std_mode)
 --- Input
 datafile: cell string for multiple MEG/EEG epoch data files
 job_mode: 
  = 0 : Plot ymax histgram
  = 1 : Plot bad trial & channel for set of thresholds
 [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 list for bad channel selection
 - thred_val_lst = [8:11 ; 7:10]; 
   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_lst  = [0.05   0.01   0.001];
   threshold probability list 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-8-10 Masa-aki Sato

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

0001 function  [fileinfo,flg,data] = ...
0002   vb_job_bad_trial_stat(datafile,job_mode,actfile,stat_file,  ...
0003     thred_val_lst,prob_val_lst,plist,thred_ch,std_mode)
0004 % Plot bad channel/trials for set of threshold
0005 % --- Usage
0006 %  vb_job_bad_trial_stat(datafile,job_mode,actfile)
0007 %  vb_job_bad_trial_stat(datafile,job_mode,actfile,stat_file)
0008 %  vb_job_bad_trial_stat(datafile,job_mode,actfile,stat_file,  ...
0009 %    thred_val_lst,prob_val_lst,plist,thred_ch,std_mode)
0010 % --- Input
0011 % datafile: cell string for multiple MEG/EEG epoch data files
0012 % job_mode:
0013 %  = 0 : Plot ymax histgram
0014 %  = 1 : Plot bad trial & channel for set of thresholds
0015 % [actfile   '.info.mat']
0016 %  : file name to save bad channel/trial info for batch mode
0017 %    or save trial statics without channel/trial rejection
0018 % [stat_file '.info.mat']
0019 %  : file name having trial statics made at the previous job
0020 %
0021 % ---- Threshold value list for bad channel selection
0022 % - thred_val_lst = [8:11 ; 7:10];
0023 %   threshold value for max(y) & diff(y) in one trial
0024 %   if (ymax/ystd) > thred_val(1) | difmax/difstd > thred_val(2),
0025 %      it is marked as bad
0026 % - prob_val_lst  = [0.05   0.01   0.001];
0027 %   threshold probability list to select bad channel
0028 %   if more than (prob_val*Ntrial) trials are bad, the channel is rejected
0029 % - thred_ch  = 1/4;
0030 %   if more than (Nch*thred_ch) channels are bad, the trial is rejected,
0031 %   before checking bad channel.
0032 % - std_mode  = 1;
0033 %   = 0: std is averaged over all trials
0034 %   = 1: std is averaged within one session (one BDF file)
0035 %
0036 %
0037 % 2009-8-10 Masa-aki Sato
0038 %
0039 % Copyright (C) 2011, ATR All Rights Reserved.
0040 % License : New BSD License(see VBMEG_LICENSE.txt)
0041 
0042 if ~exist('std_mode','var') , std_mode = 1; end;
0043 
0044 if ~isempty(datafile) && ~iscell(datafile),
0045     tmp = datafile;
0046     datafile = {tmp};
0047 end;
0048 
0049 if exist('stat_file','var') && exist([stat_file '.info.mat'],'file')
0050     % Load statics info (ratio. ymax)
0051     load([stat_file '.info.mat'])
0052 else
0053     % Get file info for multiple session MEG/EEG files
0054     % All trials are loaded
0055     fileinfo = vb_get_multi_fileinfo(datafile);
0056 end
0057 
0058 
0059 % ystd: amplitude mean calculated by histrgam for all trials in one channel
0060 % ymax : Nch x Ntrial
0061 % ystd : Nch x 1
0062 
0063 if isfield(fileinfo,'ymax1')
0064     ymax1 = fileinfo.ymax1;
0065     ymax2 = fileinfo.ymax2;
0066     ystd1 = fileinfo.ystd1;
0067     ystd2 = fileinfo.ystd2;
0068     if ~isempty(datafile) 
0069         fileinfo = vb_get_multi_fileinfo(datafile);
0070     end
0071 else
0072     ymax1 = [];
0073     ymax2 = [];
0074     ystd1 = [];
0075     ystd2 = [];
0076     fprintf('Load data ')
0077     for n=1:length(datafile)
0078         % Load MEG/EEG data
0079         fprintf('-')
0080         temp = vb_load_meg_data(datafile{n});
0081         [max1,max2,std1,std2] = vb_channel_statics(temp);
0082         ymax1 = [ymax1 , max1];
0083         ymax2 = [ymax2 , max2];
0084         ystd1 = [ystd1 , std1];
0085         ystd2 = [ystd2 , std2];
0086     end
0087     fprintf('\n')
0088 end
0089 
0090 
0091 if std_mode == 0 
0092     ystd1 = mean(ystd1,2);
0093     ystd2 = mean(ystd2,2);
0094 end
0095 
0096 ratio1 = ymax1;
0097 ratio2 = ymax2;
0098 
0099 ix1  = find(ystd1 > 0);
0100 ix2  = find(ystd2 > 0);
0101 
0102 if size(ystd1,2) == 1
0103     ratio1(ix1,:) = vb_repmultiply(ymax1(ix1,:), 1./ystd1(ix1));
0104     ratio2(ix2,:) = vb_repmultiply(ymax2(ix2,:), 1./ystd2(ix2));
0105 else
0106     ratio1(ix1) = ymax1(ix1)./ystd1(ix1);
0107     ratio2(ix2) = ymax2(ix2)./ystd2(ix2);
0108 end
0109 %if std_mode == 0
0110 %    ratio1 = vb_repmultiply(ymax1, 1./max(mean(ystd1,2), eps));
0111 %    ratio2 = vb_repmultiply(ymax2, 1./max(mean(ystd2,2), eps));
0112 %elseif any(size(ymax1) - size(ystd1))
0113 %    ratio1 = vb_repmultiply(ymax1, 1./max(ystd1, eps));
0114 %    ratio2 = vb_repmultiply(ymax2, 1./max(ystd2, eps));
0115 %else
0116 %    ratio1 = ymax1./max(ystd1, eps);
0117 %    ratio2 = ymax2./max(ystd2, eps);
0118 %end
0119 
0120 flg = [];
0121 data = [];
0122 
0123 [Nch,Ntry] = size(ratio1);
0124 thred_num  = Nch * thred_ch;
0125 
0126 if ~exist('thred_val_lst','var')||isempty(thred_val_lst)
0127     thred_val_lst1 = [8:11]; % [5:10] for ratio
0128     thred_val_lst2 = [7:10]; % for ymax1
0129 else
0130     thred_val_lst1 = thred_val_lst(1,:);
0131     thred_val_lst2 = thred_val_lst(2,:);
0132 end
0133 if ~exist('prob_val_lst','var')||isempty(prob_val_lst)
0134     prob_val_lst   = [0.02:0.01:0.05]; % [2:8]
0135 end
0136 if ~exist('plist','var')||isempty(plist)
0137     plist = [0.05   0.01   0.001];
0138 end
0139 
0140 switch job_mode
0141 case    0
0142     % Plot ymax histgram
0143     subplot 121
0144     [ylist1,h1,y1] = vb_plot_stat_hist(ratio1 ,plist);
0145     hold on
0146     title('Max/std of amplitude')
0147     subplot 122
0148     [ylist2,h2,y2] = vb_plot_stat_hist(ratio2 ,plist);
0149     hold on
0150     title('Max/std of derivative')
0151     
0152     % set output variable
0153     flg  = [ylist1; ylist2];
0154     data = [h1(:), y1(:), h2(:), y2(:)];
0155     
0156 case    1
0157     % Plot ymax histgram
0158     subplot 121
0159     [ylist1,h1,y1] = vb_plot_stat_hist(ratio1 ,plist);
0160     hold on
0161     title('Max/std of amplitude')
0162     subplot 122
0163     [ylist2,h2,y2] = vb_plot_stat_hist(ratio2 ,plist);
0164     hold on
0165     title('Max/std of derivative')
0166     
0167     % set output variable
0168     flg  = [ylist1; ylist2];
0169     data = [h1(:), y1(:), h2(:), y2(:)];
0170     
0171     %  Plot bad trial & channel for set of thresholds
0172     vb_plot_max_ratio_image(ratio1,thred_val_lst1)
0173     vb_plot_max_ratio_image(ratio2,thred_val_lst2)
0174     vb_plot_channel_trial(ratio1,ratio2, ...
0175             thred_val_lst1,thred_val_lst2,prob_val_lst,thred_num)
0176 end
0177 
0178 ch_act  = 1:Nch;
0179 try_act = 1:Ntry;
0180 try_err = [];
0181 
0182 % channel & trial statistics
0183 fileinfo.ymax1 = ymax1;
0184 fileinfo.ymax2 = ymax2;
0185 fileinfo.ystd1 = ystd1;
0186 fileinfo.ystd2 = ystd2;
0187 
0188 fileinfo.channel_id = ch_act;
0189 fileinfo.act_trial = try_act;
0190 fileinfo.err_trial = try_err;
0191 
0192 if ~exist('actfile','var')||isempty(actfile), return;end
0193 
0194 
0195 fprintf('Save trial statistics info to: %s\n',[actfile '.info.mat'])
0196 
0197 vb_fsave([actfile '.info.mat'],'fileinfo');