detect_blink_from_eog

 Detect trials in which blinks occurred for visual task

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

0001 function detect_blink_from_eog(p)
0002 % Detect trials in which blinks occurred for visual task
0003 %
0004 % Copyright (C) 2011, ATR All Rights Reserved.
0005 % License : New BSD License(see VBMEG_LICENSE.txt)
0006 
0007 disp(mfilename);
0008 
0009 % Set parameters
0010 from = -0.5;% sec
0011 to = 0.5;% sec
0012 load_spec.ChannelName = {'EOG1', 'EOG2'};
0013 lowpass = 5;% Hz
0014 th_p = 0.05;% Threshold for p
0015 
0016 file_list = {'v1', 'v2'};% File list to detect blink
0017 Nfile = length(file_list);
0018 
0019 % Determine threshold
0020 e_all = [];
0021 for file = 1:Nfile
0022     % Load EOG data
0023     trial_file = fullfile(p.proj_root, p.eeg_dirname, p.trial_dirname, [file_list{file} '.eeg.mat']);
0024     [eog, ~, time_info] = vb_load_meg_data(trial_file, load_spec);
0025     
0026     % Apply low-pass filter to EOG data
0027     for tr = 1:size(eog,3)
0028         eog(:,:,tr) = eegfilt(eog(:,:,tr), time_info.sample_frequency, 0, lowpass);
0029     end
0030     
0031     % Calculate power of EOG data
0032     e = squeeze(sum(abs(eog).^2, 1));
0033     e_all = [e_all; e(:)];
0034 end
0035 
0036 % Determine threshold
0037 ix = round(th_p*length(e_all));
0038 se=sort(e_all, 'descend');
0039 th=se(ix);
0040 
0041 [~,a1] = min(abs(time_info.time-from));
0042 [~,a2] = min(abs(time_info.time-to));
0043 
0044 % Detect blink trials for visual task
0045 for file=1:Nfile
0046     trial_file = fullfile(p.proj_root, p.eeg_dirname, p.trial_dirname, [file_list{file} '.eeg.mat']);
0047     eog = vb_load_meg_data(trial_file, load_spec);
0048     
0049     % Apply low-pass filter to EOG data
0050     Ntr = size(eog, 3);
0051     eog1 = eog;
0052     for tr = 1:Ntr
0053         eog1(:,:,tr) = eegfilt(eog(:,:,tr), time_info.sample_frequency, 0, lowpass);
0054     end
0055     
0056     % Calculate power of EOG data
0057     e = squeeze(sum(abs(eog1).^2,1));
0058     
0059     % Detect blink trials
0060     blink = zeros(1,Ntr);
0061     for tr = 1:Ntr
0062         if max(e(a1:a2, tr)) > th
0063             blink(1,tr) = 1;
0064         else
0065             blink(1,tr) = 0;
0066         end
0067     end
0068     
0069     % Dispaly and save blink information
0070     disp(['Number of blink trials = ' num2str(length(find(blink==1))) '/' num2str(length(blink)) '.'])
0071     load(trial_file, 'EEGinfo')
0072     EEGinfo.blink = blink;
0073     vb_save(trial_file, 'EEGinfo');
0074 end
0075 
0076 % For auditory and somatosensory tasks, there are no blink trials
0077 % because eyes were closed during these tasks.
0078 file_list = {'a1', 's1'};
0079 Nfile = length(file_list);
0080 for file = 1:Nfile
0081     trial_file = fullfile(p.proj_root, p.eeg_dirname, p.trial_dirname, [file_list{file} '.eeg.mat']);
0082     load(trial_file, 'EEGinfo')
0083     EEGinfo.blink = zeros(1, EEGinfo.Nrepeat);
0084     vb_save(trial_file, 'EEGinfo');
0085 end
0086