%% Extract peaks
clear all
close all

main_dir = fileparts(mfilename('fullpath')); % path of the script folder
cd(main_dir)
load('subjects.mat')
Paths_file


gen_dir = ''; %directory where TEP data is present


for iSub=1:28
    for block = 1:7
        subjdir_filedef = [gen_dir subjects{iSub} '\' subjects{iSub} '_B_' sprintf('%01d', block) '\' '6_final_preprocessing_new_interp\'];
        cd(subjdir_filedef)
        filedefdir = dir(fullfile('*def.mat'));
        filename = filedefdir.name;
        load(filename)
        if contains(filename, 'Mo_PA')
            Mo_PA{iSub} = TEP_final2;
        elseif contains(filename, 'Mo_AP')
            Mo_AP{iSub} = TEP_final2;
        elseif contains(filename, 'Mo_LM')
            Mo_LM{iSub} = TEP_final2;
        elseif contains(filename, 'Bi_PA_Contr')
            Bi_PA_Contr{iSub} = TEP_final2;
        elseif contains(filename, 'Bi_AP')
            Bi_AP{iSub} = TEP_final2;
        elseif contains(filename, 'Bi_LM')
            Bi_LM{iSub} = TEP_final2;
        else
            Bi_PA{iSub} = TEP_final2;
        end
        clear TEP_final2
    end
end

allcond = {'Bi_AP', 'Bi_LM', 'Bi_PA', 'Bi_PA_Contr', 'Mo_AP', 'Mo_LM', 'Mo_PA'};

for c = 1:7
    cond = allcond{c};
    
    if contains(cond, 'Mo_PA')
        TEP_allsubj = Mo_PA;
    elseif contains(cond, 'Mo_AP')
        TEP_allsubj = Mo_AP;
    elseif contains(cond, 'Mo_LM')
        TEP_allsubj = Mo_LM;
    elseif contains(cond, 'Bi_PA_Contr')
        TEP_allsubj = Bi_PA_Contr;
    elseif contains(cond, 'Bi_AP')
        TEP_allsubj = Bi_AP;
    elseif contains(cond, 'Bi_LM')
        TEP_allsubj = Bi_LM;
    else
        TEP_allsubj = Bi_PA;
    end
    
    cfg=[];
    TEP_gavg = ft_timelockgrandaverage(cfg, TEP_allsubj{:})
    cfg.keepindividual='yes';
    TEP_gavg_ind = ft_timelockgrandaverage(cfg, TEP_allsubj{:})
    
    cfg = [];
    cfg.channel={'F4' 'FC4'};
    cfg.avgoverchan = 'yes';
    P15 = ft_selectdata(cfg, TEP_gavg_ind);
    
    GAVG_all_peak=TEP_gavg_ind;
    GAVG_all_peak.individual=cat(2,TEP_gavg_ind.individual, ...
        zeros(length(TEP_gavg_ind.individual(:,1,1)), 1, length(TEP_gavg_ind.individual(1,1,:))), ...
        zeros(length(TEP_gavg_ind.individual(:,1,1)), 1, length(TEP_gavg_ind.individual(1,1,:)))) %one extrachannel for the marker and one for the average
    GAVG_all_peak.individual(:, (length(TEP_gavg_ind.label)+1), :) = P15.individual;
    
    GAVG_all_peak.label([1:length(TEP_gavg_ind.label)])=TEP_gavg_ind.label;
    
    GAVG_all_peak.label{length(TEP_gavg_ind.label)+1}='F4FC4';
    GAVG_all_peak.label{length(TEP_gavg_ind.label)+2}='peakP15';
    
    % P15
    channel_idx_p15 = [length(TEP_gavg_ind.label)+1];
    
    time_window = ones(length(subjects),2);
    
    for iI=1:length(TEP_gavg_ind.individual(:,1))
        
        % ADJUST WINDOWS FOR EACH SUBJ & CONDITIONS
        original_window = [0.007 0.025];
        late_window = [0.015 0.020];
        early_window = [0.007 0.018];
        verylate_window = [0.018 0.025];
        
        time_window(iI,:) = original_window;
        
        if contains(cond,'Bi_PA')
            
            if contains(cond, 'Bi_PA_Contr')
                if iI == 4 || iI == 9
                    time_window(iI,:) = late_window;
                elseif iI == 7
                    time_window(iI,:) = early_window;
                end
            else
                
                if iI == 4 || iI == 9
                    time_window(iI,:) = late_window;
                elseif iI == 23
                    time_window(iI,:) = early_window;
                end
                
            end      
            
        elseif cond == 'Bi_AP'
            if iI == 23
                time_window(iI,:) = verylate_window;
            end
            
        elseif cond == 'Bi_LM'
            if iI == 25
                time_window(iI,:) = verylate_window;
            end
        end
        
        
        samples_p15= find(GAVG_all_peak.time>=time_window(iI,1) & GAVG_all_peak.time<=time_window(iI,2));
        
        data_p15 = GAVG_all_peak.individual(iI, channel_idx_p15, samples_p15 ); % pull out the data for this channel and time window
        
        %amplitude
        p_peak_a_p15(iI)=max(data_p15(:,:));
        
        %latency
        p_peak_l_p15(iI) = GAVG_all_peak.time( samples_p15( find( data_p15(:,:)==max(data_p15(:,:)) )  ) );% get the peak latency in this window (max() can be replaced with min() for a negative peak)
        
        %calculate window around peak
        p15_peak_l_min= find(GAVG_all_peak.time==p_peak_l_p15(iI))-25;
        p15_peak_l_max= find(GAVG_all_peak.time==p_peak_l_p15(iI))+25;
        p15_peak_l_window(iI,:)=GAVG_all_peak.time(p15_peak_l_min:p15_peak_l_max);
        p15_peak_a_window(iI,:)=mean(GAVG_all_peak.individual(iI, channel_idx_p15, (p15_peak_l_min:p15_peak_l_max)),3);
        
        EEGm_max1(iI)=max(max(GAVG_all_peak.individual(iI,:,:)));
        GAVG_all_peak.individual(iI, length(TEP_gavg_ind.label)+2, [samples_p15(1) samples_p15(length(samples_p15))])= max(EEGm_max1(iI));
        GAVG_all_peak.individual(iI, length(TEP_gavg_ind.label)+2, samples_p15( find( data_p15(:,:)==max(data_p15(:,:))))-1:samples_p15( find( data_p15(:,:)==max(data_p15(:,:))))+1)= p_peak_a_p15(iI);
        
        
        
    end
    
    %plot (uncomment the following line after the time windows of all subjects have been adjusted
    cfg = [];
    cfg.latency = [-0.02 0.2];
    tmp = ft_selectdata(cfg, GAVG_all_peak);
    tmp.trial = {};
    tmp.trial = GAVG_all_peak.individual;
    tmp.dimord = 'trial_chan_time';
    cfg.highlight = 'labels';
    cfg.viewmode = 'butterfly';
    cfg.layout = 'EEG1010_AZ.lay';
    ft_databrowser(cfg, tmp);
    title(cond)
    
    
    
    outdir = ''; %output directory for saving files
    cd(outdir)
    %%uncomment the following line after the time windows of all subjects have been adjusted
    P15_amplat = cat(2, p_peak_l_p15', p15_peak_a_window);
    save([outdir '\P15_peaks_allsubj_' cond '.txt'], 'P15_amplat', '-ascii')
    
    clearvars -except Bi_AP Bi_LM Bi_PA Bi_PA_Contr Mo_AP Mo_LM Mo_PA cond c allcond subjects
    
end