% Introduction
% This script includes all steps of ERP data analysis. It helps with both batch running multiple subjects as well as documenting steps of data analysis 
%for later analysis. It runs faster withtout the matlab GUI. To run this in the terminal use "matlab -nodisplay < 109_AnalysisScript.m > 109log". You 
%need to keep a log file and inspect it to ensure that there were no errors while running the script. The log also has information about artifact 
%rejection percentages. The script was created by Firat Soylu for the ELDEN Lab (elden.ua.edu) at the University of Alabama, on October, 2016.

% Start
% Clear memory and the command window.
clear
clc
%  You might need to run EEGLAB and to get some paths in place before you run the script

% Script Settings
% Set 1 to enable, 0 to disable

clean_up      = 0;    % Cleans up all previous analysis files

import_data   = 0;    % Import raw data files and resample

run_ica       = 0;    % Run ICA, manual input for component rejection is needed after ICA

remove_ica    = 0;    % After run ICA, inspect data,d determine blink components to remove enter it in the list ICA_to_remove

process_data  = 1;    % Filtering, event list, binlister, epoching, artifact rejection, average, lp filter, diff waves

grand_avg     = 1;    % Calculate grant average and grand plots

save_files    = 1;    % Save files for each step. Essential files are saved regardless

plot_figures  = 0;    % Set to 1 to plot figures for each subject

% Subjects
subject_list = {'1160','1161','1163','1164','1165','1168','1169','1170','1172','1174','1175','1176','1177','1178','1179','1180','1181','1182','1183'};


% Channels to Interpolate
% Plot the data and determine which channels are noisy and needs to be interpolated for each subject
%chan_interpolate = {{},{},{},{'ch10'},{}};
%chan_interpolate = {{'ch10'}};
chan_interpolate = {{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{},{}};

% ICA's to Remove
%After run ICA, inspect data, determine blink components to remove enter it in the list ICA_to_remove
ICA_to_remove = {4, 3, 3, 4}; %Enter values here as well, showing subj/component match. BR_Test->4, 1160-> 3, 1161-> 3, 1163-> 4

% Global Variables
% Path to the parent folder, which contains the data folders for all subjects
home_path  = '/Path to main analysis folder/';
ALLERP = buildERPstruct([]); % Initialize the ALLERP structure and CURRENTERP
CURRENTERP = 0;
nsubj = length(subject_list); % number of subjects
figScale = [ -200.0 500.0   -200:10:500 ]; % Figure intervals

% Clean-Up
% This deletes all previous analysis files
if (clean_up)
    for s=1:nsubj % Loop through all subjects
        data_path  = [home_path subject_list{s} '/'];
        to_delete = {'*.txt', '*.set', '*.fdt', '*.erp'};
        for ext=1:length(to_delete)
            delete([data_path to_delete{ext}]);
            delete([data_path 'figures/*']);
        end
    end

fprintf('--------Deleted all previous analysis files--------')
end

% Import Data
if (import_data)
    for s=1:nsubj % Loop through all subjects
        fprintf('\n******\nProcessing subject %s\n******\n\n', subject_list{s});

        % Path to the folder containing the current subject's data
        data_path  = [home_path subject_list{s} '/'];
        disp([data_path 'rawdata/' subject_list{s} '.vhdr']);
        % import raw data
        EEG = pop_loadbv([data_path 'rawdata/'], [subject_list{s} '.vhdr']);
        %EEG = pop_loadbv([data_path 'RawData/'], [subject_list{s} '.vhdr'],[20000 50000]); %shorter segment of data for test purposes
        EEG.setname=subject_list{s};
        if (save_files)
            EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
        end

        % resample to 250
        EEG = pop_resample( EEG, 250);
        EEG.setname=[subject_list{s} '_rsmpld'];
        EEG = eeg_checkset( EEG );
        if (save_files)
            EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
        end
        
        % interpolate noisy electrodes
        chan_to_inter = chan_interpolate{s};
        if size(chan_to_inter) > 0
            for n=1:size(chan_to_inter)
                EEG = pop_eegchanoperator( EEG, { [chan_to_inter{n} '= chinterpol']} , 'ErrorMsg', 'popup', 'Warning', 'on' );
            end
        end
        EEG.setname=[EEG.setname '_intpl'];
        if (save_files)
            EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
        end
        
        %re-reference to average
        EEG = pop_chanedit(EEG, 'append',31,'changefield',{32 'labels' 'Cz'},'lookup',[home_path '../standard_BESA/standard-10-5-cap385.elp'],'setref',{'1:31' 'Cz'});
        %[ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET);
        EEG = pop_reref( EEG, [],'refloc',struct('labels',{'Cz'},'type',{''},'theta',{0},'radius',{0},'X',{5.2047e-15},'Y',{0},'Z',{85},'sph_theta',{0},'sph_phi',{90},'sph_radius',{85},'urchan',{32},'ref',{''},'datachan',{0}));
        EEG.setname=[EEG.setname  '_reref'];
        if (save_files)
            EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
        end
        
        % HP filter at 0.1
        EEG  = pop_basicfilter( EEG,  1:32 , 'Boundary', 'boundary', 'Cutoff', 0.1, 'Design', 'butter', 'Filter', 'highpass', 'Order',  2, 'RemoveDC', 'on' ); 
        EEG.setname = [EEG.setname '_hpfilt'];
        EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
        
    end
end

% Run ICA
    % you can run ICA here, getica.m removeica.m
if (run_ica)
    for s=1:nsubj % Loop through all subjects
        data_path  = [home_path subject_list{s} '/'];   % Path to the folder containing the current subject's data
        EEG = pop_loadset('filename',[subject_list{s} '_rsmpld_intpl_reref_hpfilt.set'], 'filepath', data_path);
        EEG = pop_runica(EEG, 'icatype','runica', 'dataset',1,'options',{'extended' 1}); %Run ICA, manual component rejection is required after this
        EEG.setname = [EEG.setname '_ica'];
        EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
    end
end
 
% Remove ICA
% Remove visually inspected ICA blink components.
if (remove_ica)
    for s=1:nsubj % Loop through all subjects
        data_path  = [home_path subject_list{s} '/'];   % Path to the folder containing the current subject's data
        EEG = pop_loadset('filename',[subject_list{s} '_rsmpld_intpl_reref_hpfilt_ica.set'], 'filepath', data_path);
        EEG = pop_subcomp( EEG, ICA_to_remove{s}, 0);
        EEG.setname = [subject_list{s} '_rsmpld_intpl_reref_hpfilt_ica_icapruned'];
        EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
    end
end
% Process Data

if (process_data)
    for s=1:nsubj % Loop through all subjects
        data_path  = [home_path subject_list{s} '/'];   % Path to the folder containing the current subject's data
        
        
        if exist([data_path subject_list{s} '_rsmpld_intpl_reref_hpfilt_elist_bins_be_ica_icapruned.set'],'file')
             EEG = pop_loadset('filename',[subject_list{s} '_rsmpld_intpl_reref_hpfilt_ica_icapruned.set'], 'filepath', data_path);
        else
             EEG = pop_loadset('filename',[subject_list{s} '_rsmpld_intpl_reref_hpfilt.set'], 'filepath', data_path);
        end
        
        % create event list
        EEG  = pop_creabasiceventlist( EEG , 'AlphanumericCleaning', 'on', 'BoundaryNumeric', { -99 }, 'BoundaryString', { 'boundary' }, 'Eventlist', [data_path 'elist.txt'] );
        EEG.setname = [EEG.setname '_elist'];
        if (save_files)
            EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
        end

        % assign events to bins with binlister
        EEG  = pop_binlister( EEG , 'BDF', [home_path 'BinFiles/' 'binDescriptor.txt'], 'ExportEL', [data_path 'elist.txt'], 'Ignore',  246, 'IndexEL',  1, 'SendEL2', 'EEG&Text', 'Voutput', 'EEG' );
        EEG.setname = [EEG.setname '_bins'];
        if (save_files)
            EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
        end

        % epoch data
        EEG = pop_epochbin( EEG , [-200.0  500.0],  'pre');
        EEG.setname = [EEG.setname '_be'];
        if (save_files)
            EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
        end

        % detect & reject artifacts with simple voltage threshold, Fp1 & Fp2
        EEG  = pop_artextval( EEG , 'Channel',  [1, 30], 'Flag',  1, 'Threshold', [ -100 100], 'Twindow', [ -200 496] );
        EEG.setname = [EEG.setname '_ar'];
        if (save_files)
            EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
        end
        
        % report percentage of rejected trials (collapsed across all bins)
        artifact_proportion = getardetection(EEG);
        fprintf('%s: Percentage of rejected trials was %1.2f\n', subject_list{s}, artifact_proportion);
        
        % average epochs & save
        fprintf('\n\n\n**** %s: Averaging ****\n\n\n', subject_list{s});  
        ERP = pop_averager( EEG, 'Criterion', 'good', 'SEM', 'on');
        ERP.erpname = [subject_list{s} '_ERPs'];  % name for erpset menu
        pop_savemyerp(ERP, 'erpname', ERP.erpname, 'filename', [ERP.erpname '.erp'], 'filepath', data_path);
        if (save_files)
            EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);
        end

        % create difference waves
        ERP = pop_binoperator( ERP, [home_path 'BinFiles/binEquations.txt']);
        ERP.erpname = [ERP.erpname '_diff'];  % name for erpset menu
        if (save_files)
            pop_savemyerp(ERP, 'erpname', ERP.erpname, 'filename', [ERP.erpname '.erp'], 'filepath', data_path);
        end

        % LP filter at 30 
        EEG  = pop_basicfilter( EEG,  1:32 , 'Boundary', 'boundary', 'Cutoff', 30, 'Design', 'butter', 'Filter', 'lowpass', 'Order',  2, 'RemoveDC', 'on' ); 
        EEG.setname = [EEG.setname '_lpfilt'];
        EEG = pop_saveset(EEG, 'filename', [EEG.setname '.set'], 'filepath', data_path);

        % Measure the mean amplitude for 100-150ms & 200-300ms in bins 4,5,6,
        % in channels right: 11,12,14,15, left: 23,22,19,20 & save

       % values1 = pop_geterpvalues( ERP, [100 150], 4:6, [11 12 14 15 23 22 19 20], 'Baseline', 'pre', 'Filename', [data_path 'measures1_100-150.txt'], 'Measure', 'meanbl', 'Resolution',2 );
       % values2 = pop_geterpvalues( ERP, [200 300], 4:6, [11 12 14 15 23 22 19 20], 'Baseline', 'pre', 'Filename', [data_path 'measures2_200-300.txt'], 'Measure', 'meanbl', 'Resolution',2 );

        CURRENTERP = CURRENTERP + 1;
        ALLERP(CURRENTERP) = ERP;
        if (plot_figures)
             ERP = pop_ploterps( ERP,  1:3, 1:32 , 'AutoYlim', 'on', 'Axsize', [ 0.15 0.24], 'BinNum', 'on', 'Blc', 'pre', 'Box', [ 6 6], 'ChLabel', 'on',...
            'FontSizeChan',  10, 'FontSizeLeg',  12, 'FontSizeTicks',  10, 'LegPos', 'bottom', 'Linespec', {'k-' , 'r-' , 'b-' }, 'LineWidth',  1,...
            'Maximize', 'off', 'Position', [ 103.667 29.5833 107 32], 'Style', 'Topo', 'Tag', 'ERP_figure', 'Transparency',  0, 'xscale', figScale,...
            'YDir', 'normal' );
            saveas(gcf, [data_path 'figures/' subject_list{s} '_ERP.png'],'png');
            saveas(gcf, [data_path 'figures/' subject_list{s} '_ERP.fig'],'fig');
        end
    end
end

% Grand Averages
if (grand_avg)
    ALLERP = buildERPstruct([]);
    CURRENTERP = 0;
    for s=1:nsubj % Loop through all subjects to add ERPs for each subject to ALLERP
        ERP = pop_loaderp( 'filename', [subject_list{s} '_ERPs.erp'], 'filepath', [home_path subject_list{s} '/']);
        CURRENTERP = CURRENTERP + 1;
        ALLERP(CURRENTERP) = ERP;
    end
    
    % Make a grand average using ALLERP
    ERP = pop_gaverager( ALLERP , 'Criterion', 100, 'ERPindex', 1:nsubj ); %'Criterion', 100 is left over from a previous version
    ERP.erpname = 'grand_avg';  % name for erpset menu
    ERP = pop_savemyerp(ERP, 'filename', [ERP.erpname '.erp'], 'filepath', [home_path 'GrandAvgERPs/']);
    CURRENTERP = CURRENTERP + 1;
    ALLERP(CURRENTERP) = ERP;
    
    % Create a LP 30 filtered version for plot
    ERP = pop_filterp( ERP,1:32 , 'Cutoff',30, 'Design', 'butter', 'Filter', 'lowpass');
    ERP.erpname = [ERP.erpname '_30Hz'];  % name for erpset menu
    ERP = pop_savemyerp(ERP, 'filename', [ERP.erpname '.erp'], 'filepath', [home_path 'GrandAvgERPs/']);
    CURRENTERP = CURRENTERP + 1;
    ALLERP(CURRENTERP) = ERP;

    % Plot grand figures
    ERP = pop_ploterps( ERP,  1:3, 1:32 , 'AutoYlim', 'on', 'Axsize', [ 0.15 0.24], 'BinNum', 'on', 'Blc', 'pre', 'Box', [ 6 6], 'ChLabel', 'on',...
    'FontSizeChan',  10, 'FontSizeLeg',  12, 'FontSizeTicks',  10, 'LegPos', 'bottom', 'Linespec', {'k-' , 'r-' , 'b-' }, 'LineWidth',  1,...
    'Maximize', 'off', 'Position', [ 103.667 29.5833 107 32], 'Style', 'Topo', 'Tag', 'ERP_figure', 'Transparency',  0, 'xscale', figScale,...
    'YDir', 'normal' );
    saveas(gcf, [home_path 'AveragedFigures/grandAvgERP.fig'],'fig');
    
    ERP = pop_ploterps( ERP,  4:7, 1:32 , 'AutoYlim', 'on', 'Axsize', [ 0.15 0.24], 'BinNum', 'on', 'Blc', 'pre', 'Box', [ 6 6], 'ChLabel', 'on',...
    'FontSizeChan',  10, 'FontSizeLeg',  12, 'FontSizeTicks',  10, 'LegPos', 'bottom', 'Linespec', {'k-' , 'r-' , 'b-' }, 'LineWidth',  1,...
    'Maximize', 'off', 'Position', [ 103.667 29.5833 107 32], 'Style', 'Topo', 'Tag', 'ERP_figure', 'Transparency',  0, 'xscale', figScale,...
    'YDir', 'normal' );
    saveas(gcf, [home_path 'AveragedFigures/grandAvgDiffERP.fig'],'fig');

    % Measure the mean amplitude from specific intervals and electrodes and
    % save it in a measures file
    %values1 = pop_geterpvalues( ALLERP, [100 150], 4:6, [11 12 14 15 23 22 19 20], 'Baseline', 'pre', 'Erpsets', [1:nsubj], 'Filename', [home_path 'Measures/' 'measures1_100-150.txt'], 'Measure', 'meanbl', 'Resolution',2 );
    %values2 = pop_geterpvalues( ALLERP, [200 300], 4:6, [11 12 14 15 23 22 19 20], 'Baseline', 'pre', 'Erpsets', [1:nsubj], 'Filename', [home_path 'Measures/' 'measures2_200-300.txt'], 'Measure', 'meanbl', 'Resolution',2 );
end
fprintf('\n\n\n**** FINISHED ****\n\n\n');
