opticat_script_alternative_steeperfilter.m

% Alternative version of the "OPTICAT" script that uses a "steeper" 
% high-pass filter with a custom, more narrow transition bandwidth
%
% Explanation: Procedures in the Dimigen et al. (2020) paper use 
% EEGLAB's pop_eegfiltnew() function with its default transition bandwidth 
% settings. Per default, this function filters the EEG using a relatively wide
% and frequency-dependent transition bandwidth 
% (see https://github.com/downloads/widmann/firfilt/firfilt.pdf).
%
% Removing low frequencies from the training data likely means that
% low-frequency neural activity (e.g. < 1-2 Hz) will not be properly modelled 
% in the ICA decomposition. In the paper, I did not find that this had a
% problematic effect on the stimulus-ERPs, that is, I did not see any 
% increase in overcorrection of the stimulus-ERP if the *training
% data* was high pass-filtered aggressively (e.g. at 2-2.5 Hz) with the 
% EEGLAB defaults.

% However, Supplementary Figure S7 suggests that the problematic frequency
% content in the training data is below 1 Hz, at least for Scene viewing.
% Thus, a potentially safer alternative is to filter the EEG at a lower cutoff
% (e.g. at 0.75 Hz, see Supplementary Figure S7) but with a "steeper" filter.
% 
% This should remove less of the low-frequency neural activity
% from the training data while producing comparable correction results 
% (cf. Supplementary Figure S7) as the optimal filter settings described with 
% the EEGLAB filter defaults detailed in the paper. Preserving low-frequency 
% activity in the training data would be beneficial if the user wants to work 
% with the decomposed source waveforms in source space (rather than backproject 
% the corrected data).
%
% The script below implements the OPTICAT procedures with a custom, 
% "steeper" filter setting. For more documentation on the OPTICAT procedure
% itself, see the "opticat_script.m" in this folder.
%
% I would like to thank an anonymous contributor for helpful input
% on filter design.

clear

%% add the firfilt toolbox by A. Widmann (https://github.com/widmann/firfilt)
addpath M:\Dropbox\tools\firfilt-master

%% Load your EEG dataset (can be continuous or epoched)
% Note: This dataset in EEGLAB format needs to already include 
% 'saccade' and 'fixation' events in the EEG.event structure.
% These can be added, for example, with the EYE-EEG toolbox. 
EEG = pop_loadset('filename','C:/myEEGdata.set');


OW_FACTOR        = 1    % value for overweighting of SPs (1 = add data corresponding to 100% of original data length)
REMOVE_EPOCHMEAN = true % mean-center the appended peri-saccadic epochs? (strongly recommended)
EEG_CHANNELS     = 1:45 % indices of all EEG channels (exclude any eye-tracking channels here)
                        % I recommend to also include EOG channels if they
                        % were also recorded against the common reference % (not: bipolar)

% ------------------START: STEEPER HIGH-PASS FILTER -----------------------
CUTOFF           = 0.75 % -6 dB cutoff (in Hz)
                        % Note: possibly go even lower, e.g. to 0.5 Hz if your priority is
                        % to preserve low-frequency neural activity in the decomposition. 
                        % However, compare Suppl. Figure S7 for effects on correction quality
TBW              = 0.5  % transition bandwidth (in Hz)

% compute filter order
filterorder = firwsord('hamming', EEG.srate, TBW);

% design windowed sinc filter
b = firws(filterorder, CUTOFF / (EEG.srate/2), 'high', windows('hamming', filterorder+1));

% get frequency response
[h,f] = freqz(b,1,16384,EEG.srate);

% visualize filter characteristics
figure
subplot(2,1,1); hold on; title('linear')
plot(f,abs(h))
xlim([0 2]), ylim([0 1.05])
subplot(2,1,2); hold on; title('dB')
plot(f,20*log10(abs(h)))
xlim([0 2]), ylim([-80 0])
xlabel('Hz')
% -------------------END: STEEPER HIGH-PASS FILTER ------------------------


%% Create copy of data used as training data & high pass-filter it
EEG_training = EEG;
EEG_training = firfilt(EEG_training,b,[]); % filter the training data using firfilt()

%% Cut training data into epochs, e.g. around stimulus onsets (as in Dimigen, 2020)
EEG_training = pop_epoch(EEG_training,{'S 11','S 12','S 13'},[-0.2 1.8]);

%% Overweight spike potentials
% Repeatedly append intervals around saccade onsets (-20 to +10 ms) to training data
EEG_training = pop_overweightevents(EEG_training,'saccade',[-0.02 0.01],OW_FACTOR,REMOVE_EPOCHMEAN);

%% Run ICA
fprintf('\nRunning ICA on optimized training data...')
EEG_training = pop_runica(EEG_training,'extended',1,'interupt','on','chanind',EEG_CHANNELS); % or use binary ICA for more speed

%% Remember ICA weights & sphering matrix
wts = EEG_training.icaweights;
sph = EEG_training.icasphere;

%% Remove any existing ICA solutions from your original dataset
EEG.icaact      = [];
EEG.icasphere   = [];
EEG.icaweights  = [];
EEG.icachansind = [];
EEG.icawinv     = [];

%% Transfer unmixing weights
fprintf('\nTransfering ICA weights from training data to original data...')
EEG.icasphere   = sph;
EEG.icaweights  = wts;
EEG.icachansind = EEG_CHANNELS;
EEG = eeg_checkset(EEG); % let EEGLAB re-compute EEG.icaact & EEG.icawinv


fprintf('\nIdentifying ocular ICs via saccade/fixation variance-ratio threshold...')

%% Eye-tracker-guided selection of ICs
IC_THRESHOLD     = 1.1;   % variance ratio threshold (determined as suitable in Dimigen, 2020)
SACC_WINDOW      = [5 0]; % saccade window (in samples!) to compute variance ratios (see Dimigen, 2020)
PLOTFIG          = true;  % plot a figure visualizing influence of threshold setting?
ICPLOTMODE       = 2;     % plot component topographies (inverse weights)? (2 = only plot "bad" ocular ICs)
FLAGMODE         = 3;     % overwrite existing rejection flags? (3 = yes)

%% Automatically flag ocular ICs (Plöchl et al., 2012)
[EEG, varratiotable] = pop_eyetrackerica(EEG,'saccade','fixation',SACC_WINDOW,IC_THRESHOLD,FLAGMODE,PLOTFIG,ICPLOTMODE);

%% Remove flagged ocular ICs
badcomps = EEG.reject.gcompreject;
EEG      = pop_subcomp(EEG,find(badcomps)); % remove them



%%  ----- SIMPLE CORRECTION QUALITY CHECK (TEST FOR UNDERCORRECTION) -----

%% Evaluate results: Compute saccade-locked ERPs for saccades of a single orientation
EEG_sac = pop_epoch(EEG,{'saccade'},[-0.2 0.6]); % cut epochs around saccades

%% Get saccades of certain orientations (here: plusminus 30 deg)
% (note that the origin of most eye-tracker systems is *upper* left screen corner)
% Right-going:      0 degree
% Left-going:     180 degree
% Up-going:       -90 degree (for most ET systems)
% Down-going:      90 degree (for most ET systems)

% Extract saccade angles
for e = 1:length(EEG_sac.epoch)
    ix = find([EEG_sac.epoch(e).eventlatency{:}] == 0);
    sac_angles(e) = cell2mat(EEG_sac.epoch(e).eventsac_angle(ix(1)));
end

% Plot saccade angles
[t,r] = rose(sac_angles*pi/180,36); % angle in radians, plot 10° bins
figure('name','saccade angles')
h2 = polar(t,r,'r-');
set(gca,'ydir','reverse') % origin of coordinate system = upper left corner

%% Extract left/right/up/down-going saccades (plusminus 30°)
ix_R = find( sac_angles >  -30 & sac_angles <   30);
ix_L = find( sac_angles >  150 | sac_angles < -150);
ix_U = find( sac_angles > -120 & sac_angles <  -60);
ix_D = find( sac_angles >   60 & sac_angles <  120);

%% Cut epochs around saccades (here: rightward sacc. only, adapt as you wish)
try
    fprintf('\nNumber of \"rightwards\" saccades: %i',length(ix_R))
    EEG_sac_R = pop_select(EEG_sac,'trial',ix_R); % select rightwards sacc.
    EEG_sac_R = pop_rmbase(EEG_sac_R,[-100 0],[]);   % subtract pre-sacc. baseline
catch err
    EEG_sac_R = [];
end

%% Plot saccade-ERP to check for residual artifacts
% (here: for rightward saccades)
figure('name','Control plot: Undercorrection')
plot(EEG_sac_R.times, mean(EEG_sac_R.data(EEG_CHANNELS,:,:),3))
ylabel('Saccade-related ERP')
xlabel('Time after saccade [ms]')