vb_megfile_remove_noisysensor

 Remove sensors that is possibly broken. 

 --- Syntax
 vb_megfile_remove_noisysensor(megfile,megfile_new)

 --- History
 2008-05-09 Taku Yoshioka
 2008-06-18 Taku Yoshioka
 2010-01-19 (Sako) replaced log_meg_info with log_measurement_info
 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_noisysensor(megfile,megfile_new)
0002 % Remove sensors that is possibly broken.
0003 %
0004 % --- Syntax
0005 % vb_megfile_remove_noisysensor(megfile,megfile_new)
0006 %
0007 % --- History
0008 % 2008-05-09 Taku Yoshioka
0009 % 2008-06-18 Taku Yoshioka
0010 % 2010-01-19 (Sako) replaced log_meg_info with log_measurement_info
0011 % 2010-01-22 (Matt) made more efficient
0012 % 2010-01-26 (Sako) renamed
0013 %
0014 % Copyright (C) 2011, ATR All Rights Reserved.
0015 % License : New BSD License(see VBMEG_LICENSE.txt)
0016 
0017 const = vb_define_const;
0018 
0019 % Default parameters
0020 threshold = 9;   % ratio to the median of signals
0021 NnoisyTrial = 5; % If the number of noisy trial exceeds this value, the
0022                  % sensor is removed.
0023 verbose = 2;     % verbose level
0024 if nargin<1, help vb_megfile_remove_noisysensor; end
0025 
0026 % Get information of MEG data
0027 MEGinfo = vb_load_measurement_info(megfile);
0028 ActTrial = vb_info_get_active_trial(MEGinfo);
0029 ActChannel = vb_info_get_active_channel(MEGinfo);
0030 ChannelName = vb_info_get_channel_label(MEGinfo);
0031 
0032 ch_info = vb_load_channel_info(megfile,'PLANAR');
0033 ActChannel(ch_info.ID)=0;
0034 
0035 load_spec.ChannelType = const.DATA_TYPE_MAIN;
0036 full_meg_data = vb_load_meg_data(megfile,load_spec);  
0037 
0038 % Create histogram
0039 for i=1:MEGinfo.Nrepeat
0040   if ~ActTrial(i), continue; end
0041   load_spec.TrialNumber = i;
0042   load_spec.ActiveChannel = true;
0043   
0044   meg_data = full_meg_data(ActChannel==1,:,i);
0045   
0046   if i==1, 
0047     m_max = max(abs(meg_data(:)));
0048     dx = m_max/100;
0049     x = dx:dx:m_max;
0050     y = hist(abs(meg_data(:)),x);
0051   else
0052     y = y+hist(abs(meg_data(:)),x);
0053   end
0054 end
0055 
0056 % Calculate threshold value
0057 Ndata = sum(y(:));
0058 Nmed = 0;
0059 x = x*1e15;
0060 for i=1:length(y)
0061   Nmed = Nmed+y(i);
0062   if Nmed>=Ndata/2, 
0063     th = threshold*x(i);
0064     break;
0065   end
0066 end
0067 
0068 % Plot histogram
0069 bar(x,y)
0070 yrange = ylim;
0071 line([th th],[yrange(1) yrange(2)],'LineStyle','--','Color','k');
0072 title('Histogram of signal amplitude');
0073 xlabel('MEG signal [fT]');
0074 
0075 % Show dialog to set threshold value
0076 answer = inputdlg('Input threshold value','Threshold for noisy signal', ...
0077                   1,{num2str(th,'%.1f')});
0078 close;
0079 if isempty(answer), return; end
0080 if ~isempty(answer{1}), th = str2num(answer{1}); end
0081 
0082 % Find noisy sensors
0083 ix_channel = [];
0084 clear load_spec;
0085 for i=1:MEGinfo.Nchannel
0086   if ~ActChannel(i), continue; end
0087   load_spec.ChannelName{1} = ChannelName{i};
0088   load_spec.ChannelSwitch = true;
0089   load_spec.ActiveTrial = true;  
0090   
0091   meg_data = full_meg_data(i,:,:);
0092   meg_data = meg_data*1e15;
0093   
0094   ix_trial = [];
0095   max_signal = [];
0096   for j=1:size(meg_data,3)
0097     tmp = abs(meg_data(:,:,j));
0098     if max(tmp)>th, 
0099       ix_trial = [ix_trial; j]; 
0100       max_signal = [max_signal; max(tmp)];
0101     end
0102   end
0103   
0104   [tmp,ix] = sort(-1*max_signal);
0105   ix_trial = ix_trial(ix);
0106   
0107   if length(ix_trial)>NnoisyTrial, 
0108     if verbose==2,
0109       N = min([5 length(ix_trial)]);
0110       for j=1:N
0111         subplot(N,1,j);
0112         plot(meg_data(:,:,ix_trial(j)));
0113       end
0114       title('Signal time courses');
0115       ylabel('[fT]');
0116       remove_sensor = questdlg(['Remove sensor ' ChannelName{i} '?'], ...
0117                                'Remove sensor','Yes','No','Yes');
0118       close;
0119     else
0120       remove_sensor = 'Yes';
0121     end
0122     
0123     if strcmp(remove_sensor,'Yes'), ActChannel(i) = false; end
0124   end
0125 end
0126 
0127 % Update MEGinfo struct
0128 update_spec.preprocess_parm.parm = [];
0129 update_spec.active_channel = ActChannel;
0130 vb_msrmnt_make_preprocessed_file(megfile,update_spec,megfile_new);