
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')

%% Add Util Files in to your folder
FolderName = cd;
AboveFolder = fullfile(FolderName, '..');
addpath(fullfile(AboveFolder,'Ach_Oxt_Utils'))

%% Get FluoData
FluoData = Get_OxtAch_Fluos();

%%
SleepState = Get_Aligned_Sleep();
%% OXT and Ach
figure
plot(FluoData.time,FluoData.Oxt_dFF,'color',[0 0.5 0])
hold on
plot(FluoData.time,FluoData.Ach_dFF,'color',[1 0 0])
xlabel('Time (s)')
figure;
plot(SleepState.WAKE,0.04,'k|','MarkerSize',3);
plot(SleepState.NREM,0.0425,'b|','MarkerSize',3);
plot(SleepState.REM,0.045,'m|','MarkerSize',3);
xlim([0,7500])
ylim([-0.04,0.05])

%%
%openfig('MySavedPlot.fig')
Ndown =100;
Ach_d = nanmean(reshape([FluoData.Ach_dFF(:); nan(mod(-numel(FluoData.Ach_dFF),Ndown),1)],Ndown,[]));
Oxt_d = nanmean(reshape([FluoData.Oxt_dFF(:); nan(mod(-numel(FluoData.Oxt_dFF),Ndown),1)],Ndown,[]));

nfft = 128;
[cxy,f] =mscohere(Ach_d(SleepState.NREM(1:end-20)),Oxt_d(SleepState.NREM(1:end-20)),hann(nfft/2),nfft/4,nfft,1);
figure
plot(f,cxy,'b')
ylabel('Coherence');
xlabel('Frequency [Hz]')
%savefig('Ach_OXT004_NREM_Coherence.fig')
hold on
[cxy,f] =mscohere(Ach_d(SleepState.WAKE(1:end-20)),Oxt_d(SleepState.WAKE(1:end-20)),hann(nfft/2),nfft/4,nfft,1);
plot(f,cxy,'k')

legend('NREM', 'AWAKE')

%%

Fluo_Categorized = Categorize_Spontaneous_Responses(FluoData);


%%
fs_fluo = Fluo_Categorized.fs;
Ach_Psd = Fluo_Categorized.dFF.WAKE.Ach_PSD;
Oxt_Psd = Fluo_Categorized.dFF.WAKE.Oxt_PSD;

f = linspace(0,fs_fluo/2, length(Ach_Psd(1,:)));
subplot(3,1,1)
plot(f,10*log10(mean(Ach_Psd,1)),'color',[1 0 0])
hold on
plot(f,10*log10(mean(Oxt_Psd,1)),'color',[0 0.6 0])
xlim([0,10]);
Ach_Psd = Fluo_Categorized.dFF.NREM.Ach_PSD;
Oxt_Psd = Fluo_Categorized.dFF.NREM.Oxt_PSD;
subplot(3,1,2)
plot(f,10*log10(mean(Ach_Psd,1)),'color',[1 0 0])
hold on
plot(f,10*log10(mean(Oxt_Psd,1)),'color',[0 0.6 0])
xlim([0,10]);
Ach_Psd = Fluo_Categorized.dFF.REM.Ach_PSD;
Oxt_Psd = Fluo_Categorized.dFF.REM.Oxt_PSD;
subplot(3,1,3)
plot(f,10*log10(mean(Ach_Psd,1)),'color',[1 0 0])
hold on
plot(f,10*log10(mean(Oxt_Psd,1)),'color',[0 0.6 0])
xlim([0,10]);
%%
id = 1;
figure
plot(Fluo_Categorized.dFF.REM.wbound{id}/100,Fluo_Categorized.dFF.REM.Ach_dFF{id},'color',[1 0 0])
hold on
plot(Fluo_Categorized.dFF.REM.wbound{id}/100,Fluo_Categorized.dFF.REM.Oxt_dFF{id},'color',[0 0.5 0])

%%
%%
id = 2;

plot(Fluo_Categorized.dFF.WAKE.Ach_dFF{id},'color',[1 0 0])
hold on
plot(Fluo_Categorized.dFF.WAKE.Oxt_dFF{id},'color',[0 0.6 0])

%% crosscorrelation between oxt and aCH
tr = [2000,7500]*100;
ACh_corr = FluoData.Ach_dFF(tr(1):tr(2));
OXT_corr = FluoData.Oxt_dFF(tr(1):tr(2));

[c,lags] = xcorr(OXT_corr,ACh_corr,50000,'normalized');
figure;
plot(lags/100,c)

%%%%%NREM sleep
ACh_corr = Fluo_Categorized.dFF.NREM.Ach_dFF{1,3};
OXT_corr = Fluo_Categorized.dFF.NREM.Oxt_dFF{1,3};

[c,lags] = xcorr(OXT_corr,ACh_corr,6000,'normalized');
figure;
plot(lags/100,c)

%%%%%REM sleep
ACh_corr = Fluo_Categorized.dFF.REM.Ach_dFF{1,3};
OXT_corr = Fluo_Categorized.dFF.REM.Oxt_dFF{1,3};

[c,lags] = xcorr(OXT_corr,ACh_corr,10000,'normalized');
figure;
plot(lags/100,c)

%%%%%WAKE
ACh_corr = Fluo_Categorized.dFF.WAKE.Ach_dFF{1,21};
OXT_corr = Fluo_Categorized.dFF.WAKE.Oxt_dFF{1,21};

[c,lags] = xcorr(OXT_corr,ACh_corr,15000,'normalized');
figure;
plot(lags/100,c)

%%%%%get Cross-frequency coherence
Cxy_all = [];f_all = [];
for ii = 1:length(Fluo_Categorized.dFF.NREM.Ach_dFF)
    fs = 10;
    ups = 10;
    dns = 1;
    ACh_corr = Fluo_Categorized.dFF.NREM.Ach_dFF{1,ii};
    OXT_corr = Fluo_Categorized.dFF.NREM.Oxt_dFF{1,ii};
    y = resample(OXT_corr,dns,ups);
    x = resample(ACh_corr,dns,ups);
    y = bandpass(ACh_corr,[0.005 9],fs);
    x = bandpass(OXT_corr,[0.005 9],fs);
    [Cxy,f] = mscohere(x,y,hamming(250),[],[],fs,'mimo');
    Cxy_all = [Cxy_all,Cxy];
    f_all = [f_all,f];
end

options.color_area = [128 193 219]./255;    % Blue theme
options.color_line = [ 52 148 186]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = f;
options.handle  = figure();
plot_areaerrorbar(Cxy_all', options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([0 max(f)]);title('MScoherence')
%%%%%%%%%%%%REM
Cxy_all = [];f_all = [];
for ii = 1:length(Fluo_Categorized.dFF.REM.Ach_dFF)
    fs = 10;
    ups = 10;
    dns = 1;
    ACh_corr = Fluo_Categorized.dFF.REM.Ach_dFF{1,ii};
    OXT_corr = Fluo_Categorized.dFF.REM.Oxt_dFF{1,ii};
    y = resample(OXT_corr,dns,ups);
    x = resample(ACh_corr,dns,ups);
    y = bandpass(ACh_corr,[0.005 9],fs);
    x = bandpass(OXT_corr,[0.005 9],fs);
    [Cxy,f] = mscohere(x,y,hamming(250),[],[],fs,'mimo');
    Cxy_all = [Cxy_all,Cxy];
    f_all = [f_all,f];
end

options.color_area = [128 193 219]./255;    % Blue theme
options.color_line = [ 52 148 186]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = f;
options.handle  = figure();
plot_areaerrorbar(Cxy_all', options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([0 max(f)]);title('MScoherence')

%%Cxy_all = [];f_all = [];
for ii = 1:length(Fluo_Categorized.dFF.WAKE.Ach_dFF)
    fs = 10;
    ups = 10;
    dns = 1;
    ACh_corr = Fluo_Categorized.dFF.WAKE.Ach_dFF{1,ii};
    OXT_corr = Fluo_Categorized.dFF.WAKE.Oxt_dFF{1,ii};
    y = resample(OXT_corr,dns,ups);
    x = resample(ACh_corr,dns,ups);
    y = bandpass(ACh_corr,[0.005 9],fs);
    x = bandpass(OXT_corr,[0.005 9],fs);
    [Cxy,f] = mscohere(x,y,hamming(250),[],[],fs,'mimo');
    Cxy_all = [Cxy_all,Cxy];
    f_all = [f_all,f];
end

options.color_area = [128 193 219]./255;    % Blue theme
options.color_line = [ 52 148 186]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = f;
options.handle  = figure();
plot_areaerrorbar(Cxy_all', options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([0 max(f)]);title('MScoherence')















