% Compares LFP power between conditions: ntrl and emo
% For more info about parameters look at:
% http://www.fieldtriptoolbox.org/reference/ft_freqanalysis

load Path_directory\Filename.mat;  % File saved after running Script 1.

cfg              = [];
cfg.output       = 'pow';
cfg.channel      = {'LFP*'};
cfg.method       = 'wavelet'; % Specifies the type of analysis: We selected Wavelet
cfg.width        =5;          % Window width was set up to 5 in our example.
cfg.gwidth       = 3;
cfg.foi          = 4:.5:80;   % Frequencies of interest in our case was 4-80 Hz
cfg.taper='hanning';
cfg.tapsmofrq=.3*cfg.foi;
cfg.toi          = -1:0.025:4; % Specify the time interval for calculation. In the presented case from -1 to 4s
cfg.t_ftimwin    = 5./cfg.foi;
cfg.keeptrials   = 'yes';
tfr_ntrl        = ft_freqanalysis(cfg,ntrl);
tfr_emo         = ft_freqanalysis(cfg,emo);

% freq baseline correction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
baseline=[-1 -.001];

cfg=[];
cfg.channel='LFPL23';
cfg.xlim=[-.2 4];
%cfg.ylim=[5 25];
cfg.zlim=[-30 30];
cfg.baseline=baseline;
cfg.baselinetype='absolute';%
subplot(2,2,1)
ft_singleplotTFR(cfg,tfr_ntrl);
title('left neutral')
subplot(2,2,2)
ft_singleplotTFR(cfg,tfr_emo);
title('left emo')
cfg.channel='LFPR23';
subplot(2,2,3)
ft_singleplotTFR(cfg,tfr_ntrl);
title('right neutral')
subplot(2,2,4)
ft_singleplotTFR(cfg,tfr_emo);
title('right emo')

cfg=[];
cfg.baseline=baseline;
cfg.baselinetype='relative';%'absolute';%
tfr_ntrl_bl=ft_freqbaseline(cfg,tfr_ntrl);
tfr_emo_bl=ft_freqbaseline(cfg,tfr_emo);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Statistical analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
stats=cell(1,2);
cfg = [];
cfg.channel          = {'LFPL23'};
cfg.neighbours       = [];
%cfg.avgoverfreq     = 'yes';
cfg.frequency        = [4 7];
cfg.latency          = [2 3];
cfg.method           = 'montecarlo';
cfg.statistic        = 'indepsamplesT';
cfg.correctm         = 'cluster'; % Define Statistical Analysis
cfg.clusteralpha     = 0.05;      % Define p-value
cfg.clusterstatistic = 'maxsum';
cfg.tail             = 0;
cfg.clustertail      = 0;
cfg.alpha            = 0.05;
cfg.correcttail    = 'alpha';
cfg.numrandomization = 1000;
cfg.design           = [zeros(1,numel(emo.trial))+1,zeros(1,numel(ntrl.trial))+2];
cfg.ivar             = 1;
stat_left = ft_freqstatistics(cfg, tfr_emo_bl,tfr_ntrl_bl);
stats{1}=stat_left;
cfg.channel          = {'LFPR23'};
stat_right = ft_freqstatistics(cfg, tfr_emo_bl, tfr_ntrl_bl);
stats{2}=stat_right;

figure;
for st=1:numel(stats)
    stat=stats{st};
    stat.mask=ones(size(stat.stat))*0.4+stat.mask;
    stat.powspctrm=stat.stat;
    subplot(numel(stats),1,st);
    cfg=[];
    cfg.maskparameter='mask';
    cfg.zlim='maxabs';
    ft_singleplotTFR(cfg,stat);
end



