vb_megfile_remove_noisytrial

 Remove trials that is contaminated by noise. 

 --- Syntax
 vb_megfile_remove_noisytrial(megfile,megfile_new)

 --- History
 2008-05-09 Taku Yoshioka
 2010-01-19 (Sako) modified old fashioned code (vb_load_meg_info, vb_load_meg_data)
 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_noisytrial(megfile,megfile_new)
0002 % Remove trials that is contaminated by noise.
0003 %
0004 % --- Syntax
0005 % vb_megfile_remove_noisytrial(megfile,megfile_new)
0006 %
0007 % --- History
0008 % 2008-05-09 Taku Yoshioka
0009 % 2010-01-19 (Sako) modified old fashioned code (vb_load_meg_info, vb_load_meg_data)
0010 % 2010-01-22 (Matt) made more efficient
0011 % 2010-01-26 (Sako) renamed
0012 %
0013 % Copyright (C) 2011, ATR All Rights Reserved.
0014 % License : New BSD License(see VBMEG_LICENSE.txt)
0015 
0016 % Default parameters
0017 threshold = 9; % ratio to the median of signals
0018 verbose = 2;   % verbose level
0019 if nargin<1, help vb_megfile_remove_noisytrial; end
0020 
0021 % Get information of MEG data
0022 MEGinfo = vb_load_measurement_info(megfile);
0023 ActTrial = vb_info_get_active_trial(MEGinfo);
0024 
0025 ch_info = vb_load_channel_info(megfile,'AXIAL');
0026 loadspec.ChannelName = ch_info.Name;
0027 full_meg_data = vb_load_meg_data(megfile,loadspec);
0028 
0029 % Create histogram
0030 for i=1:MEGinfo.Nrepeat
0031   if ~ActTrial(i), continue; end
0032   
0033   meg_data = full_meg_data(:,:,i);
0034   
0035   if i==1, 
0036     m_max = max(abs(meg_data(:)));
0037     dx = m_max/100;
0038     x = dx:dx:m_max;
0039     y = hist(abs(meg_data(:)),x);
0040   else
0041     y = y+hist(abs(meg_data(:)),x);
0042   end
0043 end
0044 
0045 % Calculate threshold value
0046 Ndata = sum(y(:));
0047 Nmed = 0;
0048 x = x*1e15;
0049 for i=1:length(y)
0050   Nmed = Nmed+y(i);
0051   if Nmed>=Ndata/2, 
0052     th = threshold*x(i);
0053     break;
0054   end
0055 end
0056 
0057 % Plot histogram
0058 bar(x,y)
0059 yrange = ylim;
0060 line([th th],[yrange(1) yrange(2)],'LineStyle','--','Color','k');
0061 title('Histogram of signal amplitude');
0062 xlabel('MEG signal [fT]');
0063 
0064 % Show dialog to set threshold value
0065 answer = inputdlg('Input threshold value','Threshold for noisy signal', ...
0066                   1,{num2str(th,'%.1f')});
0067 close;
0068 if isempty(answer), return; end
0069 if ~isempty(answer{1}), th = str2num(answer{1}); end
0070 close;
0071 
0072 % Find noisy trials
0073 ix_trial = [];
0074 max_signal = [];
0075 for i=1:MEGinfo.Nrepeat
0076   if ~ActTrial(i), continue; end
0077   
0078   meg_data = full_meg_data(:,:,i);
0079   meg_data = meg_data*1e15;
0080 
0081   if max(abs(meg_data(:)))>th, 
0082     max_signal = [max_signal; max(abs(meg_data(:)))];
0083     ix_trial = [ix_trial; i];
0084   end
0085 end
0086 
0087 % Remove noisy trials
0088 for i=length(ix_trial):-1:1
0089   if max_signal(i)>th, 
0090     if verbose==2, 
0091       meg_data = full_meg_data(:,:,ix_trial(i));
0092       
0093       meg_data = meg_data*1e15;
0094       plot(meg_data(:,:)');
0095       title('Signal time course');
0096       ylabel('MEG signal [fT]');
0097       remove_trial = questdlg(['Remove trial ' num2str(ix_trial(i)) '?'], ...
0098                               'Remove trial','Yes','No','Yes');
0099       close;
0100     else
0101       remove_trial = 'Yes';
0102     end
0103     if strcmp(remove_trial,'Yes'), ActTrial(ix_trial(i)) = false; end
0104   end
0105 end
0106 
0107 % Update MEGinfo struct
0108 update_spec.preprocess_parm.parm = [];
0109 update_spec.active_trial = ActTrial;
0110 vb_msrmnt_make_preprocessed_file(megfile,update_spec,megfile_new);
0111