% some preliminaries
clear 
close all
clc
% initialize fieldtrip
path_ft   = '/path/to/fieldtrip';
path_data = '/path/to/data/ASSR';

% this adds fieldtrip to the matlab pat h + the relevant subdirectories
% addpath(path_ft);
% ft_defaults;

% the variable subjects contains the list of subject-specific directories
subjects = {'45048'	'47723'	'49449'	'48288'	'45670'...
            '49353'	'47522'	'44340'	'45880'	'45611'...
            '46658'	'48993'	'45859'	'45644'	'47790'...
            '45549'	'48411'	'47371'	'48135'	'46637'...
            '45631'	'47498'	'44296'	'49121'	'46659'...
            '49152'	'44212'	'49004'	'45249'};
% this loads in a file that contains an electrode layout that can be used
% for plotting, but in the next step is used to specify the electrodes that
% will be used for the analysis
load(fullfile(path_data, 'templates', 'egi256customlay.mat'));
%% DATA EPOCHING AND GENERAL PREPROCESSING
tic 

for i=1:length(subjects)
    datapath = strcat(path_data, '/data/', subjects{i});
    dataset  = fullfile(datapath, strcat(subjects{i},'.RAW'));
    
    % in this step the epochs-of-interest are defined
    cfg           = [];
    cfg.dataset   = dataset;
    cfg.trialfun  = 'ft_trialfun_general'; 
    cfg.trialdef.eventtype = 'trigger';
    cfg.trialdef.eventvalue = {'DIN3','DIN5'};
    cfg.trialdef.prestim    = 1;
    cfg.trialdef.poststim   = 3;
    cfg           = ft_definetrial(cfg);
    trl           = cfg.trl; % this field contains the definition of the epochs of interest 
    
    % in this step the data is read from file, and basic preprocessing is
    % applied
    
    cfg.dataset    = dataset;
    cfg.channel    = egi256customlay.label;
    cfg.reref      = 'yes';
    cfg.refchannel = 'all';
    cfg.bpfilter   = 'yes';
    cfg.bpfreq     = [1 48];
    cfg.bpfilttype = 'firws';
    cfg.outputfile = fullfile(datapath, 'data');
    ft_preprocessing(cfg);
end
toc % Elapsed time is 1098.461898 seconds.

%% COMPUTATION OF EVENT-RELATED POTENTIALS AND OF THE GROUP AVERAGE
tic
for i=1:length(subjects)
    datapath = strcat(path_data, '/data/', subjects{i});
    datafile = fullfile(datapath, 'data.mat');
      
    load(datafile);
    
    % average across trials
    cfg         = [];
    cfg.channel = {'all' '-E30' '-E192'};
    erp         = ft_timelockanalysis(cfg,data);
    save(fullfile(datapath, 'erp.mat'), 'erp');
end

% load single subject ERPs in a single cell array
allerp = cell(1,length(subjects));
for i=1:length(subjects)
    datapath     = strcat(path_data, '/data/', subjects{i});
    datafile     = fullfile(datapath, 'erp.mat');
    load(datafile);
    allerp{i} = erp;
end

% compute grand average across subjects
cfg         = [];
cfg.latency = [-.5 2];
gaerp       = ft_timelockgrandaverage(cfg,allerp{:});
toc

%% VISUALIZATION OF THE GROUP AVERAGE ERP

% plot a sensor layout with the ERP time course at each electrode
cfg            = [];
cfg.layout     = egi256customlay;
cfg.graphcolor = 'r';

figure;
ft_multiplotER(cfg,gaerp);

% plot mean ERP over central electrodes
cfg            = [];
cfg.ylim       = [-1 1];
cfg.channel    = {'E6', 'E7', 'E8', 'E9', 'E15', 'E16', 'E17', 'E23', 'E186', 'E198'};
cfg.linewidth  = 2;
cfg.graphcolor = 'r';
figure;
ft_singleplotER(cfg,gaerp);

% below are some 'low-level' matlab commands that do some window dressing
% on the figure
set(gca, 'Fontsize',20);
title('Mean over central electrodes');
set(gca,'box','on');
xlabel('time [sec]');
ylabel('amplitude (\muV)');

% plot mean ERP over central electrodes toghether with STD 
% compute grand average over subjects but keep individual representations
% in the output
cfg         = [];
cfg.latency = [-.5 2];
cfg.keepindividual = 'yes';
gaerpind        = ft_timelockgrandaverage(cfg,allerp{:});
% and average over channels
cfg = [];
cfg.channel     = {'E6', 'E7', 'E8', 'E9', 'E15', 'E16', 'E17', 'E23', 'E186', 'E198'};
cfg.avgoverchan = 'yes';
gaerpind        = ft_selectdata(cfg,gaerpind);
x = gaerpind.time'; % x-axis def
y = mean(squeeze(gaerpind.individual),1)'; % y-axis def
e = std(squeeze(gaerpind.individual),1)';  % std
low = y - e;% lower bound
high = y + e;% upper bound
figure;
hp = patch([x; x(end:-1:1); x(1)], [low; high(end:-1:1); low(1)], 'r');
hold on;
hl = line(x,y);
set(hp, 'facecolor', [1 0.8 0.8], 'edgecolor', 'none');
set(hl, 'color', 'r','linewidth',2);
set(gca,'FontSize',20);
title ('Mean and STD over central electrodes');
set(gca,'box','on');
xlabel('time [sec]');
ylabel('amplitude (\muV)');

% plot topography of the N1 and P2 components
cfg          = [];
cfg.layout   = egi256customlay;
cfg.xlim     = [0.17 0.23];
cfg.zlim     = [-1 1];
cfg.colormap = parula;
cfg.marker   = 'off';
cfg.style    = 'fill';
cfg.comment  = 'off';
cfg.colorbar = 'yes';

figure;
subplot(2,1,1);ft_topoplotER(cfg,gaerp);
title('N1 topography')
c = colorbar;
c.LineWidth = 1;
c.FontSize = 18;
title(c,'\muV')

cfg.xlim = [0.25 0.3];
subplot(2,1,2);ft_topoplotER(cfg,gaerp);
title('P2 topography');
c = colorbar;
c.LineWidth = 1;
c.FontSize = 18;
title(c,'\muV');
%% STATISTICAL EVALUATION OF TIME DOMAIN DATA

% create a 'neighbourhood-structure', that defines how spatial clusters are
% formed
tic
[~,ftpath] = ft_version;
elec       = ft_read_sens(strcat(ftpath,'/template/electrode/GSN-HydroCel-256.sfp'));

cfg               = [];
cfg.method        = 'distance';
cfg.neighbourdist = 3;
cfg.feedback      = 'yes';
neighbours        = ft_prepare_neighbours(cfg,elec);

% split data into activity and baseline
act = cell(1,length(subjects));
bl  = cell(1,length(subjects));
for i=1:length(subjects)
    cfg = [];
    cfg.latency = [-.5 0];
    bl{i} = ft_selectdata(cfg,allerp{i});
    cfg.latency = [0 .5];
    act{i} = ft_selectdata(cfg,allerp{i});
    bl{i}.time = act{i}.time;
end

% create a FieldTrip style 'design'-matrix, conditions: 1 = act, 2 = bl

subj   = length(subjects);
design = zeros(2,2*subj);
for i = 1:subj
  design(1,i) = i;
end
for i = 1:subj
  design(1,subj+i) = i;
end
design(2,1:subj)        = 1;
design(2,subj+1:2*subj) = 2;

% compute the difference as a T-statistic, and do an inferential test
cfg                  = [];
cfg.channel          = {'all' '-E30'};% E30 broken
cfg.method           = 'montecarlo';
cfg.statistic        = 'ft_statfun_depsamplesT';
cfg.clusterthreshold = 'nonparametric_common';
cfg.correctm         = 'cluster';
cfg.alpha            = 0.025;% .05/2 because two-sided
cfg.numrandomization = 1000;
cfg.neighbours       = neighbours;
cfg.design           = design;
cfg.uvar             = 1;
cfg.ivar             = 2;

staterp = ft_timelockstatistics(cfg, act{:},bl{:});

% plot stat output
cfg               = [];
cfg.layout        = egi256customlay;
cfg.parameter     = 'stat';
cfg.maskparameter = 'mask';
cfg.graphcolor    = 'r';

figure;
ft_multiplotER(cfg,staterp);

% plot stat difference ERP over central electrodes
cfg            = [];
cfg.ylim       = [-10 6];
cfg.channel    = {'E6', 'E7', 'E8', 'E9', 'E15', 'E16', 'E17', 'E23', 'E186', 'E198'};
cfg.parameter  = 'stat';
cfg.maskparameter = 'mask';
cfg.linewidth  = 2;
cfg.graphcolor = 'r';

figure;
ft_singleplotER(cfg,staterp);
set(gca, 'Fontsize',20);
title ('Mean over central electrodes');
set(gca,'box','on');
xlabel('time [sec]');
ylabel('act vs. bl (t-val)');

% plot N1 and P2 topography in units of t-value
cfg           = [];
cfg.layout    = egi256customlay;
cfg.xlim      = [0.17 0.23];
cfg.parameter = 'stat';
cfg.zlim      = [-5 5];
cfg.colormap  = parula;
cfg.marker    = 'off';
cfg.style     = 'fill';
cfg.comment   = 'off';
cfg.colorbar  = 'yes';
figure;
subplot(2,1,1);ft_topoplotER(cfg,staterp);
title('N1 topography')
c = colorbar;
c.LineWidth = 1;
c.FontSize = 18;
title(c,'t-val')
cfg.xlim = [0.25 0.3];
subplot(2,1,2);ft_topoplotER(cfg,staterp);
title('P2 topography');
c = colorbar;
c.LineWidth = 1;
c.FontSize = 18;
title(c,'t-val');
toc
%% plot analysis pipeline in browser
cfg = [];
cfg.filename = 'staterp';
cfg.filetype = 'html';
ft_analysispipeline(cfg, staterp);

!open staterp.html