addpath(genpath('humans/functions/external/fooof_mat-main'))
addpath('humans/functions')


human_datadir = 'humans/data';

%% Fig 3b

% Parameters for time-frequency spectrogram
frex = 5:5:200;
baselinetime = [-0.5 -0.1];
srLFP = 2048;

% ER10 macro1
load(fullfile(human_datadir, 'macrocontacts/ER10_macro1_HFO_snippets.mat'))

figure(1), clf
subplot(121)
plot(tt, mean(snippets,1))
xlim ([-0.1 0.1])
axis square

subplot(122)
TF = time_freq_wavelet(snippets', frex, srLFP, tt, 1, baselinetime, 1);
xlim ([-0.1 0.1])
clim([0 25])
axis square


% ER5 macro1
load(fullfile(human_datadir, 'macrocontacts/ER05_macro1_HFO_snippets.mat'))

figure(2), clf
subplot(121)
plot(tt, mean(snippets,1))
xlim ([-0.1 0.1])
axis square

subplot(122)
TF = time_freq_wavelet(snippets', frex, srLFP, tt, 1, baselinetime, 1);
xlim ([-0.1 0.1])
clim([0 25])
axis square

%% Fig 3c
% Parameters for PSD
freqs = 1:1:200;

% ER10 macro1
load(fullfile(human_datadir, 'macrocontacts/ER10_macro1_HFO_snippets.mat'))

[powspctrm] = pwelch(snippets', 100, [], freqs, srLFP);

corrected = fooof_correct_powspctrm(powspctrm, freqs);

figure(3), clf
plot(freqs, mean(log10(powspctrm),2))
hold on
plot(freqs, mean(corrected,2))
xlim('Frequency (Hz)')
ylim('Log(power)')

% ER5 macro1
load(fullfile(human_datadir, 'macrocontacts/ER05_macro1_HFO_snippets.mat'))

[powspctrm] = pwelch(snippets', 100, [], freqs, srLFP);

corrected = fooof_correct_powspctrm(powspctrm, freqs);

figure(4),clf
plot(freqs, mean(log10(powspctrm),2))
hold on
plot(freqs, mean(corrected,2))
xlim('Frequency (Hz)')
ylim('Log(power)')
%% Fig 3d Histogram
load('Figure 3/AllChannels_PSD.mat')

micro = logical([tbl.micro]);
tbl_macro = tbl(~micro,:);

rip_freqs = tbl_macro.Peak_Frequency;

figure (5), clf
histogram(rip_freqs, 'BinWidth',1)
hold on
xline(60, 'r--')
title('Macro')
xlabel('Frequency (Hz)')
ylabel('Channel counts')

