Human_Amygdala_MUA_sEEG_FearVideo/code_MATLAB/Figure_2.m

function fig = Figure_2(file_path)

% TO BE INCLUDED IN THE PATH
% fieldtrip: fieldtrip-20191028 www.fieldtrip.org
% shadeErrobar: https://nl.mathworks.com/matlabcentral/fileexchange/26311-raacampbell-shadederrorbar
%
% include plotting toolbox gramm: https://github.com/piermorel/gramm
%
chs = [1 3 5 6 8 10 12]; % healthy amygdalae
chtot = 0;
for nSubject=1:9
    file_name = sprintf('Data_Subject_%.2d_Session_01.h5',nSubject);
    f = nix.File([file_path,file_name],nix.FileMode.ReadOnly);
    
    sectionSession = f.openSection('Session');
    all_trials = sectionSession.openProperty('Number of trials').values{1}.value;
    
    block = f.blocks{1};
    group_iEEG = block.openGroup('iEEG data');

    for nTrial = 1:all_trials
        dataArray_iEEG = group_iEEG.dataArrays{nTrial};
        striEEGLabels = dataArray_iEEG.dimensions{1}.labels;
        tiEEG = (0:(double(max(dataArray_iEEG.dataExtent))-1))* ...
                    dataArray_iEEG.dimensions{2}.samplingInterval ...
                    + dataArray_iEEG.dimensions{2}.offset;
        data_iEEG = dataArray_iEEG.readAllData;
        dataMacro.time{1, nTrial} = tiEEG;
        dataMacro.trial{1, nTrial} = data_iEEG;
        
    end
    dataMacro.label = striEEGLabels';
    dataMacro.fsample = 2000;   

    %% Power spectrum with fieldtrip
    
    cfg              = [];
    cfg.output       = 'pow';
    cfg.method       = 'mtmconvol';
    cfg.taper        = 'hanning';
    cfg.keeptrials   = 'yes';
    cfg.foi          = logspace(log10(1),log10(100),30);  
        cfg.foi          = linspace( 1, 100 ,50);                         % analysis 2 to 30 Hz in steps of 2 Hz
% analysis 2 to 30 Hz in steps of 2 Hz
    cfg.t_ftimwin    = 5./cfg.foi;   % length of time window = 0.5 sec
    cfg.toi          = -5:0.1:24;
    cfg.trials       = 2:2:17;
    TFR_face         = ft_freqanalysis(cfg, dataMacro);
    cfg.trials       = 1:2:17;
    TFR_land         = ft_freqanalysis(cfg, dataMacro );
    
     
    ch = 1:length(TFR_land.label);
    tt = TFR_land.time;
    ff = TFR_land.freq;
     
    pwr_face  = squeeze(nanmean(TFR_face.powspctrm(:,ch,:,:)));
    pwr_land  = squeeze(nanmean(TFR_land.powspctrm(:,ch,:,:)));
    
    psd_face  = nanmean(pwr_face,length(size(pwr_face)));
    psd_land  = nanmean(pwr_land,length(size(pwr_land)));
    PSD_face(chtot+ch,:) = psd_face;
    PSD_land(chtot+ch,:) = psd_land;
    chtot = chtot+ch(end);
    
    
end


%% plot Figure 2
clear g
%figure('Position',[100 100 800 350]);
fig = figure;

data = [PSD_face(chs,:); PSD_land(chs,:)];
cval={'Aversive' 'Neutral'};
cind= [ones(chtot/2, 1)*1; ones(chtot/2, 1)*2];
c=cval(cind);
g=gramm( 'x', ff, 'y', data, 'color',c);
custom_statfun = @(y)([10*log10(nanmean(y));
    10*log10(nanmean(y)) - nanstd(10*log10(nanmean(y)))/sqrt(8*6);
    10*log10(nanmean(y)) + nanstd(10*log10(nanmean(y)))/sqrt(8*6)]);
g.stat_summary('setylim', true, 'type', custom_statfun);
g.set_names( 'x', 'Frequency [Hz]', 'y', 'Power [dB]' );
g.draw();

% shadedErrorBar(ff,  10*log10(nanmean(PSD_face(chs,:))), ...
%                 (nanstd(10*log10(PSD_face(chs,:)))/sqrt(8*6)),'r');
% hold on;
% shadedErrorBar(ff,  10*log10(nanmean(PSD_land(chs,:))), ...
%                 (nanstd(10*log10(PSD_land(chs,:)))/sqrt(8*6)),'b');
% 
% xlabel('Frequency [Hz]')
% ylabel('Power  [dB]')
% 
% text(60,14,'Neutral','fontsize',16,'color','b')
% text(60,12,'Aversive','fontsize',16,'color','r')

print('-dtiff','-r300','fig_2.tiff')