vb_test_bad_trial_stat

 vb_test_bad_trial
 check bad trials

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

0001 % vb_test_bad_trial
0002 % check bad trials
0003 %
0004 % Copyright (C) 2011, ATR All Rights Reserved.
0005 % License : New BSD License(see VBMEG_LICENSE.txt)
0006 clear all
0007 
0008 % [actfile   '.info.mat']
0009 %  : file name to save bad channel/trial info for batch mode
0010 %    or save trial statics without channel/trial rejection
0011 % [stat_file '.info.mat']
0012 %  : file name having trial statics made at the previous job
0013 %
0014 % job_mode for 'vb_job_bad_trial_stat'
0015 %  = 0 : Plot ymax histgram
0016 %  = 1 : Plot bad trial & channel for set of thresholds
0017 %  trial statics are saved to actfile
0018 
0019 job_mode  = 1;
0020 data_mode = 2;
0021 
0022 switch    data_mode
0023 case    1
0024     % MEG epoch data
0025     root_dir  = [getenv('MATHOME') '/Sekiyama'];
0026     data_file = {'AM090311a_voice_mik_B'};
0027     ext = '_epoch.meg.mat';
0028     stat_file = [];
0029     actfile   = [root_dir '/' 'test' ];
0030 case    2
0031     % EEG epoch data
0032     root_dir  = [getenv('MATHOME') '/BCI_SK/epoch_data01/080903'];
0033     data_file = {'Sako_20080903_01'; 'Sako_20080903_02'};
0034     ext = '_epoch.eeg.mat';
0035     stat_file = [];
0036     actfile   = [root_dir '/' 'test' ];
0037 case    3
0038     % Statics data
0039     root_dir  = [getenv('MATHOME') '/BCI_SK/epoch_data01'];
0040     data_file = [];
0041     ext = '_epoch.eeg.mat';
0042     stat_file = {...
0043                  [root_dir '/080903/SK_all'];
0044                  [root_dir '/080904/SK_all'];
0045                  [root_dir '/080909/SK_all'];
0046                  [root_dir '/080910/SK_all'];
0047                  [root_dir '/080917/SK_all'];
0048                  [root_dir '/080918/SK_all'];
0049                  [root_dir '/080924/SK_all']; };
0050     actfile   = [];
0051 end
0052 
0053 Nfile = length(data_file);
0054 
0055 for n=1:Nfile
0056     data_file{n} = [root_dir '/' data_file{n}  ext];
0057 end
0058 
0059 % ---- Threshold value list for bad channel selection
0060 % - thred_val_lst = [8:11 ; 7:10];
0061 %   threshold value for max(y) & diff(y) in one trial
0062 %   if (ymax/ystd) > thred_val(1) | difmax/difstd > thred_val(2),
0063 %      it is marked as bad
0064 % - prob_val_lst  = [0.05   0.01   0.001];
0065 %   threshold probability list to select bad channel
0066 %   if more than (prob_val*Ntrial) trials are bad, the channel is rejected
0067 % - thred_ch  = 1/4;
0068 %   if more than (Nch*thred_ch) channels are bad, the trial is rejected,
0069 %   before checking bad channel.
0070 
0071 thred_val_lst = [8:11 ; 7:10]; 
0072 prob_val_lst  = [0.05   0.01   0.001];
0073 thred_ch = 1/4;
0074 
0075 % P-value to determine threshold value from Max-value histgram
0076 plist = [0.05   0.01   0.005];
0077 
0078 xmax = 20;
0079 xd = 0.25;
0080 
0081 if isempty(stat_file)
0082     % check trial statics for data_file
0083     [fileinfo,y,p] = ...
0084         vb_job_bad_trial_stat(data_file,job_mode,actfile,stat_file,  ...
0085         thred_val_lst,prob_val_lst,plist,thred_ch);
0086     return
0087 else
0088     % check trial statics for multiple stat_file
0089     Nfile = length(stat_file);
0090     Nlist = 1:Nfile;
0091     %Nlist = 1;
0092     x1 = (0:xd:xmax)';
0093     x2 = x1;
0094     h1 = zeros(length(x1),1);
0095     h2 = zeros(length(x2),1);
0096     ysum = zeros(2,length(plist)+1);
0097     for n= Nlist
0098         [fileinfo,y,hn] = ...
0099             vb_job_bad_trial_stat(data_file,job_mode,actfile,stat_file{n},  ...
0100             thred_val_lst,prob_val_lst,plist,thred_ch);
0101         if size(y)==size(ysum)
0102             ysum = ysum + y;
0103         end
0104         if ~isempty(hn)
0105             [h1 ,x1]= vb_merge_histgram(h1,x1,hn(:,1),hn(:,2));
0106             [h2 ,x2]= vb_merge_histgram(h2,x2,hn(:,3),hn(:,4));
0107         end
0108     end
0109     ysum = ysum/Nfile;
0110 end
0111 
0112 if isempty(hn), return; end;
0113 
0114 % summary of multiple stat_file info
0115 y1 = vb_find_hist_threshold(h1,x1,plist);
0116 y2 = vb_find_hist_threshold(h2,x2,plist);
0117 
0118 figure
0119 
0120 subplot 121
0121 plot(x1,h1)
0122 hold on
0123 for n=1:length(y1)
0124     plot([y1(n) y1(n)], [0 max(h1)],'r--')
0125 end
0126 
0127 subplot 122
0128 plot(x2,h2)
0129 hold on
0130 for n=1:length(y2)
0131     plot([y2(n) y2(n)], [0 max(h2)],'r--')
0132 end
0133 
0134 return
0135 
0136 fprintf('P=[')
0137 fprintf('%g ,' ,p)
0138 fprintf(']\n')
0139 fprintf('Y =[')
0140 fprintf('%g ,' ,ysum(1,:))
0141 fprintf(']\n')
0142 fprintf('Yd=[')
0143 fprintf('%g ,' ,ysum(2,:))
0144 fprintf(']\n')
0145