%% Setup

% paths to data
eeg_path = fullfile('raw_data', 'eeg-pc', 'localiser');
beh_path = fullfile('raw_data', 'stim-pc', 'data', 'localiser');

% import eeglab (assumes eeglab has been added to path), e.g.
addpath('C:/EEGLAB/eeglab2020_0')
[ALLEEG, EEG, CURRENTSET, ALLCOM] = eeglab;

% As this uses fastica algorithm for ICA, FastICA needs to be on the path, e.g.
addpath('C:/EEGLAB/FastICA_25')

% region of interest for finding maximal electrodes
roi = {'TP7', 'CP5', 'P7', 'P5', 'P9', 'PO7', 'PO3', 'O1'};

% cutoff probability for identifying eye and muscle related ICA components with ICLabel
icl_cutoff = 0.85;

% sigma parameter for ASR
asr_sigma = 20;

%% Clear Output Folders

delete(fullfile('localiser_sample_data', '*.csv'))

%% Import lab book

% handle commas in vectors
lab_book_file = fullfile('raw_data', 'stim-pc', 'participants.csv');
lab_book_raw_dat = fileread(lab_book_file);

[regstart, regend] = regexp(lab_book_raw_dat, '\[.*?\]');

for regmatch_i = 1:numel(regstart)
    str_i = lab_book_raw_dat(regstart(regmatch_i):regend(regmatch_i));
    str_i(str_i==',') = '.';
    lab_book_raw_dat(regstart(regmatch_i):regend(regmatch_i)) = str_i;
end

lab_book_fixed_file = fullfile('raw_data', 'stim-pc', 'participants_tmp.csv');
lab_book_fixed_conn = fopen(lab_book_fixed_file, 'w');
fprintf(lab_book_fixed_conn, lab_book_raw_dat);
fclose(lab_book_fixed_conn);

lab_book_readopts = detectImportOptions(lab_book_fixed_file, 'VariableNamesLine', 1, 'Delimiter', ',');
% read subject ids as class character
lab_book_readopts.VariableTypes{strcmp(lab_book_readopts.SelectedVariableNames, 'subj_id')} = 'char';
lab_book = readtable(lab_book_fixed_file, lab_book_readopts);

delete(lab_book_fixed_file)

%% Count the total number of excluded electrodes

n_bads = 0;
n_bads_per_s = zeros(size(lab_book, 1), 0);

for subject_nr = 1:size(lab_book, 1)
    bad_channels = eval(strrep(strrep(strrep(lab_book.loc_bad_channels{subject_nr}, '[', '{'), ']', '}'), '.', ','));
    n_bads_per_s(subject_nr) = numel(bad_channels);
    n_bads = n_bads + numel(bad_channels);
end

perc_bads = n_bads / (64 * size(lab_book, 1)) * 100;

%% Set up results table

max_elec_columns = {'subject_id',...
    'max_elec_bacs', 'max_time_bacs', 'max_diff_bacs',...
    'max_elec_noise', 'max_time_noise', 'max_diff_noise'};
empty_tablecells = cell(size(lab_book, 1), numel(max_elec_columns));
max_elecs = cell2table(empty_tablecells);
max_elecs.Properties.VariableNames = max_elec_columns;

%% Iterate over subjects

% record trial exclusions
total_excl_trials_incorr = zeros(1, size(lab_book, 1));
total_excl_trials_rt = zeros(1, size(lab_book, 1));

n_bad_ica = zeros(size(lab_book, 1), 0);

for subject_nr = 1:size(lab_book, 1)
    
    subject_id = lab_book.subj_id{subject_nr};
    fprintf('\n\n Subject Iteration %g/%g, ID: %s\n', subject_nr, size(lab_book, 1), subject_id)
    
    %% get subject-specific info from lab book
    exclude = lab_book.exclude(subject_nr);
    bad_channels = eval(strrep(strrep(strrep(lab_book.loc_bad_channels{subject_nr}, '[', '{'), ']', '}'), '.', ','));
    bad_channels_pictureword = eval(strrep(strrep(strrep(lab_book.pw_bad_channels{subject_nr}, '[', '{'), ']', '}'), '.', ','));
    bad_trigger_indices = eval(strrep(lab_book.loc_bad_trigger_indices{subject_nr}, '.', ','));

    % add PO4 to bad channels, which seems to be consistently noisy, even when not marked as bad
    if sum(strcmp('PO4', bad_channels))==0
        bad_channels(numel(bad_channels)+1) = {'PO4'};
    end
    
    %% abort if excluded
    
    if exclude
        % this is not planned to be used, but will be an easy way for other
        % researchers to see the effect of excluding specific participants
        % by editing the participants.csv file
        fprintf('Subject %s excluded. Preprocessing aborted.\n', subject_id)
        fprintf('Lab book note: %s\n', lab_book.note{subject_nr})
        continue
    end
    
    if (numel(bad_channels) >= 10) || (numel(bad_channels_pictureword) >= 10)
        fprintf('Subject %s excluded as >=10 electrodes marked as bad in either task. Preprocessing aborted.\n', subject_id)
        fprintf('Lab book note: %s\n', lab_book.note{subject_nr})
        continue
    end
    
    %% load participant's data
    
    % load raw eeg
    raw_datapath = fullfile(eeg_path, append(subject_id, '.bdf'));
    
    % abort if no EEG data collected yet
    if ~isfile(raw_datapath)
        fprintf('Subject %s skipped: no EEG data found\n', subject_id)
        continue
    end
    
    EEG = pop_biosig(raw_datapath, 'importevent', 'on', 'rmeventchan', 'off');
    
    % load behavioural
    all_beh_files = dir(beh_path);
    beh_regex_matches = regexpi({all_beh_files.name}, append('^', subject_id, '_.+\.csv$'), 'match');
    regex_emptymask = cellfun('isempty', beh_regex_matches);
    beh_regex_matches(regex_emptymask) = [];
    subj_beh_files = cellfun(@(x) x{:}, beh_regex_matches, 'UniformOutput', false);
    
    if numel(subj_beh_files)>1
        fprintf('%g behavioural files found?\n', size(subj_beh_files))
        break
    end
    
    beh_datapath = fullfile(beh_path, subj_beh_files{1});
    beh = readtable(beh_datapath);
    
    %% Set data features
    
    % set channel locations
    
    orig_locs = EEG.chanlocs;
    EEG.chanlocs = pop_chanedit(EEG.chanlocs, 'load', {'BioSemi64.loc', 'filetype', 'loc'});  % doesn't match order for the data
    
    % set channel types
    for ch_nr = 1:64
        EEG.chanlocs(ch_nr).type = 'EEG';
    end
    
    for ch_nr = 65:72
        EEG.chanlocs(ch_nr).type = 'EOG';
    end
    
    for ch_nr = 73:79
        EEG.chanlocs(ch_nr).type = 'MISC';
    end
    
    for ch_nr = 65:79
        EEG.chanlocs(ch_nr).theta = [];
        EEG.chanlocs(ch_nr).radius = [];
        EEG.chanlocs(ch_nr).sph_theta = [];
        EEG.chanlocs(ch_nr).sph_phi = [];
        EEG.chanlocs(ch_nr).X = [];
        EEG.chanlocs(ch_nr).Y = [];
        EEG.chanlocs(ch_nr).Z = [];
    end
    
    % change the order of channels in EEG.data to match the new order in chanlocs
    data_reordered = EEG.data;
    for ch_nr = 1:64        
        % make sure the new eeg data array matches the listed order
        ch_lab = EEG.chanlocs(ch_nr).labels;
        orig_locs_idx = find(strcmp(lower({orig_locs.labels}), lower(ch_lab)));
        data_reordered(ch_nr, :) = EEG.data(orig_locs_idx, :);
    end
    EEG.data = data_reordered;
    
    % remove unused channels
    EEG = pop_select(EEG, 'nochannel', 69:79);
    
    % plot the ROI for the paper
    if strcmp(subject_id, '1')

        roi_fig = figure;
        roi_idx = find(ismember({EEG.chanlocs.labels}, roi));

        hold on;
        topoplot(zeros(64, 0), EEG.chanlocs, 'electrodes', 'off');
        % set line width
        set(findall(gca, 'Type', 'Line'), 'LineWidth', 1);
        for i = 1:64
            if ismember(i, roi_idx)
                markcol = [1, 0, 0];
            else
                markcol = [0.75, 0.75, 0.75];
            end
            topoplot(zeros(64, 0), EEG.chanlocs, 'colormap', [0,0,0], 'emarker', {'.', markcol, 15, 1}, 'plotchans', i, 'headrad', 0);
        end
        hold off

        set(roi_fig, 'Units', 'Inches', 'Position', [0, 0, 1.5, 1.5], 'PaperUnits', 'Inches', 'PaperSize', [1.5, 1.5])
        exportgraphics(roi_fig, 'figs/roi_channels.pdf', 'BackgroundColor','none')

        close all

    end

    % remove bad channels
    ur_chanlocs = EEG.chanlocs;  % store a copy of the full channel locations before removing (for later interpolation)
    bad_channels_indices = find(ismember(lower({EEG.chanlocs.labels}), lower(bad_channels)));
    EEG = pop_select(EEG, 'nochannel', bad_channels_indices);
    
    %% Identify events (trials)
    
    % make the sopen function happy
    x = fileparts( which('sopen') );
    rmpath(x);
    addpath(x,'-begin');
    
    % build the events manually from the raw eeg file (pop_biosig removes event offsets)
    % NB: this assumes no resampling between reading the BDF file and now
    bdf_dat = sopen(raw_datapath, 'r', [0, Inf], 'OVERFLOWDETECTION:OFF');
    event_types = bdf_dat.BDF.Trigger.TYP;
    event_pos = bdf_dat.BDF.Trigger.POS;
    event_time = EEG.times(event_pos);
    sclose(bdf_dat);
    clear bdf_dat;
    
    triggers = struct(...
        'off', 0,...
        'word', 101,...
        'bacs', 102,...
        'noise', 103,...
        'practice', 25);
    
    % add 61440 to each trigger value (because of number of bits in pp)
    trigger_labels = fieldnames(triggers);
    for field_nr = 1:numel(trigger_labels)
        triggers.(trigger_labels{field_nr}) = triggers.(trigger_labels{field_nr}) + 61440;
    end
        
    % remove the first trigger if it is at time 0 and has a value which isn't a recognised trigger
    if (event_time(1)==0 && ~ismember(event_types(1), [triggers.off, triggers.word, triggers.bacs, triggers.noise, triggers.practice]))
        event_types(1) = [];
        event_pos(1) = [];
        event_time(1) = [];
    end
    
    % remove the new first trigger if it has a value of off
    if (event_types(1)==triggers.off)
        event_types(1) = [];
        event_pos(1) = [];
        event_time(1) = [];
    end
    
    % check every second trigger is an offset
    offset_locs = find(event_types==triggers.off);
    if any(offset_locs' ~= 2:2:numel(event_types))
        fprintf('Expected each second trigger to be an off?')
        break
    end
    
    % check every first trigger is non-zero
    onset_locs = find(event_types~=triggers.off);
    if any(onset_locs' ~= 1:2:numel(event_types))
        fprintf('Expected each first trigger to be an event?')
        break
    end
    
    % create the events struct manually    
    events_onset_types = event_types(onset_locs);
    events_onsets = event_pos(onset_locs);
    events_offsets = event_pos(offset_locs);
    events_durations = events_offsets - events_onsets;
    
    EEG.event = struct();
    for event_nr = 1:numel(events_onsets)
        EEG.event(event_nr).type = events_onset_types(event_nr);
        EEG.event(event_nr).latency = events_onsets(event_nr);
        EEG.event(event_nr).offset = events_offsets(event_nr);
        EEG.event(event_nr).duration = events_durations(event_nr);
    end
    
    % copy the details over to urevent
    EEG.urevent = EEG.event;
    
    % record the urevent in event, for reference if they change
    for event_nr = 1:numel(events_onsets)
        EEG.event(event_nr).urevent = event_nr;
    end
    
    % remove bad events recorded in lab book (misfired triggers)
    EEG = pop_editeventvals(EEG, 'delete', find(ismember([EEG.event.urevent], bad_trigger_indices)));
    
    % remove practice trials
    EEG = pop_editeventvals(EEG, 'delete', find(ismember([EEG.event.type], triggers.practice)));
    
    % check the events make sense
    if sum(~ismember([EEG.event.type], [triggers.word, triggers.bacs, triggers.noise])) > 0
        fprintf('Unexpected trial types?\n')
        break
    end
    
    if numel({EEG.event.type})~=300
        fprintf('%g trial triggers detected (expected 300)?\n',  numel({EEG.event.type}))
        break
    end
    
    if sum(ismember([EEG.event.type], [triggers.word])) ~= sum(ismember([EEG.event.type], [triggers.bacs]))
        fprintf('Unequal number of word and BACS trials?\n')
        break
    end
    
    if sum(ismember([EEG.event.type], [triggers.word])) ~= sum(ismember([EEG.event.type], [triggers.noise]))
        fprintf('Unequal number of word and noise trials?\n')
        break
    end
    
    % add the trials' onsets, offsets, durations, and triggers to the behavioural data
    beh.event = zeros(size(beh, 1), 1);
    beh.latency = zeros(size(beh, 1), 1);
    for row_nr = 1:size(beh, 1)
        cond_i = beh.condition(row_nr);
        beh.event(row_nr) = triggers.(cond_i{:});  % look up the trial's expected trigger
        beh.latency(row_nr) = EEG.event(row_nr).latency;
        beh.offset(row_nr) = EEG.event(row_nr).offset;
        beh.duration(row_nr) = EEG.event(row_nr).duration;
        beh.duration_ms(row_nr) = (EEG.event(row_nr).duration * 1000/EEG.srate) - 500;  % minus 500 as event timer starts at word presentation, but rt timer starts once word turns green
    end
    
    % check events expected in beh are same as those in the events struct
    if any(beh.event' ~= [EEG.event.type])
        fprintf('%g mismatches between behavioural data and triggers?\n', sum(beh.event' ~= [EEG.event.type]))
        break
    end
    
    % check the difference between the durations and the response times (should be very small)
    % hist(beh.rt - beh.duration_ms, 100)
    
    % record trial numbers in EEG.event
    for row_nr = 1:size(beh, 1)
        EEG.event(row_nr).trl_nr = beh.trl_nr(row_nr);
    end
    
    %% Remove segments of data that fall outside of blocks
    
    % record block starts
    beh.is_block_start(1) = 1;
    for row_nr = 2:size(beh, 1)
        beh.is_block_start(row_nr) = beh.block_nr(row_nr) - beh.block_nr(row_nr-1) == 1;
    end
    % record block ends
    beh.is_block_end(size(beh, 1)) = 1;
    for row_nr = 1:(size(beh, 1)-1)
        beh.is_block_end(row_nr) = beh.block_nr(row_nr+1) - beh.block_nr(row_nr) == 1;
    end
    % record block boundaries (first start and last end point of each block, with 0.75 seconds buffer)
    beh.block_boundary = zeros(size(beh, 1), 1);
    for row_nr = 1:size(beh, 1)
        if beh.is_block_start(row_nr)
            beh.block_boundary(row_nr) = beh.latency(row_nr) - (EEG.srate * 0.75);
        elseif beh.is_block_end(row_nr)
            beh.block_boundary(row_nr) = beh.offset(row_nr) + (EEG.srate * 0.75);
        end
    end
    
    % get the boundary indices in required format (start1, end1; start2, end2; start3, end3)
    block_boundaries = reshape(beh.block_boundary(beh.block_boundary~=0), 2, [])';
    
    % remove anything outside of blocks
    EEG = pop_select(EEG, 'time', (block_boundaries / EEG.srate));
    
    %% Trial selection
    
    % include only correct responses
    beh_filt_acc_only = beh(beh.acc==1, :);
    excl_trials_incorr = size(beh, 1)-size(beh_filt_acc_only, 1);
    total_excl_trials_incorr(subject_nr) = excl_trials_incorr;
    fprintf('Lost %g trials to incorrect responses\n', excl_trials_incorr)
    
    % include only responses faster than 1500 ms
    beh_filt = beh_filt_acc_only(beh_filt_acc_only.rt<=1500, :);
    excl_trials_rt = size(beh_filt_acc_only, 1)-size(beh_filt, 1);
    total_excl_trials_rt(subject_nr) = excl_trials_rt;
    fprintf('Lost %g trials to RTs above 1500\n', excl_trials_rt)
    
    fprintf('Lost %g trials in total to behavioural data\n', size(beh, 1)-size(beh_filt, 1))
    
    % filter the events structure
    discarded_trls = beh.trl_nr(~ismember(beh.trl_nr, beh_filt.trl_nr));
    discarded_events_indices = [];  % (collect in a for loop, as [EEG.event.trl_nr] would remove missing data)
    for event_nr = 1:size(EEG.event, 2)
        if ismember(EEG.event(event_nr).trl_nr, discarded_trls)
            discarded_events_indices = [discarded_events_indices, event_nr];
        end
    end
    EEG = pop_editeventvals(EEG, 'delete', discarded_events_indices);
    
    % check the discarded trials are the expected length
    if numel(discarded_trls) ~= size(beh, 1)-size(beh_filt, 1)
        fprintf('Mismatch between behavioural data and EEG events in the number of trials to discard?')
        break
    end
    
    % check the sizes match
    if numel([EEG.event.trl_nr]) ~= size(beh_filt, 1)
        fprintf('Inconsistent numbers of trials between events structure and behavioural data after discarding trials?')
        break
    end
    
    % check the trl numbers match
    if any([EEG.event.trl_nr]' ~= beh_filt.trl_nr)
        fprintf('Trial IDs mmismatch between events structure and behavioural data after discarding trials?')
        break
    end
    
    %% Rereference, downsample, and filter
    
    % rereference
    EEG = pop_reref(EEG, []);
    
    % downsample if necessary
    if EEG.srate ~= 512
        EEG = pop_resample(EEG, 512);
    end
    
    % filter
    % EEG = eeglab_butterworth(EEG, 0.5, 40, 4, 1:size(EEG.chanlocs, 2));  % preregistered filter
    EEG = eeglab_butterworth(EEG, 0.1, 40, 4, 1:size(EEG.chanlocs, 2));  % filter with lower highpass

    %% ICA
    
    % apply ASR
    %EEG_no_asr = EEG;
    EEG = clean_asr(EEG, asr_sigma, [], [], [], [], [], [], [], [], 1024);  % The last number is available memory in mb, needed for reproducibility

    rng(3101)  % set seed for reproducibility
    EEG = pop_runica(EEG, 'icatype', 'fastica', 'approach', 'symm');

    % classify components with ICLabel
    EEG = iclabel(EEG);

    % store results for easy indexing
    icl_res = EEG.etc.ic_classification.ICLabel.classifications;
    icl_classes = EEG.etc.ic_classification.ICLabel.classes;
    
    % identify and remove artefact components
    artefact_comps = find(icl_res(:, strcmp(icl_classes, 'Eye')) >= icl_cutoff | icl_res(:, strcmp(icl_classes, 'Muscle')) >= icl_cutoff);
    fprintf('Removing %g artefact-related ICA components\n', numel(artefact_comps))
    n_bad_ica(subject_nr) = numel(artefact_comps);
    %EEG_no_iclabel = EEG;
    EEG = pop_subcomp(EEG, artefact_comps);
            
    %% Interpolate bad channels
    
    % give the original chanlocs structure so EEGLAB interpolates the missing electrode(s)
    if numel(bad_channels)>0
        EEG = pop_interp(EEG, ur_chanlocs);
    end
    
    %% Epoch the data
    
    % identify and separate into epochs
    EEG_epo = struct();
    EEG_epo.word = pop_epoch(EEG, {triggers.word}, [-0.25, 1]);
    EEG_epo.bacs = pop_epoch(EEG, {triggers.bacs}, [-0.25, 1]);
    EEG_epo.noise = pop_epoch(EEG, {triggers.noise}, [-0.25, 1]);
    
    % remove baseline
    EEG_epo.word = pop_rmbase(EEG_epo.word, [-200, 0]);
    EEG_epo.bacs = pop_rmbase(EEG_epo.bacs, [-200, 0]);
    EEG_epo.noise = pop_rmbase(EEG_epo.noise, [-200, 0]);
    
    % check times vectors are identical
    if ~isequal(EEG_epo.word.times, EEG_epo.bacs.times, EEG_epo.noise.times)
        fprintf('The times vectors in the epoch structures are not identical!')
        break
    end
    
    %% Get the maximal electrode
    % (word Vs. BACS for main analysis, but word Vs. noise also found)
    
    fprintf('Getting maximal electrodes...\n')
    
    % get channel means for each condition
    ch_avg = struct();
    ch_avg.word = mean(EEG_epo.word.data, 3);
    ch_avg.bacs = mean(EEG_epo.bacs.data, 3);
    ch_avg.noise = mean(EEG_epo.noise.data, 3);
    
    % get index of time window
    targ_window = [120, 200];
    targ_window_idx = EEG_epo.word.times >= targ_window(1) & EEG_epo.word.times <= targ_window(2);
    
    % get index of roi channels
    eeg_chan_idx = ismember({EEG.chanlocs.labels}, roi);  % EEG chanlocs same as chanlocs in ch_avg structs as they're copied over
    
    % store vectors of times and channels in ch_avg
    ch_avg.times = EEG_epo.word.times(targ_window_idx);  % taken from word condition but identical across conditions
    ch_avg.chanlocs = EEG.chanlocs(eeg_chan_idx);
    
    % get only roi electrode data in target window
    ch_avg.word = ch_avg.word(eeg_chan_idx, targ_window_idx);
    ch_avg.bacs = ch_avg.bacs(eeg_chan_idx, targ_window_idx);
    ch_avg.noise = ch_avg.noise(eeg_chan_idx, targ_window_idx);
    
    % get differences of interest
    % - directional, so find max of these
    ch_avg.diff_word_bacs = ch_avg.bacs - ch_avg.word;
    ch_avg.diff_word_noise = ch_avg.noise - ch_avg.word;
    
    % find the maximum difference indices
    mean_bacs_diff_perchan = mean(ch_avg.diff_word_bacs, 2);
    max_bacs_ch_idx = mean_bacs_diff_perchan == max(mean_bacs_diff_perchan);
    
    mean_noise_diff_perchan = mean(ch_avg.diff_word_noise, 2);
    max_noise_ch_idx = mean_noise_diff_perchan == max(mean_noise_diff_perchan);
    
    % if multiple channels have an equal mean difference, select one randomly (but reprocubily)
    if sum(max_bacs_ch_idx) > 1
        rng(42 + subject_nr)
        perm_idx = randperm(sum(max_bacs_ch_idx));
        maxes_idx = find(max_bacs_ch_idx);
        max_bacs_ch_idx = maxes_idx(perm_idx(numel(maxes_idx)));
    end
    
    if sum(max_noise_ch_idx) > 1
        rng(42 + subject_nr)
        perm_idx = randperm(sum(max_noise_ch_idx));
        maxes_idx = find(max_noise_ch_idx);
        max_noise_ch_idx = maxes_idx(perm_idx(numel(maxes_idx)));
    end
    
    % get the channel names
    chan_names = {ch_avg.chanlocs.labels};
    max_chan_bacs = chan_names{max_bacs_ch_idx};
    max_chan_noise = chan_names{max_noise_ch_idx};
    
    % get the timepoint and value of the maximum difference for the max channels
    [max_chan_bacs_peak_diff, max_chan_bacs_peak_diff_idx] = max(ch_avg.diff_word_bacs(max_bacs_ch_idx, :));
    max_chan_bacs_peak_diff_signed = ch_avg.diff_word_bacs(max_bacs_ch_idx, max_chan_bacs_peak_diff_idx);
    max_chan_bacs_peak_time = ch_avg.times(max_chan_bacs_peak_diff_idx);
    
    [max_chan_noise_peak_diff, max_chan_noise_peak_diff_idx] = max(ch_avg.diff_word_noise(max_noise_ch_idx, :));
    max_chan_noise_peak_diff_signed = ch_avg.diff_word_noise(max_noise_ch_idx, max_chan_noise_peak_diff_idx);
    max_chan_noise_peak_time = ch_avg.times(max_chan_noise_peak_diff_idx);
    
    % store the values in the table
    max_elecs.subject_id(subject_nr) = {subject_id};
    
    max_elecs.max_elec_bacs(subject_nr) = {max_chan_bacs};
    max_elecs.max_time_bacs(subject_nr) = {max_chan_bacs_peak_time};
    max_elecs.max_diff_bacs(subject_nr) = {max_chan_bacs_peak_diff_signed};
    
    max_elecs.max_elec_noise(subject_nr) = {max_chan_noise};
    max_elecs.max_time_noise(subject_nr) = {max_chan_noise_peak_time};
    max_elecs.max_diff_noise(subject_nr) = {max_chan_noise_peak_diff_signed};
    
    %% Save sample-level data for all electrodes
    
    disp('Getting sample-level localiser results...')
    
    % resample to 256 Hz
    EEG_256 = pop_resample(EEG, 256);
    
    % get epochs of low-srate data
    EEG_epo_256 = pop_epoch(EEG_256, {triggers.word, triggers.bacs, triggers.noise}, [-0.25, 0.5]);
    % remove baseline
    EEG_epo_256 = pop_rmbase(EEG_epo_256, [-200, 0]);
    
    % pre-allocate the table
    var_names = {'subj_id', 'stim_grp', 'resp_grp', 'trl_nr', 'ch_name', 'time', 'uV'};
    var_types = {'string', 'string', 'string', 'double', 'string', 'double', 'double'};
    nrows = 64 * size(EEG_epo_256.times, 2) * size(beh_filt, 1);
    sample_res = table('Size',[nrows, numel(var_names)], 'VariableTypes',var_types, 'VariableNames',var_names);
    
    sample_res.subj_id = repmat(beh_filt.subj_id, 64*size(EEG_epo_256.times, 2), 1);
    sample_res.stim_grp = repmat(beh_filt.stim_grp, 64*size(EEG_epo_256.times, 2), 1);
    sample_res.resp_grp = repmat(beh_filt.resp_grp, 64*size(EEG_epo_256.times, 2), 1);
    
    % get the 64 channel eeg data as an array
    eeg_arr = EEG_epo_256.data(1:64, :, :);
    
    % a vector of all eeg data
    eeg_vec = squeeze(reshape(eeg_arr, 1, 1, []));
    
    % array and vector of the channel labels for each value in EEG.data
    channel_labels_arr = cell(size(eeg_arr));
    channel_label_lookup = {EEG_epo_256.chanlocs.labels};
    for chan_nr = 1:size(eeg_arr, 1)
        channel_labels_arr(chan_nr, :, :) = repmat(channel_label_lookup(chan_nr), size(channel_labels_arr, 2), size(channel_labels_arr, 3));
    end
    
    channel_labels_vec = squeeze(reshape(channel_labels_arr, 1, 1, []));
    
    % array and vector of the item numbers for each value in EEG.data
    times_arr = zeros(size(eeg_arr));
    times_lookup = EEG_epo_256.times;
    for time_idx = 1:size(eeg_arr, 2)
        times_arr(:, time_idx, :) = repmat(times_lookup(time_idx), size(times_arr, 1), size(times_arr, 3));
    end
    
    times_vec = squeeze(reshape(times_arr, 1, 1, []));
    
    % array and vector of the trial numbers
    trials_arr = zeros(size(eeg_arr));
    trials_lookup = beh_filt.trl_nr;
    for trl_idx = 1:size(eeg_arr, 3)
        trials_arr(:, :, trl_idx) = repmat(trials_lookup(trl_idx), size(trials_arr, 1), size(trials_arr, 2));
    end
    
    trials_vec = squeeze(reshape(trials_arr, 1, 1, []));
        
    % store sample-level results in the table
    sample_res.ch_name = channel_labels_vec;
    sample_res.trl_nr = trials_vec;
    sample_res.time = times_vec;
    sample_res.uV = eeg_vec;
    
    % look up and store some info about the trials
    trial_info_lookup = beh_filt(:, {'trl_nr', 'condition', 'string', 'item_nr'});
    sample_res = outerjoin(sample_res, trial_info_lookup, 'MergeKeys', true);
    
    % sort by time, channel, item_nr
    sample_res = sortrows(sample_res, {'time', 'ch_name', 'trl_nr'});
    
    % Save the sample-level results
    disp('Saving sample-level localiser results...')
    writetable(sample_res, fullfile('localiser_sample_data', [subject_id, '.csv']));
    
end

%% save the results

fprintf('\nSaving results...\n')
writetable(max_elecs, 'max_elecs.csv');

fprintf('Finished preprocessing localiser data!\n')

%% Functions

% custom function for applying a Butterworth filter to EEGLAB data
function EEG = eeglab_butterworth(EEG, low, high, order, chanind)
    fprintf('Applying Butterworth filter between %g and %g Hz (order of %g)\n', low, high, order)
    % create filter
    [b, a] = butter(order, [low, high]/(EEG.srate/2));
    % apply to data (requires transposition for filtfilt)
    data_trans = single(filtfilt(b, a, double(EEG.data(chanind, :)')));
    EEG.data(chanind, :) = data_trans';
end

% custom function for finding the closest timepoint in an EEG dataset
function [idx, closesttime] = eeglab_closest_time(EEG, time)
    dists = abs(EEG.times - time);
    idx = find(dists == min(dists));
    % in the unlikely case there are two equidistant times, select one randomly
    if numel(idx) > 1
        fprintf('Two equidistant times! Selecting one randomly.')
        idx = idx(randperm(numel(idx)));
        idx = idx(1);
    end
    closesttime = EEG.times(idx);
end