%% 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);

%% FREQUENCY DOMAIN ANALYSIS
tic
freqact = cell(1, nsubj);
freqbl  = cell(1, nsubj);
for i=1:length(subjects)
    datafile = fullfile(path_data, 'data', subjects{i}, 'erp.mat');
    
    % load the ERP data of subject i from disk
    load(datafile);
    
    % split the data into a pre-stimulus and post-stimulus window
    cfg         = [];
    cfg.latency = [-1 0-1./300]; % subtract 1 sample at the end to get a
                  % nice number of samples, and nice frequency spacing
    bl          = ft_selectdata(cfg, erp);
    cfg.latency = [0.5 1.5-1./300];
    act         = ft_selectdata(cfg, erp);

    % compute power spectrum of the ERP in the specified latency windows
    cfg         = [];
    cfg.foilim  = [0 45]; 
    cfg.channel = {'all' '-E30' '-E3', '-E4', '-E11', '-E12', '-E192'}; 
    cfg.method  = 'mtmfft';
    cfg.taper   = 'hanning';
    freqact{i}  = ft_freqanalysis(cfg, act);
    freqbl{i}   = ft_freqanalysis(cfg, bl);
end

% compute grandaverage spectra
gafreqact = ft_freqgrandaverage([], freqact{:});
gafreqbl  = ft_freqgrandaverage([], freqbl{:});

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

% plot multiple channels representation
cfg            = [];
cfg.layout     = egi256customlay;
cfg.graphcolor = 'rk';
cfg.xlim       = [20 45];

figure;
ft_multiplotER(cfg, gafreqact, gafreqbl);

% plot single channel representation
cfg            = [];
cfg.ylim       = [0 13e-3];
cfg.channel    = {'E6', 'E7', 'E8', 'E9', 'E15', 'E16', 'E17', 'E23', 'E186', 'E198'};
cfg.xlim       = [20 45];
cfg.linewidth  = 2;
cfg.graphcolor = 'rk';


figure;
ft_singleplotER(cfg, gafreqact, gafreqbl);

% 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('frequency [Hz]');
ylabel('power (\muV^2/Hz)');
legend('act','bl');
toc

%% STATISTICAL EVALUATION OF FREQUENCY DOMAIN DATA
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];

cfg                  = [];
cfg.channel          = {'all' '-E30'};% E30 bad
cfg.frequency        = 40;
cfg.method           = 'montecarlo';
cfg.clusterthreshold = 'nonparametric_common';
cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.0001;
cfg.alpha            = 0.025;
cfg.numrandomization = 1000;
cfg.neighbours       = neighbours;
cfg.statistic        = 'ft_statfun_depsamplesT';
cfg.design           = design;
cfg.uvar             = 1;
cfg.ivar             = 2;

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

% plot stat output
cfg           = [];
cfg.layout    = egi256customlay;
cfg.parameter = 'stat';
cfg.highlight = 'on';
cfg.highlightchannel = statfreq.label(statfreq.mask==1);
cfg.colormap  = parula;
cfg.marker    = 'off';
cfg.style     = 'fill';
cfg.comment   = 'off';
cfg.colorbar  = 'yes';
cfg.zlim      = [-6 6];

figure;
subplot(2,1,1);ft_topoplotER(cfg,statfreq);c = colorbar;
c.LineWidth = 1;
c.FontSize = 18;
title(c,'t-val');
toc