eegmp_analysis/code/a_eegmp_calculate_erp_tfr.m

function a_eegmp_calculate_erp_tfr(rootpath, id)

%% paths
if ismac
    currentFile = mfilename('fullpath');
    [pathstr,~,~] = fileparts(currentFile);
    cd(fullfile(pathstr,'..'))
    rootpath = pwd;
    id = 'sub-011';
end

pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
pn.data_erp = fullfile(rootpath, 'data', 'erp'); mkdir(pn.data_erp);
pn.data_erf = fullfile(rootpath, 'data', 'erf'); mkdir(pn.data_erf);
pn.tools = fullfile(rootpath, '..', 'eegmp_preproc', 'tools');
    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults

%% load data_eeg

% load preprocessed eeg data_eeg
load(fullfile(pn.data_eeg, [id,'_task-xxxx_eeg_art.mat']), 'data_eeg', 'events');

%% further preprocessing

% apply notch filter

cfg = [];
cfg.dftfilter = 'yes';
data_eeg = ft_preprocessing(cfg, data_eeg);

%% CSD transform
    
csd_cfg = [];
csd_cfg.elecfile = 'standard_1005.elc';
csd_cfg.method = 'spline';
data_eeg = ft_scalpcurrentdensity(csd_cfg, data_eeg);

%% time-frequency transform using wavelets

freq_cfg = [];
freq_cfg.channel     = 'all';
freq_cfg.method      = 'wavelet';
freq_cfg.width       = 5;
freq_cfg.keeptrials  = 'yes';
freq_cfg.output      = 'pow';
freq_cfg.foi         = 1:1:40;
freq_cfg.toi         = -1:0.05:2;
tfr = ft_freqanalysis(freq_cfg, data_eeg);

% single-trial log10-transform

tfr.powspctrm = log10(tfr.powspctrm);

%% single-trial baseline-correction for ERPs

time = data_eeg.time{1};

data_eeg_bl = data_eeg;

for indTrial = 1:numel(data_eeg_bl.trial)
    curbl = squeeze(nanmean(data_eeg_bl.trial{indTrial}(:,time>-.2 & time<0),2));
    data_eeg_bl.trial{indTrial} = data_eeg_bl.trial{indTrial}-repmat(curbl,1,numel(time));
end

%% split tfr and time series data by condition, average across trials

parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
for ind_param = 1:numel(parameter)
    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
    for ind_cond = 1:numel(conds.(parameter{ind_param}))
        cfg = [];
        cfg.trials = ismember(events.(parameter{ind_param}), ...
            conds.(parameter{ind_param}){ind_cond});
        tfravg.(parameter{ind_param}){ind_cond} = ft_freqdescriptives(cfg, tfr);
        erp.(parameter{ind_param}){ind_cond} = ft_timelockanalysis(cfg, data_eeg);
        erp_bl.(parameter{ind_param}){ind_cond} = ft_timelockanalysis(cfg, data_eeg_bl);
    end
end

erf = tfravg; clear tfravg; 

%% save data

save(fullfile(pn.data_erp, [id, '_erp.mat']), 'erp', 'conds');
save(fullfile(pn.data_erp, [id, '_erp_bl.mat']), 'erp_bl', 'conds');
save(fullfile(pn.data_erf, [id, '_erf.mat']), 'erf', 'conds');