
set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultTextFontSize', 12);
set(0,'defaultAxesFontSize',12);
set(0,'DefaultFigureColormap',jet)
set(0,'defaultfigurecolor',[1 1 1])
set(0,'DefaultAxesTitleFontWeight', 'normal')
set(0, 'DefaultAxesFontName', 'Arial');

%%
clear all;
name = 'LIMP014_014_211007_122021';
load([name,'ACh_LFP.mat'])
load([name,'.SleepState.states.mat'])
%
fs_ripple = RipLFPs.samplingRate;

N1 = ceil(lfp_fbtime(1)/2e4*fs_ripple);
N2 = floor(lfp_fbtime(2)/2e4*fs_ripple);
ripple_ = RipLFPs.data(N1:N2);
Theta_ = THLFPs.data(N1:N2);

%% NREM
fs1 = 100;
tstart = ceil(lfp_fbtime(1)/2e4);
tend = floor(lfp_fbtime(2)/2e4);
sleep_ints = bz_INTtoIDX(SleepState.ints,statenames ={'WAKE','','NREM','','REM'});
sleeps = sleep_ints.states;
sleep_ = round(round(interp(sleeps(tstart:end),fs1)*5)/5);
N_dff = length(dFoverF);
if N_dff< length(sleep_)
    sleep_ = sleep_(1:N_dff);
end
dFoverF_ = zeros(size(dFoverF));
dFoverF_(sleep_ == 1) = dFoverF(sleep_ == 1);%
peak_mean_diff = 20;
threshold_ = 0.05;
[pks, locs, width] = findpeaks(dFoverF_,'MinPeakProminence',threshold_,...
    'WidthReference','halfprom','MinPeakDistance',peak_mean_diff*fs1);

dF_event_mean = [];
dFoverF_ = zeros(size(dFoverF));
dFoverF_(sleep_ == 1) = dFoverF(sleep_ ==1);
dFoverF_(dFoverF_ ==0) = NaN;
locs_min =[];
count_ = 0;
for i = 1:length(locs)
    if (locs(i)-peak_mean_diff*fs1 >0) && (locs(i)+peak_mean_diff*fs1<N_dff)
        a = dFoverF_(locs(i)-peak_mean_diff*fs1:locs(i)+peak_mean_diff*fs1);
        if ~anynan(a)
            count_= count_ + 1;
            dF_event_mean = [dF_event_mean a];
            locs_min = [locs_min locs(i)];
        end
    end
end
%
options.handle     =  figure();
set(gcf,'Position',[100 100 500 250])
options.color_area = [255 0 0]./255;
options.color_line = [255 0 0]./255;
options.alpha      = 0.2;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = -peak_mean_diff:1/fs1:peak_mean_diff;
plot_areaerrorbar(dF_event_mean',options)
box off

%%
locs_insec = (locs_min/fs1 + lfp_fbtime(1)/2e4)';
[ wavAvg, lfpAvg ] = bz_eventWavelet (RipLFPs, locs_insec,'twin', [20,20],'plotWave',false);
%save([name,'_Peak_wavelet_ripChan.mat'],'wavAvg','lfpAvg')
%%
figure;
set(gcf,'Position',[100 100 500 300])
h = pcolor(downsample(wavAvg.timestamps,50),wavAvg.freqs,downsample(wavAvg.data,50)');hold on;
set(h, 'EdgeColor', 'none');
shading interp;
set(gca,'YDir','normal')
set(gca,'YScale','log');ylim([30 250]);
% ylim([wavespec.freqs(1) wavespec.freqs(end)]);
colormap('jet')
xlabel('time (ms)'); ylabel('frequency (Hz)');
colorbar
clim([40, 160])
%xline(0,'m','LineWidth',2)
yyaxis right
df = diff(wavAvg.freqs);
wavAvg2 = 0.5*(wavAvg.data(:,2:end) + wavAvg.data(:,1:end-1));
wavAvg2_band = wavAvg2.*df;
Fband = wavAvg.freqs(99)-wavAvg.freqs(87);
wavAvg_smoothed = smoothdata(mean(wavAvg2_band(:,87:99),2),'sgolay',1e3);
plot(wavAvg.timestamps,wavAvg_smoothed,'k')

fs_ripple = 1250;
a = (-mean(wavAvg_smoothed(1:5*fs_ripple)) +  mean(wavAvg_smoothed(15*fs_ripple:20*fs_ripple)))/mean(wavAvg_smoothed(1:5*fs_ripple))
%set(gcf,'renderer','Painters')
%exportgraphics(gcf,'LIMP016_004_220101_144210_AchNREM_Peak.pdf','BackgroundColor','none','ContentType','vector')
%%
window = 2048/8;
noverlap = window*3/4;
fs = 1250;
figure;
set(gcf,'Position',[100 100 500 300])
[pxx,f] = pwelch(lfpAvg.data,window,noverlap,window,fs);
plot(f,10*log10(pxx),'k');xlim([4,250]);
%semilogx(wavAvg.freqs,mean(wavAvg.data,1));xlim([30,250]);


%%
[spec_clean, pink_spectrum, mean_spectrum] = spectrogram_regress_pink( pxx, f );
loglog(f,(spec_clean),'k');xlim([30,250]);
%% AWAKE
fs1 = 100;
tstart = ceil(lfp_fbtime(1)/2e4);
tend = floor(lfp_fbtime(2)/2e4);
sleep_ints = bz_INTtoIDX(SleepState.ints);
sleeps = sleep_ints.states;
sleeps(sleeps ==3) =5;
sleeps(sleeps ==2) =3;
sleep_ = round(interp(sleeps(tstart:tend),fs1)*2)/2;
N_dff = length(dFoverF);
if N_dff< length(sleep_)
    sleep_ = sleep_(1:N_dff);
end
dFoverF_ = zeros(size(dFoverF));
dFoverF_(sleep_ == 1) = dFoverF(sleep_ == 1);
peak_mean_diff = 20;
[pks, locs, width] = findpeaks(dFoverF_,'MinPeakProminence',threshold_,...
    'WidthReference','halfprom','MinPeakDistance',peak_mean_diff*fs1);

dF_event_mean = [];
dFoverF_ = zeros(size(dFoverF));
dFoverF_(sleep_ == 1) = dFoverF(sleep_ == 1);
dFoverF_(dFoverF_ ==0) = NaN;
locs_min =[];
count_ = 0;
for i = 1:length(locs)
    if (locs(i)-peak_mean_diff*fs1 >0) && (locs(i)+peak_mean_diff*fs1<N_dff)
        a = dFoverF_(locs(i)-peak_mean_diff*fs1:locs(i)+peak_mean_diff*fs1);
        if ~anynan(a)
            count_= count_ + 1;
            dF_event_mean = [dF_event_mean a];
            locs_min = [locs_min locs(i)];
        end
    end
end
%
options.handle     =  figure();
set(gcf,'Position',[100 100 500 250])
options.color_area = [255 0 0]./255;
options.color_line = [255 0 0]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = -peak_mean_diff:1/fs1:peak_mean_diff;
plot_areaerrorbar(dF_event_mean',options)
box off
%%
locs_insec = (locs_min/fs1 + lfp_fbtime(1)/2e4)';
[ wavAvg, lfpAvg ] = bz_eventWavelet (THLFPs, locs_insec,'twin', [20,20]);
%save([name,'Through_wavelet_LFPChan.mat'],'wavAvg','lfpAvg')

%%
figure;
set(gcf,'Position',[100 100 500 300])
h = pcolor(downsample(wavAvg.timestamps,50),wavAvg.freqs,downsample(wavAvg.data,50)');hold on;
set(h, 'EdgeColor', 'none');
shading interp;
set(gca,'YDir','normal')
set(gca,'YScale','log');ylim([30 250]);
% ylim([wavespec.freqs(1) wavespec.freqs(end)]);
colormap('jet')
xlabel('time (ms)'); ylabel('frequency (Hz)');
colorbar
clim([70, 270])
%xline(0,'m','LineWidth',2)
yyaxis right
df = diff(wavAvg.freqs);
wavAvg2 = 0.5*(wavAvg.data(:,2:end) + wavAvg.data(:,1:end-1));
wavAvg2_band = wavAvg2.*df;
Fband = wavAvg.freqs(86)-wavAvg.freqs(72);
plot(wavAvg.timestamps,smoothdata(mean(wavAvg2_band(:,72:86),2),'sgolay',1e3)/Fband,'k','LineWidth',1)

%set(gcf,'renderer','Painters')
%exportgraphics(gcf,'LIMP016_004_220101_144210_AchNREM_Peak.pdf','BackgroundColor','none','ContentType','vector')
%%
window = 2048/8;
noverlap = window*3/4;
fs = 1250;
figure;
set(gcf,'Position',[100 100 500 300])
[pxx,f] = pwelch(lfpAvg.data,window,noverlap,window,fs);
plot(f,10*log10(pxx),'k');xlim([4,250]);
%semilogx(wavAvg.freqs,mean(wavAvg.data,1));xlim([30,250]);

%% REM
fs1 = 100;
tstart = ceil(lfp_fbtime(1)/2e4);
tend = floor(lfp_fbtime(2)/2e4);
sleep_ints = bz_INTtoIDX(SleepState.ints,statenames ={'WAKE','','NREM','','REM'});
sleeps = sleep_ints.states;
sleep_ = round(round(interp(sleeps(tstart:tend),fs1)*5)/5);
N_dff = length(dFoverF);
if N_dff< length(sleep_)
    sleep_ = sleep_(1:N_dff);
end
dFoverF_ = zeros(size(dFoverF));
dFoverF_(sleep_ == 5) = dFoverF(sleep_ == 5);%
peak_mean_diff = 10;
threshold_ = 0.05;
[pks, locs, width] = findpeaks(dFoverF_,'MinPeakProminence',threshold_,...
    'WidthReference','halfprom','MinPeakDistance',peak_mean_diff*fs1);
figure;findpeaks(dFoverF_,'MinPeakProminence',threshold_,...
    'WidthReference','halfprom','MinPeakDistance',peak_mean_diff*fs1);

dF_event_mean = [];
dFoverF_ = zeros(size(dFoverF));
dFoverF_(sleep_ == 5) = dFoverF(sleep_ ==5);
dFoverF_(dFoverF_ == 0) = NaN;
locs_min =[];
count_ = 0;
for i = 1:length(locs)
    if (locs(i)-peak_mean_diff*fs1 >0) && (locs(i)+peak_mean_diff*fs1<N_dff)
        a = dFoverF_(locs(i)-peak_mean_diff*fs1:locs(i)+peak_mean_diff*fs1);
        if ~anynan(a)
            count_= count_ + 1;
            dF_event_mean = [dF_event_mean a];
            locs_min = [locs_min locs(i)];
        end
    end
end
%
options.handle     =  figure();
set(gcf,'Position',[100 100 500 250])
options.color_area = [255 0 0]./255;
options.color_line = [255 0 0]./255;
options.alpha      = 0.2;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = -peak_mean_diff:1/fs1:peak_mean_diff;
plot_areaerrorbar(dF_event_mean',options)
box off

%%
locs_insec = (locs_min/fs1 + lfp_fbtime(1)/2e4)';
[ wavAvg, lfpAvg ] = bz_eventWavelet (THLFPs, locs_insec,'twin', [20,20],'plotWave',false);
%save([name,'_Peak_wavelet_ripChan.mat'],'wavAvg','lfpAvg')
%%
figure;
set(gcf,'Position',[100 100 500 300])
h = pcolor(downsample(wavAvg.timestamps,50),wavAvg.freqs,downsample(wavAvg.data,50)');hold on;
set(h, 'EdgeColor', 'none');
shading interp;
set(gca,'YDir','normal')
set(gca,'YScale','log');ylim([30 250]);
% ylim([wavespec.freqs(1) wavespec.freqs(end)]);
colormap('jet')
xlabel('time (ms)'); ylabel('frequency (Hz)');
colorbar
clim([70, 270])
%xline(0,'m','LineWidth',2)
yyaxis right
df = diff(wavAvg.freqs);
wavAvg2 = 0.5*(wavAvg.data(:,2:end) + wavAvg.data(:,1:end-1));
wavAvg2_band = wavAvg2.*df;
Fband = wavAvg.freqs(86)-wavAvg.freqs(72);
plot(wavAvg.timestamps,smoothdata(mean(wavAvg2_band(:,72:86),2),'sgolay',1e3)/Fband,'k','LineWidth',1)

%set(gcf,'renderer','Painters')
%exportgraphics(gcf,'LIMP016_004_220101_144210_AchNREM_Peak.pdf','BackgroundColor','none','ContentType','vector')
%%
window = 2048/8;
noverlap = window*3/4;
fs = 1250;
figure;
set(gcf,'Position',[100 100 500 300])
[pxx,f] = pwelch(lfpAvg.data,window,noverlap,window,fs);
plot(f,10*log10(pxx),'k');xlim([4,250]);