%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%   SSVEP SYL MAIN   %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Veronica Montani

% this script accompanies the manuscript titled 
%    "Frequency-tagged visual evoked responses track syllable effects in visual word recognition"
%   from Veronica Montani, Valerie Chanoine, Jonathan Garinger, Johannes
%   Ziegler.

% You will need the following additional files: 
%              proc_data_syl_100mkv_Int_thr25_filt5_100_trig.mat
%              Fieldtrip toolbox
% question -> veronica.montani@live.com

clear all;close all; %clc;

ft_defaults

%%% Define work folder
% cd ../../../


%% load data
my_data = load('proc_data_syl_100mkv_Int_thr25_filt5_100_trig.mat');


EXP.subj=1:17;
nSubj=numel(EXP.subj);
sR = 512; % sampling rate, (samples/unit time)
cond = {'Cong','Inc'}; % conditions: congruent vs incongruent

% convert eeglab format + addiging events
data = [];

for isbj = EXP.subj
    data(isbj).EGB = eeglab2fieldtrip(my_data.EXP.EEG{1,isbj},'preprocessing');
    data(isbj).EGB.event = {my_data.EXP.EEG{1,isbj}.event(1:end).type}; % events (epochs labels)
end



%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% STATIC SSVEP 18.75 Hz %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1) frequency decomposition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% redefine trials

cyl_sec = 1/18.75;                                % 1 second diveded for the frequency of interest = 1 cycle
start_18 = ((floor(0.55/cyl_sec)) * cyl_sec);     % 0.55 = 533
end_18 = 1.76; % hard-coded, final freq 18.73 


% cut into start_18 to end_18 window
cfg = [];
cfg.toilim = [start_18 end_18]; % time window in seconds 
data(nSubj).tr =[];
for isbj = EXP.subj
    data(isbj).tr = ft_redefinetrial(cfg, data(isbj).EGB);
end

%% detrend
cfg = [];

cfg.detrend = 'yes';
for isbj = EXP.subj
    data(isbj).tr = ft_preprocessing(cfg, data(isbj).tr);
end


%% freq decomposition
% 
cfg                 = [];
cfg.method          = 'mtmfft';
cfg.output          = 'powandcsd';
cfg.taper           = 'hanning';
cfg.channel         = (1:64)';

cfg.tapsmofrq = 4; % 4
cfg.foi    = 18.75;

fft_syl =[];
for isbj = EXP.subj
    cTr = 1:numel(data(isbj).tr.trial);        % vectors all trials
    for iType=1:numel(cond)
        cTy_ind = strcmp(cond{iType},data(isbj).EGB.event);
        cfg.trials = cTr(cTy_ind);             % selection of trials accordingly with cond
        fft_syl.(cond{iType}){isbj} = ft_freqanalysis(cfg, data(isbj).tr);
    end
end


% save('fft_syl.mat', '-struct', 'fft_syl'); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2) cluster based permutation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%  FT_STATFUN_DIFF_ 
addpath /Applications/MATLAB_R2014a/FIELDTRIP/fieldtrip-20160211/template/neighbours

cfg                 = [];   
cfg.channel = 'all';

% FREQ
sel_frq = 18.75;

% CONFIG           
cfg.method    = 'montecarlo'; 
cfg.statistic = 'ft_statfun_depsamplesT';
cfg.tail = 0;
cfg.alpha = 0.025; 

% cluster
cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.05;
cfg.clusterstatistic = 'wcm';
cfg.wcm_weight = 1.5;
cfg.minnbchan        = 2;
cfg.clustertail      = 0;
% prepare_neighbours determines what sensors may form clusters
cfg_neighb.method    = 'triangulation';
cfg_neighb.channel = 'all';
cfg.neighbours       = ft_prepare_neighbours(cfg_neighb, fft_syl.Cong{1,1});

cfg.numrandomization = 2000;
cfg.frequency = sel_frq;

% design
Nsubj = length(fft_syl.Cong);
nCond = 2;
design1 = []; for n=1:nCond, design1 = [design1 ones(1,Nsubj)*n]; end
design2 = repmat(1:Nsubj, [1,nCond]);
cfg.design(1,:) = design1;
cfg.design(2,:) = design2;
cfg.ivar  = 1;
cfg.uvar  = 2;


cfg.parameter = 'powspctrm';


% FREQSTAT                                                                 
stat_F_syl= ft_freqstatistics(cfg,fft_syl.Cong{:},fft_syl.Inc{:}); 


%% topoplot stat
% find cluster labels
cls_ind = find(stat_F_syl.mask);
cls_lab = stat_F_syl.label(cls_ind);
% cfg
cfg = [];
cfg.alpha  = 0.025;
cfg.parameter = 'stat';
cfg.highlightsymbolseries     = ['*','*','+','o','.'] %1x5 vector, highlight marker symbol series (default ['*','x','+','o','.'] for p < [0.01 0.05 0.1 0.2 0.3]
cfg.highlightsizeseries = [6 14 6 6 6]; % 1x5 vector, highlight marker size series   (default [6 6 6 6 6] for p < [0.01 0.05 0.1 0.2 0.3])
cfg.layout = 'biosemi64.lay';
cfg.subplotsize = [1 1];% layout of subplots ([h w], default [3 5])
cfg.colorbar = 'yes';
cfg.gridscale = 1000;
% plot
ft_clusterplot(cfg, stat_F_syl)




%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% DYNAMIC SSVEP 18.75 Hz %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1) frequency decomposition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% redefinetrial for freq-time
cyl_sec = 1/18.75;                                % 1 second diveded for the frequency of interest = 1 cycle
start_18 = -0.1992;                                    % time x freq
end_18 = 1.9341;          % hard coded 

% cut into start_18 end_18 window
cfg = [];
cfg.toilim = [start_18 end_18]; % time window in seconds 
data(nSubj).tr =[];
for isbj = EXP.subj
    data(isbj).tr = ft_redefinetrial(cfg, data(isbj).EGB);
end

%% detrend
cfg = [];

cfg.detrend = 'yes';
for isbj = EXP.subj
    data(isbj).tr = ft_preprocessing(cfg, data(isbj).tr);
end



%% time frequency decomposition

cfg                 = [];
cfg.method          = 'mtmconvol';
cfg.output          = 'powandcsd';
%cfg.output          = 'fourier';
cfg.taper           = 'hanning';
cfg.channel         = (1:64)';

cfg.tapsmofrq = 4; % 4
cfg.foi    = [18.75 25];

% set the time window 
cyWin = cyl_sec * 9; % hard coded, 9 cycles= 0.480 msec
cfg.t_ftimwin       = cyWin * ones(length(cfg.foi), 1);

% set the timepoints of interest:
step = 0.02; 
cfg.toi             = start_18:step:end_18;


fft_tsyl =[];
for isbj = EXP.subj
    cTr = 1:numel(data(isbj).tr.trial);        % vectors all trials
    for iType=1:numel(cond) % condition
        cTy_ind = strcmp(cond{iType},data(isbj).EGB.event);
        cfg.trials = cTr(cTy_ind);             % selection of trials accordingly with cond
        fft_tsyl.(cond{iType}){isbj} = ft_freqanalysis(cfg, data(isbj).tr);
    end
end


save('fft_tsyl.mat', '-struct', 'fft_tsyl'); % this is also for the source


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2) cluster based permutation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%  FT_STATFUN_DIFF (time frequency)
addpath /Applications/MATLAB_R2014a/FIELDTRIP/fieldtrip-20160211/template/neighbours

cfg                 = [];   
cfg.channel = 'all';

% FREQ
sel_frq = 18.75;

% CONFIG           
cfg.method    = 'montecarlo'; 
cfg.statistic = 'ft_statfun_depsamplesT';

cfg.tail = 0;
cfg.alpha = 0.025;

cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.05;
cfg.clusterstatistic = 'wcm';
cfg.wcm_weight = 1.5;
cfg.minnbchan        = 2;
cfg.clustertail      = 0;
% prepare_neighbours determines what sensors may form clusters
cfg_neighb.method    = 'triangulation';
cfg_neighb.channel = 'all';
cfg.neighbours       = ft_prepare_neighbours(cfg_neighb, fft_tsyl.Cong{1,1});

cfg.numrandomization = 2000;
cfg.frequency = sel_frq;
cfg.latency = [0.05 1.2];

% design
Nsubj = length(fft_tsyl.Cong);

nCond = 2;

design1 = []; for n=1:nCond, design1 = [design1 ones(1,Nsubj)*n]; end
design2 = repmat(1:Nsubj, [1,nCond]);
cfg.design(1,:) = design1;
cfg.design(2,:) = design2;
cfg.ivar  = 1;
cfg.uvar  = 2;

cfg.parameter = 'powspctrm';


% FREQSTAT                                                                 
stat_TF_syl= ft_freqstatistics(cfg,fft_tsyl.Cong{:},fft_tsyl.Inc{:}); % Cong vs Inc

%save('stat_TF_syl.mat', '-struct', 'stat_TF_syl');


%% plot stat
%%%% options,cfg.comment and cfg.commentpos
cfg = [];
cfg.alpha  = 0.025;
cfg.parameter = 'stat';
cfg.highlightsymbolseries     = ['*','*','+','o','.'] %1x5 vector, highlight marker symbol series (default ['*','x','+','o','.'] for p < [0.01 0.05 0.1 0.2 0.3]
cfg.layout = 'biosemi64.lay';
cfg.gridscale = 1000;
cfg.subplotsize = [2 4];% layout of subplots ([h w], default [3 5])
ft_clusterplot(cfg, stat_TF_syl);

% just to have the colorbar...
% cfg.subplotsize = [1 1];
% ft_clusterplot(cfg, stat_TF_syl); colorbar;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3) common filter for the source
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% time freq decomposition two conditions (to compute a common filter in source...)

cfg                 = [];
cfg.method          = 'mtmconvol';
cfg.output          = 'powandcsd';
%cfg.output          = 'fourier';
cfg.taper           = 'hanning';
cfg.channel         = (1:64)';
%cfg.keeptrials = 'yes';

cfg.tapsmofrq = 4; % 4
cfg.foi    = [18.75 25];

% set the time window 
cyWin = cyl_sec * 9;
cfg.t_ftimwin       = cyWin * ones(length(cfg.foi), 1);

% set the timepoints of interest:
step = 0.02;
cfg.toi             = start_18:step:end_18;


fft_tsyl_cmb =[];

sel_type = [1 2];   % selection of the condition
for isbj = EXP.subj
    cTr = 1:numel(data(isbj).tr.trial);        % vectors all trials
    cTy_ind = ismember(data(isbj).EGB.event,cond(sel_type));
    cfg.trials = cTr(cTy_ind);             % selection of trials accordingly with cond
    fft_tsyl_cmb.CI{isbj} = ft_freqanalysis(cfg, data(isbj).tr);
end

 
save('fft_tsyl_cmb.mat', '-struct', 'fft_tsyl_cmb');


%%%%%%%%%%%%%%%
%%%%% END %%%%%
%%%%%%%%%%%%%%%