%%
cd('G:\My Drive\Paper_ACh_OXT\Fig2\nonNREMpackets_data')
clear all;
%% figure2B
files = dir('*.mat');
% GET ALL SPECTRUM
NREM_spectrum_all =[];
for ii = 1:length(files)
    
    load(files(ii).name)
    
    % exclude data shorter than 6 seconds
    idx = find(diff(nNREM_Data.nNREMpac_Data.ints,1,2)>6);
    idx = idx(idx<length(nNREM_Data.nNREMpac_Data.LFP));
    nREM_ACh = nNREM_Data.nNREMpac_Data.ACh(idx,1);
    nREM_LFP = nNREM_Data.nNREMpac_Data.LFP(idx,1);
    
    % get spectrum NREM packects
    THLFPs.samplingRate = 1250;
    thetarange = linspace(1,100,250);
    NREM_spectrum= {};
    for jj = 1:length(nREM_LFP)
        if ~isempty(nREM_LFP{jj,1})>0
            lfp = nREM_LFP{jj,1};
            [ss,ft,tt] = spectrogram(lfp,2*THLFPs.samplingRate,0.6*THLFPs.samplingRate,thetarange,THLFPs.samplingRate);
            [spectrogrammt,tm,fm] = MTSpectrogram(lfp,'window',1);
            bands = SpectrogramBands(spectrogrammt,fm);
        %     figure;
        %     colormap('jet')
        %     imagesc(tt,log2(thetarange),abs(ss));clim([0,8e4]);
        %     yticks([1,2,3,4,5]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
        %     set(gca,'YDir','normal');
            NREM_spectrum{jj,1} = ss;
            NREM_spectrum{jj,2} = tt;
        end
    end
    ii
    NREM_spectrum_all = [NREM_spectrum_all;NREM_spectrum ];
end
NREM_spectrum_all = NREM_spectrum_all(~all(cellfun(@isempty,NREM_spectrum_all),2),:)

%% get time-normalized 

p =  cell2mat(cellfun(@size,NREM_spectrum_all,'UniformOutput',false));
F_idx = max(p(:,2));
B = NREM_spectrum_all{find(max(p(:,2))),1}
%%
Nsize = size(NREM_spectrum_all,1);
Nf = size(NREM_spectrum_all{1},1);
NREMSpectrumNormalized = zeros([Nsize,Nf,F_idx]);
tt= linspace(0,1,F_idx);
ff = 1:Nf;
[Tq,Fq] = meshgrid(tt,ff);
for i=1:Nsize
    sp = NREM_spectrum_all{i};
    t_= linspace(0,1,size(sp,2));
    [T,F] = meshgrid(t_,ff);
    NREMSpectrumNormalized(i,:,:) =abs(interp2(T,F,sp,Tq,Fq));
end
%% Normal version
figure
set(gcf,'Position',[100 100 600 300])
subplot(2,1,1)
spPower = squeeze(mean((abs(NREMSpectrumNormalized)),1));
h = pcolor(tt,log2(thetarange),(spPower));colormap("jet");
clim(([0.2e5 1.5e5]));
set(h, 'EdgeColor', 'none'); shading interp;
set(gca,'YDir','normal');
ylabel('Frequency (Hz)'); xlabel('Time (s)');ylim(log2([10 20]));
yticks(log2([10,20]))
yt = get(gca, 'YTick');
set (gca, 'YTickLabel', 2.^yt);
box off
colorbar
subplot(2,1,2)
h = pcolor(tt,log2(thetarange),(spPower));colormap("jet");
clim(([0.2e5 2.6e5]));
set(h, 'EdgeColor', 'none');shading interp;
set(gca,'YDir','normal');
ylabel('Frequency (Hz)'); xlabel('Time (s)');ylim(log2([1 10]));
yticks(log2([1,10]))
yt = get(gca, 'YTick');
set (gca, 'YTickLabel', 2.^yt);
box off
colorbar
%exportgraphics(gcf,'Fig2B_NREMPAcket_Spectrogram.pdf','BackgroundColor','none','ContentType','vector')

%% Zscored version
figure
set(gcf,'Position',[100 100 600 300])
subplot(2,1,1)
spPower = zscore(squeeze(mean((abs(NREMSpectrumNormalized)),1)),[],2);
h = pcolor(tt,log2(thetarange),(spPower));colormap("jet");
%clim(([0.2e5 1.5e5]));
set(h, 'EdgeColor', 'none');
set(gca,'YDir','normal');
ylabel('Frequency (Hz)'); xlabel('Time (s)');ylim(log2([10 20]));
yticks(log2([10,20]))
yt = get(gca, 'YTick');
set (gca, 'YTickLabel', 2.^yt);
clim([-3,3])
box off
colorbar
subplot(2,1,2)
h = pcolor(tt,log2(thetarange),(spPower));colormap("jet");
%clim(([0.2e5 2.6e5]));
set(h, 'EdgeColor', 'none');
set(gca,'YDir','normal');
ylabel('Frequency (Hz)'); xlabel('Time (s)');ylim(log2([1 10]));
yticks(log2([1,10]))
yt = get(gca, 'YTick');
set (gca, 'YTickLabel', 2.^yt);
box off
colorbar
clim([-3,3])
exportgraphics(gcf,'Fig2B_NREMPAcket_Spectrogram_Zscored.pdf','BackgroundColor','none','ContentType','vector')
%%

figure

plot(tt,smoothdata(mean(zscore(spPower(25:50,:),[],2),1),'gaussian',200),'r')
hold on
plot(tt,smoothdata(mean(zscore(spPower(2:10,:),[],2),1),'gaussian',200),'k')
plot(tt,smoothdata(mean(zscore(spPower(75:200,:),[],2),1),'gaussian',200),'b');
%% get ACh and normalize the time
NREM_ACh_all  = {};
for ii = 1:5
    
    load(files(ii).name)
    
    % exclude data shorter than 6 seconds
    idx = find(diff(nNREM_Data.nNREMpac_Data.ints,1,2)>6);
    idx = idx(idx<length(nNREM_Data.nNREMpac_Data.LFP));
    nREM_ACh = nNREM_Data.nNREMpac_Data.ACh(idx,1);
    
    NREM_ACh_all = [NREM_ACh_all;nREM_ACh ];
end
p =  cell2mat(cellfun(@size,NREM_ACh_all,'UniformOutput',false));
F_idx = max(p(:,1));
B = NREM_ACh_all{find(max(p(:,2))),1}

Nsize = size(NREM_ACh_all,1);
NREM_AChNormalized = [];
NREM_AChNormalized_dd = [];
tt= linspace(0,1,F_idx);
for i=1:Nsize
    ach = NREM_ACh_all{i};
    t_= linspace(0,1,size(ach,1));
    NREM_AChNormalized(i,:) =interp1(t_,ach,tt);
    NREM_AChNormalized_dd(i,:) = downsample(NREM_AChNormalized(i,:),10);
end

tt_dd = downsample(tt,10);
options.color_area = [255 204 204]./255;    % Blue theme
options.color_line = [255 0 0]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = tt_dd;
% figure;plot_areaerrorbar(zscore(NREM_AChNormalized,[],2), options)
figure;plot_areaerrorbar(zscore(NREM_AChNormalized_dd,[],2), options)

%% get OXT and normalize the time
NREM_OXT_all  = {};
for ii = 6:9
    
    load(files(ii).name)
    
    % exclude data shorter than 6 seconds
    idx = find(diff(nNREM_Data.nNREMpac_Data.ints,1,2)>6);
    idx = idx(idx<length(nNREM_Data.nNREMpac_Data.LFP));
    nREM_OXT = nNREM_Data.nNREMpac_Data.ACh(idx,1);
    
    NREM_OXT_all = [NREM_OXT_all;nREM_OXT ];
end

NREM_OXT_all = NREM_OXT_all(~all(cellfun(@isempty,NREM_OXT_all),2),:)

p =  cell2mat(cellfun(@size,NREM_OXT_all,'UniformOutput',false));
F_idx = max(p(:,1));
B = NREM_OXT_all{find(max(p(:,2))),1}

Nsize = size(NREM_OXT_all,1);
NREM_OXTNormalized = [];
NREM_OXTNormalized_dd = [];
tt= linspace(0,1,F_idx);
for i=1:Nsize
    ach = NREM_OXT_all{i};
    t_= linspace(0,1,size(ach,1));
    NREM_OXTNormalized(i,:) =interp1(t_,ach,tt);
    NREM_OXTNormalized_dd(i,:) = downsample(NREM_OXTNormalized(i,:),10);
end

hold on;
tt_dd = downsample(tt,10);
options.color_area = [229 255 204]./255;    % Blue theme
options.color_line = [ 0 173 76]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = tt_dd;
% plot_areaerrorbar(zscore(NREM_OXTNormalized,[],2), options)
plot_areaerrorbar(zscore(NREM_OXTNormalized_dd,[],2), options)

exportgraphics(gcf,'Fig2C_ACh_OXT.pdf','BackgroundColor','none','ContentType','vector')
%% get NREMepisodes data
NREMepi_all = {};
NREMepi_all.ACh = [];
NREMepi_all.LFP = [];
for ii = 1:length(files)
    
    load(files(ii).name)
    
    % exclude data shorter than 6 seconds
    idx = find(diff(nNREM_Data.nNREMpac_Data.ints,1,2)>6);
    idx = idx(idx<length(nNREM_Data.nNREMpac_Data.LFP));
    nREM_ACh = nNREM_Data.nNREMpac_Data.ACh(idx,1);
    nREM_LFP = nNREM_Data.nNREMpac_Data.LFP(idx,1);
    
    % get spectrum NREM packects
    THLFPs.samplingRate = 1250;
    thetarange = linspace(1,100,250);
    NREM_spectrum= {};
    for jj = 1:length(nREM_LFP)
        if ~isempty(nREM_LFP{jj,1})>0
            lfp = nREM_LFP{jj,1};
            [ss,ft,tt] = spectrogram(lfp,2*THLFPs.samplingRate,0.6*THLFPs.samplingRate,thetarange,THLFPs.samplingRate);
            [spectrogrammt,tm,fm] = MTSpectrogram(lfp,'window',1);
            bands = SpectrogramBands(spectrogrammt,fm);
        %     figure;
        %     colormap('jet')
        %     imagesc(tt,log2(thetarange),abs(ss));clim([0,8e4]);
        %     yticks([1,2,3,4,5]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
        %     set(gca,'YDir','normal');
            NREM_spectrum{jj,1} = ss;
            NREM_spectrum{jj,2} = tt;
        end
    end
    ii
    NREM_spectrum_all = [NREM_spectrum_all;NREM_spectrum ];
end

