vb_megfile_remove_trial_ext

 Remove trials based on external signals

 --- Syntax
 vb_megfile_remove_trial_ext(megfile, megfile_new, parm)

 --- History
 2008-05-09 Taku Yoshioka
 2010-01-22 (Matt) made more efficient
 2010-01-26 (Sako) renamed

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

0001 function vb_megfile_remove_trial_ext(megfile,megfile_new)
0002 % Remove trials based on external signals
0003 %
0004 % --- Syntax
0005 % vb_megfile_remove_trial_ext(megfile, megfile_new, parm)
0006 %
0007 % --- History
0008 % 2008-05-09 Taku Yoshioka
0009 % 2010-01-22 (Matt) made more efficient
0010 % 2010-01-26 (Sako) renamed
0011 %
0012 % Copyright (C) 2011, ATR All Rights Reserved.
0013 % License : New BSD License(see VBMEG_LICENSE.txt)
0014 
0015 % Default parameters
0016 threshold = 8;
0017 verbose = 2;   % verbose level
0018 if nargin<1, help vb_megfile_remove_trial_ext; end
0019 
0020 % Get information of MEG data
0021 MEGinfo = vb_load_measurement_info(megfile);
0022 ActTrial = vb_info_get_active_trial(MEGinfo);
0023 ExtChannelName = vb_megfile_get_channel_label_extra(megfile);
0024 
0025 clear load_spec;
0026 load_spec.ActiveChannel = false;
0027 load_spec.ChannelType = 'EXTRA';
0028 [full_ext_data] = vb_load_meg_data(megfile,load_spec);
0029 
0030 ext_data = full_ext_data(:,:,ActTrial==1);
0031 
0032 clear load_spec;
0033 
0034 % Find external channels having no information
0035 Nc = size(ext_data,1);
0036 ix_channel = [];
0037 for i=1:Nc
0038   tmp = sum(abs(ext_data(i,:,:)),3);
0039   if std(tmp)~=0, ix_channel = [ix_channel; i]; end
0040 end
0041 
0042 % Plot signal time courses of all external channels
0043 Nxmax = 3;
0044 Nymax = 4;
0045 j = 1;
0046 while j>0
0047   figure;
0048   ixx_channel ...
0049       = (1+(j-1)*Nxmax*Nymax):min([length(ix_channel) j*Nxmax*Nymax]);
0050   Nx = min([length(ixx_channel) Nxmax]);
0051   Ny = ceil(length(ixx_channel)/Nx);
0052   for i=1:length(ixx_channel)
0053     subplot(Ny,Nx,i);
0054     ix_channel_org = ix_channel(ixx_channel(i));
0055     plot(ext_data(ix_channel_org,:,1));
0056     title(['Channel ' num2str(ix_channel_org) ': ' ...
0057            ExtChannelName{ix_channel_org}]);
0058   end
0059   if length(ix_channel)<=j*(Nxmax*Nymax), j=0;
0060   else j = j+1; end
0061 end
0062 
0063 % Select external channels
0064 answer = inputdlg(['Input channel number(s) for choosing trials' ...
0065                    ' removed\n (e.g., [1 4 5])'], ...
0066                   'Input external channels');
0067 close all;
0068 if isempty(answer), return; 
0069 elseif ~isempty(answer{1}), ix_channel = str2num(answer{1}); end
0070 
0071 % Input threshold
0072 for i=1:length(ix_channel)
0073   tmp = abs(ext_data(ix_channel(i),:,:));
0074   hist(tmp,100);
0075   th(i) = min([threshold*median(tmp(:)) max(tmp(:))]);
0076   yrange = ylim;
0077   line([th(i) th(i)],[yrange(1) yrange(2)],'LineStyle','--','Color','k');
0078   answer = inputdlg(['Input threshold value for channel' ...
0079                      num2str(ix_channel(i))],'Threshold for noisy signal', ...
0080                     1,{num2str(th(i),'%.1f')});
0081   close;
0082   if ~isempty(answer) && ~isempty(answer{1}), 
0083     th(i) = str2num(answer{1}); 
0084   end
0085 end
0086 
0087 
0088 % Find noisy trials
0089 ActTrial = vb_info_get_active_trial(MEGinfo);
0090 for i=1:MEGinfo.Nrepeat
0091   if ~ActTrial(i), continue; end 
0092   
0093   ext_data = full_ext_data(:,:,i);
0094   
0095   remove_trial = false;
0096   for j=1:length(ix_channel)
0097     if max(abs(ext_data(ix_channel(j),:)))>th(j), remove_trial=true; end
0098   end
0099   
0100   if remove_trial, 
0101     if verbose==2, 
0102       
0103       figure;
0104       for j=1:length(ix_channel)
0105         title(['Channel ' num2str(ix_channel(j))]);
0106         subplot(length(ix_channel),1,j);
0107         plot(ext_data(ix_channel(j),:));
0108       end
0109       remove_trial = questdlg(['Remove trial ' num2str(i) '?'], ...
0110                               'Remove trial','Yes','No','Yes');
0111       close; close;
0112     else
0113       remove_trial = 'Yes';
0114     end
0115   end
0116   if strcmp(remove_trial,'Yes'), ActTrial(i) = false; end
0117 end
0118 
0119 % Update MEGinfo struct
0120 update_spec.preprocess_parm.parm = [];
0121 update_spec.active_trial = ActTrial;
0122 vb_msrmnt_make_preprocessed_file(megfile,update_spec,megfile_new);
0123