%%%%%%%%%%%%%%
%This script starts from raw BP files
%This version is set for time-frequency -- no low pass filtering
%%%%%%%%%%%%%%

%% create a subject list
subjectList={'1',
'2',
'3',
'4',
'5',
'6',
'7',
'8',
'9',
'10',
'11',
'12',
'13',
'14',
'15',
'16',
'17',
'18',
'19',
'20',
'21',
'22',
'23',
'24',
'25',
'26',
'27',
'28',
'29',
'30',
'31',
'32',
'33',
'34',
'35',
'36'}


thresholdList=[100,
100,
110,
80,
120,
100,
75,
150,
120,
75,
120,
100,
75,
74,
80,
100,
100
100,
80,
100,
100,
100,
100,
140,
100,
100,
100,
100,
100,
100,
60,
100,
100,
75,
80,
80]

%% run a loop
openPath                            = '/Users/jsprouse/Desktop/class data/'

for i = 14:length(subjectList)
	%clear whatever happened in the last loop, keep the subjectList and the variable i

	clearvars -except subjectList i openPath thresholdList
	
	%choose a subject from the list
	subject                         = char(subjectList(i));
    threshold = thresholdList(i);
    
%%%%%%%%%%%
%Use ft_definetrial to read in unaccusative and unergative trials
%%%%%%%%%%%

	%specify the dataset
	cfg 							= [];
	cfg.dataset             		= [openPath subject '.vhdr'];
	cfg.headerfile 					= [openPath subject '.vhdr'];

	%define the trials of interest
	cfg.trialdef.eventtype 			= 'Stimulus';
	cfg.trialdef.prestim 			= 1;
	cfg.trialdef.poststim 			= 1;
	cfg.trialdef.eventvalue         = {'55' '66' '77' '88' '155' '166' '177' '188'};

	%choose a trial function (I modified this one to not look for "S" and "R")
	cfg.trialfun					='ft_trialfun_general_modified';
	
	%run the definetrial function
	cfg_N400          				= ft_definetrial(cfg);

%%%%%%%%%%%%%%%
%Use ft_preprocessing to do some preprocessing
%%%%%%%%%%%%%%%

	%reference
	cfg_N400.channel                = 'all';	
	cfg_N400.reref                  = 'yes';
	cfg_N400.refchannel             = {'M1' 'M2'};
	
	%high-pass filter (no need for low-pass because this is a time-frequency pipeline)
	cfg_N400.hpfilter               = 'yes';
	cfg_N400.hpfreq                 = .1;
	%cfg_N400.hpfiltord				= 2;
	cfg_N400.hpfilttype             = 'firws';
	%cfg_N400.lpfilter              	= 'yes';
	%cfg_N400.lpfreq              	= 30;
	cfg_N400.hpinstabilityfix 		= 'reduce';		%This tends to fail without reducing the order (6 is default, 3 worked here)
	
	%remove line noise (just in case we look this high)?
	cfg_N400.bsfilter               = 'yes';
	cfg_N400.bsfreq                 = [59 61]	%60hz
	
	%demean and detrend? (detrending really helps with drifts; not sure what, if any, effect it would have on sustained ERPs)
	cfg_N400.demean                 ='yes';
	cfg_N400.detrend                ='yes';

	%specify the channel locations
	cfg_N400.layout                 = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
	
	%run the preprocessing function
	N400                            = ft_preprocessing(cfg_N400)

    %do we need to interpolate any channels?
    if i==14
        %define neighbors manually
        neighbors                   = [];
        neighbors.label             = {'P7'};
        neighbors.neighblabel       = {'O1' 'P3' 'CP5'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'P7'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
        
        neighbors                   = [];
        neighbors.label             = {'T7'};
        neighbors.neighblabel       = {'FT9' 'F7' 'C3' 'FC5' 'CP5' 'P7'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'T7'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
    end
    
    if i==16
        %define neighbors manually
        neighbors                   = [];
        neighbors.label             = {'T8'};
        neighbors.neighblabel       = {'FT10' 'F8' 'C4' 'FC6' 'CP6' 'P8'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'T8'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
        
        neighbors                   = [];
        neighbors.label             = {'T7'};
        neighbors.neighblabel       = {'FT9' 'F7' 'C3' 'FC5' 'CP5' 'P7'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'T7'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
    end
    
    if i==19
        %define neighbors manually
        neighbors                   = [];
        neighbors.label             = {'F7'};
        neighbors.neighblabel       = {'F3' 'FC5'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'F7'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
    end
    
    if i==22
        %define neighbors manually
        neighbors                   = [];
        neighbors.label             = {'O2'};
        neighbors.neighblabel       = {'Pz' 'P4' 'P8'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'O2'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
        
        neighbors                   = [];
        neighbors.label             = {'O1'};
        neighbors.neighblabel       = {'Pz' 'P3' 'P7'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'O1'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
        
        neighbors                   = [];
        neighbors.label             = {'Oz'};
        neighbors.neighblabel       = {'Pz' 'P3' 'P4' 'O1' 'O2'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'Oz'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
    end
    
    if i==30
        %define neighbors manually
        neighbors                   = [];
        neighbors.label             = {'T7'};
        neighbors.neighblabel       = {'F7' 'P7' 'CP5' 'FC5' 'C3'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'T7'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
        
        neighbors                   = [];
        neighbors.label             = {'CP5'};
        neighbors.neighblabel       = {'P3' 'P7' 'C3' 'CP1' 'T7'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'CP5'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
    end
    
    if i==33
        %define neighbors manually
        neighbors                   = [];
        neighbors.label             = {'T7'};
        neighbors.neighblabel       = {'F7' 'P7' 'CP5' 'FC5' 'C3'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'T7'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
        
        neighbors                   = [];
        neighbors.label             = {'T8'};
        neighbors.neighblabel       = {'F8' 'P8' 'CP6' 'FC6' 'C4'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'T8'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
        
        neighbors                   = [];
        neighbors.label             = {'Fp2'};
        neighbors.neighblabel       = {'Fp1' 'Fz' 'F4' 'F8'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'Fp2'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
        
        neighbors                   = [];
        neighbors.label             = {'FT10'};
        neighbors.neighblabel       = {'F8' 'FC6' 'T8'}
        
        cfg_repair                  = [];
        cfg_repair.layout           = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20181210/template/layout/acticap-64ch-standard2.mat';
        cfg_repair.badchannel       = {'FT10'};
        cfg_repair.neighbours       = neighbors
        N400                        = ft_channelrepair(cfg_repair, N400)
    end
    
%%%%%%%%%%%%%%
%Artifact detection
%%%%%%%%%%%%%%

%Voltage Threshold detection

%it seems worthwhile to filter for threshold-based detection, since EMG might trip the threshold
%Since noisy channels can trip this, you may need to specify different thresholds for different participants; or only look at certain channels.
%Filtering is strange... it seems to still mark based on the raw values. Not sure why.
	
	cfg_art                                 = [];
	cfg_art.trl                             = [N400.sampleinfo N400.trialinfo];
	cfg_art.continuous                      = 'no';
	cfg_art.artfctdef.threshold.bpfilter    = 'no';
	cfg_art.artfctdef.threshold.bpfiltord   = 3;
	cfg_art.artfctdef.threshold.bpfreq      = [.1 30];
	cfg_art.artfctdef.threshold.min         = -threshold;
	cfg_art.artfctdef.threshold.max         = threshold;
	cfg_art.artfctdef.threshold.range       = 2*threshold;
	cfg_art.artfctdef.threshold.channel     = {'all'};

	%run detection
	[cfg_art]                               = ft_artifact_threshold(cfg_art, N400);

	%data browser to see the trials
	%cfg_art.ylim                           = [-200 200];
	%cfg                                    = ft_databrowser(cfg_art, N400);

	%run rejection and save
	N400_clean                              = ft_rejectartifact(cfg_art, N400);


%%%%%%%%%%%%%%%%%%%
%save for future processing
%%%%%%%%%%%%%%%%%%%

	%save the results
	saveFileName                            = [subject '.N400_epoch.mat'];
	savePath                                = '/Users/jsprouse/Desktop/preprocessed/';
	saveFile                                = [savePath saveFileName];
	
	save(saveFile,'N400_clean');

end	%this ends the loop

%voltage thresholds by subject
%1 100
%2 100 - a couple of jump-like artifacts weren't caught by the threshold...  
%3 110
%4 80
%5 120
%6 100
%7 75
%8 150 - TOO MANY ARTIFACTS
%9 120 - borderline artifacts (~20 trials left per condition)
%10 75
%11 120
%12 100 - TOO MANY ARTIFACTS
%13 75
%14 74 - interpolated T7 and P7
%15 80 - borderline artifacts (~20 trials left per condition)
%16 100 - interpolated T7 and T8
%17 100
%18 100 - borderline artifacts (~20 trials left per condition)
%19 80 - interpolated F7
%20 100
%21 100
%22 100 - interpolated O1, O2, Oz
%23 100
%24 140
%25 100 - TOO MANY ARTIFACTS
%26 100
%27 100
%28 100
%29 100
%30 100 - interpolated T7 and CP5
%31 60
%32 100 - borderline too many artifacts (~20 trials left in each condition)
%33 100 - interpolated T7, T8, Fp2, FT10
%34 75
%35 80
%36 80

%% Figure out how many trials are remaining for each participant?

openPath                                    = '/Users/jsprouse/Desktop/1. N400 epochs/';
trialCount                                  = zeros(1,length(subjectList))
for i = 1:length(subjectList)
	%clear whatever happened in the last loop, keep the subjectList and the variable i

	clearvars -except subjectList i openPath thresholdList trialCount
	
	%choose a subject from the list
	subject                                 = char(subjectList(i));
    
    openFileName                            = [subject '.N400_epoch.mat'];
	openFile                                = [openPath openFileName];

    load(openFile);

    trialCount(i)                           = size(N400_clean.trialinfo,1);
end