%% 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 path + 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'};
nsubj = numel(subjects);
        
%% TIME-FREQUENCY DOMAIN ANALYSIS
tfreq  = cell(1,nsubj);
btfreq = cell(1,nsubj);
for i = 1:nsubj
    datafile = fullfile(path_data, 'data', subjects{i}, 'erp.mat');
    
    % load the ERP data of subject i from disk
    load(datafile);
    
    % compute time-freq representation TFR
    cfg              = [];
    cfg.output       = 'pow';
    cfg.channel      = {'all','-E30'};
    cfg.method       = 'mtmconvol';
    cfg.taper        = 'hanning';
    cfg.foi          = 0:2:46;                          
    cfg.toi          = -.5:0.05:2;               
    cfg.t_ftimwin    = ones(length(cfg.foi),1).*0.5;  
    tfreq{i}         = ft_freqanalysis(cfg, erp);

    cfg              = [];
    cfg.baselinetype = 'relchange';
    cfg.baseline     = [-0.5 -0.25];
    btfreq{i}        = ft_freqbaseline(cfg, tfreq{i});
end

% grand average time freq data
gatfreq = ft_freqgrandaverage([],tfreq{:});

% this loads in a file that contains an electrode layout that can be used
% for plotting
load(fullfile(path_data, 'templates', 'egi256customlay.mat'));

% plot grand average time-freq data
cfg          = [];
cfg.channel  = {'all','-E3', '-E11', '-E12', '-E19', '-E26', '-E33', '-E192'};
cfg.layout   = egi256customlay;
cfg.zlim     = [0 16e-3];
cfg.ylim     = [30 45];
cfg.colormap = parula;
figure;
ft_multiplotTFR(cfg,gatfreq);


% plot single TFR average across electrodes of interest
cfg          = [];
cfg.zlim     = [0 16e-3];
cfg.ylim     = [30 45];
cfg.channel  = {'E6', 'E7', 'E8', 'E9', 'E15', 'E16', 'E17', 'E23', 'E186', 'E198'};
cfg.colormap = parula;

figure;
ft_singleplotTFR(cfg,gatfreq);

% some MATLAB low-level manipulation to improve the figure's aesthetics
set(gca, 'Fontsize',20);
title ('Mean over central electrodes');
set(gca,'box','on');
xlabel('time [sec]');
ylabel('frequency [Hz]');
c = colorbar;
c.LineWidth = 1;
c.FontSize  = 18;
title(c,'\muV^2/Hz');

%% STATISTICAL COMPARISON STIMULATION AGAINST BASELINE

% split data into activity and baseline
freqbl  = cell(1,nsubj);
freqact = cell(1,nsubj);
for i = 1:nsubj
    cfg            = [];
    cfg.latency    = [-.5 0];
    freqbl{i}      = ft_selectdata(cfg,tfreq{i});
    cfg.latency    = [0 .5];
    freqact{i}     = ft_selectdata(cfg,tfreq{i});
    freqbl{i}.time = freqact{i}.time; % FieldTrip only allows for a direct
                     %comparison if the time-axes are equal
end

tic

% prepare electrode neighbours needed for spatial clustering
elec = ft_read_sens(fullfile(path_ft, 'template', 'electrode/GSN-HydroCel-256.sfp'));

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

% prepare the 'design' matrix for the statistical evaluation
design = [1:nsubj 1:nsubj; ones(1,nsubj) ones(1,nsubj)*2];

% test the null-hypothesis of exchangeability between conditions
cfg                  = [];
cfg.channel          = {'all' '-E30'};% E30 broken
cfg.avgoverfreq      = 'no';
cfg.method           = 'montecarlo';
cfg.clusterthreshold = 'nonparametric_common';
cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan        = 2;
cfg.neighbours       = neighbours;
cfg.tail             = 0;
cfg.clustertail      = 0;
cfg.alpha            = 0.025;
cfg.numrandomization = 1000;
cfg.statistic        = 'ft_statfun_depsamplesT';
cfg.design           = design;
cfg.uvar             = 1; % row in design indicating subjects, repeated measure
cfg.ivar             = 2; % row in design indicating condition for contrast

statfreq = ft_freqstatistics(cfg, freqact{:}, freqbl{:});

% plot stat output
cfg           = [];
cfg.layout    = egi256customlay;
cfg.channel   = {'E6', 'E7', 'E8', 'E9', 'E15', 'E16', 'E17', 'E23', 'E186', 'E198'};
cfg.parameter = 'stat';
cfg.colormap  = parula;
cfg.ylim      = [30 45];
cfg.marker    ='off';
cfg.style     = 'fill';
cfg.comment   = 'off';
cfg.maskparameter = 'mask';
cfg.maskstyle = 'outline';
cfg.colorbar  = 'yes';
cfg.zlim      = [-5 5];

figure;
ft_singleplotTFR(cfg,statfreq);

% some MATLAB low-level manipulation to improve the figure's aesthetics
set(gca, 'Fontsize',20);
title ('Mean over central electrodes');
set(gca,'box','on');
xlabel('time [sec]');
ylabel('frequency [Hz]');
c = colorbar;
c.LineWidth = 1;
c.FontSize  = 18;
title(c,'t-val');
toc