     
set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultTextFontSize', 10);
set(0,'defaultAxesFontSize',10);
set(0,'DefaultFigureColormap',jet)
set(0,'defaultfigurecolor',[1 1 1])
set(0,'DefaultAxesTitleFontWeight', 'normal')
set(0, 'defaulttextinterpreter', 'none') % or 'tex'

%% Load the data
clear all;clc
filename = dir('*_OXT_Arti_ripples.mat');

for ii = 2:7
    %% load data
    basename = filename(ii).name;
    load(basename);
    %%
    fs = 1e4;
    Vscale = 1/(2.0^15)*10;
    LaserShutter = double(Sp2.LaserShutter)*Vscale;
    % Define your threshold
    threshold = 2.5;                                                        
    % Calculate the difference between consecutive elements
    diff_signal = diff(LaserShutter > threshold);
    % Find the start times (where diff_signal is 1)
    start_times = find(diff_signal == 1) + 1;
    % Find the end times (where diff_signal is -1)
    end_times = find(diff_signal == -1);
    % Display the start and end times
    disp('Start times:');
    disp(start_times);
    disp('End times:');
    disp(end_times);
    Npulses = min(length(start_times),length(end_times));
    pulse_times= [start_times(1:Npulses) end_times(1:Npulses)+300 end_times(1:Npulses)-start_times(1:Npulses)]/10000 ;
    %%
    OXT_EM = double(Sp2.OXT_Em)*Vscale;
    peakLocation = round(0.5*(pulse_times(:,1)+pulse_times(:,2))*10000);
    samplesToRemove = 1000;
    for i = 1:length(pulse_times)
        % Calculate the range of samples to remove
        rangeStart = max(1, peakLocation(i) - samplesToRemove);
        rangeEnd = min(length(OXT_EM), peakLocation(i) + round(1.5*samplesToRemove));
        
        % Remove the samples within the defined range
        OXT_EM(rangeStart:rangeEnd) = NaN;
    end
    OXT_EM = fillmissing(OXT_EM, 'nearest');
    
%     plot(downsample(OXT_EM,10),'k');xlim([0 inf]);
    
    %%
    p = [];
    p = polyfit(Sp2.times,OXT_EM,5);
    
    Ybaseline = polyval(p,Sp2.times);
    figure
    plot(downsample(Sp2.times,10),downsample(OXT_EM - Ybaseline,10)); hold on;
    dFoverF = smoothdata((OXT_EM - Ybaseline)./Ybaseline,'movmean',10);
    hold on;
    plot(downsample(Sp2.times,10), downsample(dFoverF,10),'r');xlim([0 inf]);
    
    pre = 5;
    post = 5;
    eventAveragedDF = []; % Initialize array for event-averaged dFoverF
    Ns = (round(length(dFoverF)/100))*100;
    dFoverF_matrix = reshape(dFoverF(1:Ns), 100, []);
    dFoverF_ds = mean(dFoverF_matrix, 1);
    fs_down = 100;
    windowSize = fs_down; % Define the window size around each event time
    dFoverFWindow = [];
    stim_locs = round(0.5*(pulse_times(:,1)+pulse_times(:,2))*100);
    for i = 1:length(stim_locs)
        eventTime = stim_locs(i); % Get the current event time
        
        % Extract dFoverF values within the window around the event time
        windowStart = max(1, eventTime - pre*windowSize);
        windowEnd = min(length(dFoverF_ds), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds)
            dFoverFWindow(i,:) = dFoverF_ds(windowStart:windowEnd) - mean(dFoverF_ds(windowStart+fs_down*3:windowStart+fs_down*5));
        end
        % Calculate the average of the extracted dFoverF values
    end
    
    % plot the event-averaged dFoverF
    figure
    time_window = (-pre*windowSize:post*windowSize)/100;
    plot(time_window,smoothdata(mean(dFoverFWindow),"gaussian",20),'r');
    %%
    % Plot the event-averaged dFoverF
%     figure
%     plot((1:size(dFoverFWindow,2))/fs_down, mean(dFoverFWindow,1),'k');hold on; xlim([0 10]);
%     plot([3 3],[mean(dFoverFWindow(:,3*fs_down-5:3*fs_down+5),2)-std(dFoverFWindow(:,3*fs_down-5:3*fs_down+5),[],2) mean(dFoverFWindow(:,3*fs_down-5:3*fs_down+5),2)+std(dFoverFWindow(:,3*fs_down-5:3*fs_down+5),[],2)],'r');
%     plot([5 5],[mean(dFoverFWindow(:,5*fs_down-5:5*fs_down+5),2)-std(dFoverFWindow(:,5*fs_down-5:5*fs_down+5),[],2) mean(dFoverFWindow(:,5*fs_down-5:5*fs_down+5),2)+std(dFoverFWindow(:,5*fs_down-5:5*fs_down+5),[],2)],'r');
%     title('OXT 1mg/kg');
%     hold on;
%     xlabel('Time (s)');
%     ylabel('dF/F');
%     legend('dF/F','Stimulus','Location','northeast');
%     
%     % plot the event-averaged dFoverF
%     figure
%     time_window = (-pre*windowSize:post*windowSize)/100;
%     plot(time_window,smoothdata(mean(dFoverFWindow),"gaussian",20),'b');
    %% aritifical ripples dFoverF
    Npp = 1;
    NREM_flag = 3;
    WAKE_flag = 1;
    REM_flag = 5;
    pre = 5; % in seconds
    post = 5;  % in seconds
    smoothing_sigma = 15; % smoothing sigma in samples
    lfp_time = Sp2.lfp_time;
    pulsetimes_synched= round(0.5*(pulse_times(:,1)+pulse_times(:,2)) + lfp_time(1));
    pulsetimes_synched = pulsetimes_synched(pulsetimes_synched>0);
    SleepState = Sp2.SleepState;
    pulse_sleepstate = SleepState.idx.states(pulsetimes_synched);
    
    stim_NREM  = stim_locs(pulse_sleepstate==NREM_flag & circshift(pulse_sleepstate,5)==NREM_flag & circshift(pulse_sleepstate,-5)==NREM_flag);
    dFoverF_NREM = zeros(length(stim_NREM), pre*windowSize+post*windowSize+1);
    stim_locs = round(0.5*(pulse_times(:,1)+pulse_times(:,2))*100);
    
    for i = 1:Npp:length(stim_NREM)-1
        eventTime = stim_NREM(i); % Get the current event time
        % Extract dFoverF values within the window around the event time
        windowStart = max(1, eventTime - pre*windowSize);
        windowEnd = min(length(dFoverF_ds), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds)
            dFoverF_NREM(i,:) = dFoverF_ds(windowStart:windowEnd) - mean(dFoverF_ds(windowStart+fs_down*(pre-3):windowStart+fs_down*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    figure
    set(gcf,'Position',[100 100 700 400])
    time_window = (-pre*windowSize:post*windowSize)/100;
    plot(time_window,zscore(smoothdata(mean(dFoverF_NREM),"gaussian",smoothing_sigma)),'b');
    hold
    
    stim_WAKE = stim_locs(pulse_sleepstate==WAKE_flag & circshift(pulse_sleepstate,5)==WAKE_flag & circshift(pulse_sleepstate,-5)==WAKE_flag);
    dFoverF_WAKE = [];
    stim_locs = round(0.5*(pulse_times(:,1)+pulse_times(:,2))*100);
    for j = 1:Npp:length(stim_WAKE)
        eventTime = stim_WAKE(j); % Get the current event time
        % Extract dFoverF values within the window around the event time
        windowStart = max(1, eventTime - pre*windowSize);
        windowEnd = min(length(dFoverF_ds), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds)
            dFoverF_WAKE(j,:) = dFoverF_ds(windowStart:windowEnd) -mean(dFoverF_ds(windowStart+fs_down*(pre-3):windowStart+fs_down*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    time_window = (-pre*windowSize:post*windowSize)/100;
    plot(time_window,zscore(smoothdata(mean(dFoverF_WAKE),"gaussian",smoothing_sigma)),'k');
    
    stim_REM = stim_locs(pulse_sleepstate==REM_flag & circshift(pulse_sleepstate,5)==REM_flag & circshift(pulse_sleepstate,-5)==REM_flag);
    dFoverF_REM = [];
    stim_locs = round(0.5*(pulse_times(:,1)+pulse_times(:,2))*100);
    for k = 1:Npp:length(stim_REM)
        eventTime = stim_REM(k); % Get the current event time
        % Extract dFoverF values within the window around the event time
        windowStart = max(1, eventTime - pre*windowSize);
        windowEnd = min(length(dFoverF_ds), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds)
            dFoverF_REM(k,:) = dFoverF_ds(windowStart:windowEnd) - mean(dFoverF_ds(windowStart+fs_down*(pre-3):windowStart+fs_down*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    time_window = (-pre*windowSize:post*windowSize)/100;
    plot(time_window,zscore(smoothdata(mean(dFoverF_REM),"gaussian",smoothing_sigma)),'m');
    ylabel('dF/F (zscored)');xlabel('Time (s)');xlim([-pre post]);ylim([-3 3]);
    box off
    legend('NREM','WAKE','REM',Location='best')
    legend('boxoff') 
    yline(0,'r--')
    xline(0,'r--')
    %exportgraphics(gcf,[basename(1:45),'_dFoverF_ArtRipple_sleepstate.pdf'],'BackgroundColor','none','ContentType','vector');

    dFoverF_SleepState.NREM_ArtRipple = dFoverF_NREM;
    dFoverF_SleepState.WAKE_ArtRipple = dFoverF_WAKE;
    dFoverF_SleepState.REM_ArtRipple = dFoverF_REM;
    %% Spontaneous ripples dFoverF

    ripple_times = Sp2.ripples.timestamps(:,1:2); 
    lfp_time = Sp2.lfp_time;
    ripples_synched= round(0.5*(ripple_times(:,1)+ripple_times(:,2)));
    ripples_synched = ripples_synched(ripples_synched>0);
    SleepState = Sp2.SleepState;
    ripples_sleepstate = SleepState.idx.states(ripples_synched);
    
    ripples_locs = round((Sp2.ripples.peaks-lfp_time(1))*100);
    stim_NREM  = ripples_locs(ripples_sleepstate==NREM_flag & circshift(ripples_sleepstate,5)==NREM_flag & circshift(ripples_sleepstate,-5)==NREM_flag);
    dFoverF_NREM = zeros(length(stim_NREM), pre*windowSize+post*windowSize+1);
    stim_locs = round(0.5*(pulse_times(:,1)+pulse_times(:,2))*100);
    stim_NREM = stim_NREM((stim_NREM-pre*windowSize)>1);

    for i = 1:Npp:length(stim_NREM)-1
        eventTime = stim_NREM(i); % Get the current event time
        % Extract dFoverF values within the window around the event time
        windowStart = max(1, eventTime - pre*windowSize);
        windowEnd = min(length(dFoverF_ds), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds)
            dFoverF_NREM(i,:) = dFoverF_ds(windowStart:windowEnd) - mean(dFoverF_ds(windowStart+fs_down*(pre-3):windowStart+fs_down*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    figure
    set(gcf,'Position',[100 100 700 400])
    time_window = (-pre*windowSize:post*windowSize)/100;
    plot(time_window,zscore(smoothdata(mean(dFoverF_NREM),"gaussian",smoothing_sigma)),'b');
    hold
    
    stim_WAKE = ripples_locs(ripples_sleepstate==WAKE_flag & circshift(ripples_sleepstate,5)==WAKE_flag & circshift(ripples_sleepstate,-5)==WAKE_flag);
    dFoverF_WAKE = [];
    stim_WAKE = stim_WAKE((stim_WAKE-pre*windowSize)>1);
    for j = 1:Npp:length(stim_WAKE)
        eventTime = stim_WAKE(j); % Get the current event time
        % Extract dFoverF values within the window around the event time
        windowStart = max(1, eventTime - pre*windowSize);
        windowEnd = min(length(dFoverF_ds), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds)
            dFoverF_WAKE(j,:) = dFoverF_ds(windowStart:windowEnd) -mean(dFoverF_ds(windowStart+fs_down*(pre-3):windowStart+fs_down*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    time_window = (-pre*windowSize:post*windowSize)/100;
    plot(time_window,zscore(smoothdata(mean(dFoverF_WAKE),"gaussian",smoothing_sigma)),'k');
    
    stim_REM = ripples_locs(ripples_sleepstate==REM_flag & circshift(ripples_sleepstate,5)==REM_flag & circshift(ripples_sleepstate,-5)==REM_flag);
    dFoverF_REM = [];
    stim_REM = stim_REM((stim_REM-pre*windowSize)>1);
    for k = 1:Npp:length(stim_REM)
        eventTime = stim_REM(k); % Get the current event time
        % Extract dFoverF values within the window around the event time
        windowStart = max(1, eventTime - pre*windowSize);
        windowEnd = min(length(dFoverF_ds), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds)
            dFoverF_REM(k,:) = dFoverF_ds(windowStart:windowEnd) - mean(dFoverF_ds(windowStart+fs_down*(pre-3):windowStart+fs_down*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    time_window = (-pre*windowSize:post*windowSize)/100;
    plot(time_window,zscore(smoothdata(mean(dFoverF_REM),"gaussian",smoothing_sigma)),'m');
    ylabel('dF/F (zscored)');xlabel('Time (s)');xlim([-pre post]);ylim([-3 3]);
    box off
    legend('NREM','WAKE','REM',Location='best')
    legend('boxoff') 
    yline(0,'r--')
    xline(0,'r--')
    %exportgraphics(gcf,[basename(1:45),'dFoverF_Ripple_sleepstate.pdf'],'BackgroundColor','none','ContentType','vector');
    
    dFoverF_SleepState.NREM_Ripple = dFoverF_NREM;
    dFoverF_SleepState.WAKE_Ripple = dFoverF_WAKE;
    dFoverF_SleepState.REM_Ripple = dFoverF_REM;

    %%Save the data
    save([basename(1:45),'_dFoverF_Ripple_sleepstate.mat'],'dFoverF_SleepState'); 
     
    %% save fig
% Set figure size and position
figure('Position', [100, 100, 1200, 600]);
movmean_size = 30;
% Set plot options
options.color_area = [0 150 0]./255;    % Green theme
options.color_line = [0 150 0]./255;
options.alpha      = 0.1;
options.line_width = 2;
options.error      = 'sem';

% Define common x-axis
options.x_axis = [1:size(dFoverF_SleepState.NREM_ArtRipple,2)]./100-5;

% Plot NREM_ArtRipple_OXT
subplot(2,3,1);
plot_areaerrorbar(smoothdata(dFoverF_SleepState.NREM_ArtRipple,2,'movmean',movmean_size), options);
xlim([-5 5]); box off;
title('NREM_ArtRipple_OXT');
xline(0,'r--');
yline(0,'r--');

% Plot WAKE_ArtRipple_OXT
subplot(2,3,2);
plot_areaerrorbar(smoothdata(dFoverF_SleepState.WAKE_ArtRipple,2,'movmean',movmean_size), options);
xlim([-5 5]); box off;
title('WAKE-ArtRipple-OXT');
xline(0,'r--');

% Plot REM_ArtRipple_OXT
subplot(2,3,3);
plot_areaerrorbar(smoothdata(dFoverF_SleepState.REM_ArtRipple,2,'movmean',movmean_size), options);
xlim([-5 5]); box off;
title('REM-ArtRipple-OXT');
xline(0,'r--');
yline(0,'r--');

% Plot NREM-Ripple-OXT
subplot(2,3,4);
plot_areaerrorbar(smoothdata(dFoverF_SleepState.NREM_Ripple,2,'movmean',movmean_size), options);
xlim([-5 5]); box off;
title('NREM-Ripple-OXT');
xline(0,'r--');
yline(0,'r--');

% Plot WAKE-Ripple-OXT
subplot(2,3,5);
plot_areaerrorbar(smoothdata(dFoverF_SleepState.WAKE_Ripple,2,'movmean',movmean_size), options);
xlim([-5 5]); box off;
title('WAKE-Ripple-OXT');
xline(0,'r--');
yline(0,'r--');

% Export figure
exportgraphics(gcf,[basename(1:45),'_Ripple_OXT.pdf'],'BackgroundColor','none','ContentType','vector');

end



%% Find the peaks of the signal
[pks,locs] = findpeaks(OXT.data-Ybaseline ,'MinPeakHeight',0.1,'MinPeakDistance',2000);

figure
% plot the dFoverF and peaks on top with xline
plot(OXT.time, dFoverF,'k');xlim([500 1500]); hold on;
xline(OXT.time(locs),'r','LineWidth',1.5,'LabelHorizontalAlignment','center','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom')
xlim([500 1500]);




%% Load RAW data
Fluo_data = DemodFluo('basename','ACh_OXT011_004_231209_122948_OXT_Arti_ripples', 'mod_freq', 400);

%%

[pks,locs] = findpeaks(Fluo_data.Fluo_dFF ,'MinPeakHeight',0.004,'MinPeakDistance',2000);
figure
% plot the dFoverF and peaks on top with xline
plot(Fluo_data.time, Fluo_data.Fluo_dFF,'k'); hold on;
xline(Fluo_data.time(locs),'r','LineWidth',1.5,'LabelHorizontalAlignment','center','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom','LabelHorizontalAlignment','center','LabelVerticalAlignment','bottom','LabelOrientation','horizontal','LabelVerticalAlignment','bottom')


%%
figure;
fs = 10000;
plot(Fluo_data.time, 10*(Fluo_data.Fluo_ds-mean(Fluo_data.Fluo_ds)),'k');
hold on
tt = downsample(OXT_Em.time(2*fs+1:end-2*fs),100) -OXT_Em.time(2*fs) ;
pp = downsample(double(OXT_Em.data(2*fs+1:end-2*fs))*1/(2.0^15)*10,100);
plot(tt,pp-mean(pp),'r')


%%
plot(downsample(tt,40), downsample(Vraw,40));
hold on
xline(pulses_peak_times(1:100),'r')
