function fig = Figure_3( strNIXFolderPath )

% Figure_3.m reproduces Figure 3 in the dataset publication
% fig = Figure_3(strNIXFolderPath) reproduces Figure 3
% strNIXFolderPath is the path of the folder with NIX files
% fig is the figure handle
%
% Add the toolbox 'gramm' to the MATLAB path
% The toolbox can be downloaded from https://github.com/piermorel/gramm
%
% Add the toolbox 'FieldTrip' to the MATLAB path
% The toolbox can be downloaded from http://www.fieldtriptoolbox.org/download

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([strNIXFolderPath  filesep 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');
    
    dataMacro = [];
    for nTrial = 1:all_trials
        dataArray_iEEG = group_iEEG.dataArrays{nTrial};
        striEEGLabels = dataArray_iEEG.dimensions{1}.labels;
        if(length(dataArray_iEEG.dataExtent)==2)
            tiEEG = (0:(double(dataArray_iEEG.dataExtent(2))-1))*dataArray_iEEG.dimensions{2}.samplingInterval+dataArray_iEEG.dimensions{2}.offset;
        else
            tiEEG = (0:(double(dataArray_iEEG.dataExtent(1))-1))*dataArray_iEEG.dimensions{2}.samplingInterval+dataArray_iEEG.dimensions{2}.offset;
        end 
        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);
    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 3
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();

% Save figure
print('-dtiff','-r300','fig_3.tiff')

end