% If required, for the string variable root, you must supply the directory 
% in which you can find the data that are used in the tutorial.

root='/home/common/courses/toolkit2005/EEGMEG2005/twente/';
addpath(root);
cd(root);

% 1. Randomization tests for between-trial experiments.

% 1.1 Compare event-related averages

% 1.1.1 Preprocessing the data of subject 01

cfg=[];
cfg.trialfile='Subject01.ds';
cfg.trialfun='maketrlSemCon';
cfg.trialdef.eventtype='backpanel trigger';
cfg.trialdef.eventvalue=3;
cfg.channel={'MEG'};
cfg.dataset='Subject01.ds';
cfg.blc='yes';
cfg.blcwindow=[-0.2 0];
cfg.lpfilter='yes';
cfg.lpfreq=35;
cfg.padding=0;

% Store the output of preprocessing.m in FICpreproc:

[FICpreproc] = preprocessing(cfg);

% We want to compare the fully congruent and fully incongruent condition. Therefore, we also have to extract 
% the trials in the fully congruent condition. These trials were associated the trigger number 9. 
% Thus, in our new configuration we must have cfg.trialdef.eventvalue=9. All other fields of cfg remain the same. 
% Store the output of preprocessing.m in FCpreproc:

cfg.trialdef.eventvalue=9;
[FCpreproc] = preprocessing(cfg);

% Save the preprocessed data in the file preprocdata:

save preprocdata FICpreproc FCpreproc;

% If you want to skip preprocessing.m, load the output from the file preprocdataorig:

load preprocdataorig;

% 1.1.2 Run timelockanalysis

cfg = [];
cfg.keeptrials = 'yes';
FICdata = timelockanalysis(cfg, FICpreproc);
FCdata = timelockanalysis(cfg, FCpreproc);
    
% Save the output of timelockanalysis.m in the file timelockdata:

save timelockdata FICdata FCdata;

% If you want to skip timelockanalysis.m, load the output from the file timelockdataorig:

load timelockdataorig;

% 1.1.3 Perform the randomization test

cfg.statistic  = 'indepsamplesT';
cfg.alphathresh = 0.05; 
cfg.makeclusters = 'yes'; 
cfg.minnbchan = 2;
cfg.clusterteststat = 'maxsum';
cfg.onetwo = 'twosided';
cfg.alpha = 0.05;
cfg.nranddraws = 100;
cfg.channel = {'MEG'};
cfg.latency = [0 1];

[clusrand] = clusterrandanalysis(cfg, FCdata, FICdata)

% Save the output of clusterrandanalysis.m in clusrandFCvsFIC:

save clusrandFCvsFIC clusrand;

% If you want to skip clusterrandanalysis.m, load the output from the file clusrandFCvsFICorig:

load clusrandFCvsFICorig;

% 1.1.4 Plot the results

% construct a mask

mask=(clusrand.posclusterslabelmat==1) | (clusrand.negclusterslabelmat==1);

% All the information topoplotER.m needs will be stored in plotdata:

plotdata.maskedraweffect=clusrand.raweffect.*mask;
plotdata.label=clusrand.label;
plotdata.time=clusrand.time;

% To plot a sequence of ten topographic plots, equally spaced between 0 and 1 second:
 
cfg=[];
cfg.yparam='maskedraweffect';
cfg.xlim=[0.1:0.1:1.0];
cfg.colorbar='no';
clf;
topoplotER(cfg,plotdata);


% 1.2. Compare Time-frequency representations

% 1.2.1 TFR-specific preprocessing (without low-pass filtering)

cfg=[];
cfg.trialfile='Subject01.ds';
cfg.trialfun='maketrlSemCon';
cfg.trialdef.eventtype='backpanel trigger';
cfg.channel={'MEG'};
cfg.dataset='Subject01.ds';
cfg.blc='yes';
cfg.blcwindow=[-0.2 0];

% For the FIC condition, store the output of preprocessing.m in FICpreproc:

cfg.trialdef.eventvalue=3;
[FICpreproc] = preprocessing(cfg);

% For the FC condition, store the output in FCpreproc:

cfg.trialdef.eventvalue=9;
[FCpreproc] = preprocessing(cfg);

% Save the preprocessed data in the file preprocTFRdata:

save preprocTFRdata FICpreproc FCpreproc;

% If you want to skip preprocessing.m, load the output from the file preprocTFRdataorig:

load preprocTFRdataorig;


% 1.2.3. Calculate the TFRs for the two experimental conditions (FC and FIC)

% Use freqanalysis.m with the following configuration:

cfg = [];
cfg.output       = 'pow';
cfg.sgn          = 'MEG';
cfg.method       = 'mtmconvol'; 
cfg.foi          = [5 10 20 40 80];
cfg.t_ftimwin    = 5./cfg.foi;
cfg.tapsmofrq    = 0.5*cfg.foi;						            
cfg.toi          = -0.5:0.05:1.4; 
cfg.pad          = 'maxperlen';
cfg.keeptrials   = 'yes';	           
cfg.keeptapers   = 'no';

% Calculate the TFRs for the two experimental conditions:
 
TFRFC = freqanalysis(cfg, FCpreproc);
TFRFIC = freqanalysis(cfg, FICpreproc);

% To save the TFRs in the file TFR:

save TFR TFRFC TFRFIC;

% If you want to skip freqanalysis.m, load the output from the file TFRorig:

load TFRorig;

% Run clusterrandanalysis.m to compare TFRFC and TFRFIC. 
% Note that, except for the field cfg.latency, the following configuration is identical to the configuration that was used for 
% comparing event-related averages:

cfg=[];
cfg.statistic  = 'indepsamplesT';
cfg.alphathresh = 0.05; 
cfg.makeclusters = 'yes'; 
cfg.minnbchan = 2;
cfg.clusterteststat = 'maxsum';
cfg.onetwo = 'twosided';
cfg.alpha = 0.05;
cfg.nranddraws = 100;
cfg.channel = {'MEG'};
cfg.latency = [0 1.4];

[clusrand] = clusterrandanalysis(cfg, TFRFC, TFRFIC);

% Save the output of clusterrandanalysis.m in clusrandTFRFCvsFIC:

save clusrandTFRFCvsFIC clusrand;

% If you want to skip clusterrandanalysis.m, load the output from the file clusrandTFRFCvsFICorig:

load clusrandTFRFCvsFICorig;

% To show the topography of the first positive cluster, we use topoplotTFR.m:

mask=(clusrand.posclusterslabelmat==1);
plotdata.powspctrm=clusrand.raweffect.*mask;
plotdata.label=clusrand.label;
plotdata.toi=clusrand.time;
plotdata.foi=clusrand.freq;

% To plot a sequence of seven topographic plots for the lowest frequency (5 Hertz), equally spaced between 0 and 1.4 seconds, 
% use the following configuration:
 
cfg=[];
cfg.xlim=[0:0.2:1.4];
cfg.ylim=[5 5];    % Frequency interval for which to plot
cfg.colorbar='no';
clf;
topoplotTFR(cfg,plotdata);

% Also try the other frequencies

% 2. Randomization test for a within-trial experiment.

% 2.1. Calculate the TFRs for the activation and the baseline period

cfg = [];
cfg.output       = 'pow';
cfg.sgn          = 'MEG';
cfg.method       = 'mtmconvol'; 
cfg.foi          = [5 10 20 40 80];
cfg.t_ftimwin    = 5./cfg.foi;
cfg.tapsmofrq    = 0.5*cfg.foi;	            
cfg.pad          = 'maxperlen';
cfg.keeptrials   = 'yes';	           
cfg.keeptapers   = 'no';

%To calculate the TFR for the baseline condition (the epoch preceeding the FIC sentence ending):
 
cfg.toi          = -0.5:0.05:-0.05; 
TFRFICbaseline = freqanalysis(cfg, FICpreproc);

% And to calculate the TFR for the activation condition (the epoch following the FIC sentence ending):
 
cfg.toi          = 0.95:0.05:1.4; 
TFRFICactivation = freqanalysis(cfg, FICpreproc);

% To save the TFRs in the file TFRbasact:

save TFRbasact TFRFICbaseline TFRFICactivation;

% If you want to skip freqanalysis.m, load the output from the file TFRbasactorig:

load TFRbasactorig;


% 2.2. Perform the randomization test

cfg=[];
cfg.statistic  = 'actvsblT';
cfg.alphathresh = 0.05; 
cfg.makeclusters = 'yes'; 
cfg.minnbchan = 2;
cfg.clusterteststat = 'maxsum';
cfg.onetwo = 'twosided';
cfg.alpha = 0.05;
cfg.nranddraws = 100;
cfg.channel = {'MEG'};
cfg.latency = 'all';

[clusrand] = clusterrandanalysis(cfg, TFRFICactivation, TFRFICbaseline);

% Save the output of clusterrandanalysis.m in clusrandTFRFICactvsbas:

save clusrandTFRFICactvsbas clusrand;

% If you want to skip clusterrandanalysis.m, load the output from the file clusrandTFRFICactvsbasorig:

load clusrandTFRFICactvsbasorig;

% 2.3. Plot the results

mask=(clusrand.posclusterslabelmat==1);

% All the information needed by topoplotTFR.m will be stored in plotdata:

plotdata.powspctrm=clusrand.raweffect.*mask;
plotdata.label=clusrand.label;
plotdata.toi=clusrand.time;
plotdata.foi=clusrand.freq;

% To plot a sequence of four topographic plots for the lowest frequency (5 Hertz), 
% equally spaced between 0.95 and 1.35 seconds:
 
cfg=[];
cfg.xlim=[0.95:0.1:1.35];
cfg.ylim=[5 5];    % Frequency interval for which to plot
cfg.colorbar='no';
clf;
topoplotTFR(cfg,plotdata);

% Also try the other frequencies

% 3. Randomization test for a within-subjects experiment

% 3.1. Load the data

% For an explanation about how the data were preprocessed (i.c., averaging,
% planar gradient computation, combined planar computation) see the
% tutorial.

load gravgerfcporig;

% 3.2. Perform the statistical test

cfg=[];
cfg.statistic  = 'depsamplesT';
cfg.alphathresh = 0.05; 
cfg.makeclusters = 'yes'; 
cfg.minnbchan = 2;
cfg.clusterteststat = 'maxsum';
cfg.onetwo = 'twosided';
cfg.alpha = 0.05;
cfg.nranddraws = 100;
cfg.channel = {'MEG'};
cfg.latency = [0 1];

[clusrand] = clusterrandanalysis(cfg, gravg_erf_cp_FIC, gravg_erf_cp_FC);

% If you want to skip clusterrandanalysis.m, load the output from the file clusranderfcpFICvsFCorig:

load clusranderfcpFICvsFCorig;

% 3.3. Plot the results

mask=clusrand.posclusterslabelmat==1;

plotdata.label=clusrand.label;
plotdata.time=clusrand.time;
plotdata.maskedraweffect=clusrand.raweffect.*mask;
cfg=[];
cfg.yparam='maskedraweffect';
cfg.colorbar='no';
cfg.xlim=[0:0.1:1];
clf;
topoplotER(cfg,plotdata);






















