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');
    %%
    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');
    
%% random chosse time points
    data = Sp2.times; % replace with your 1xN data
    numPoints = length(peakLocation); % number of random time points to select
    numIterations = 1000; % number of times to select random time points
    
    % Initialize a cell array to hold the results
    dFoverF_NREM_Session = [];
    % Loop over the number of iterations
    for i = 1:numIterations
       
        % Generate random indices
        indices = randi([1 length(data)], 1, numPoints);
    
        % Select random time points and store in the cell array
        randomTimePoints = data(indices);
        randomTimePoints_aligned = randomTimePoints;

        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
        fs = 100;

        lfp_time = Sp2.lfp_time;
        randomTimePoints_synched= round(randomTimePoints_aligned + lfp_time(1));
        randomTimePoints_synched = randomTimePoints_synched(randomTimePoints_synched>0);
        SleepState = Sp2.SleepState;
        randomTimePoints_sleepstate = SleepState.idx.states(randomTimePoints_synched);
        stim_locs = randomTimePoints_synched*fs;
        %% get data of NREM
        pre = 5;
        post = 5;
        eventAveragedDF = []; % Initialize array for event-averaged dFoverF
        Ns = length(dFoverF_ds);
        windowSize = fs; % Define the window size around each event time
        dFoverFWindow = [];
        stim_NREM  = stim_locs(randomTimePoints_sleepstate==NREM_flag & circshift(randomTimePoints_sleepstate,5)==NREM_flag & circshift(randomTimePoints_sleepstate,-5)==NREM_flag);
        dFoverF_NREM = zeros(length(stim_NREM), pre*windowSize+post*windowSize+1);
        %stim_locs = round(0.5*(pulse_times_aligned(:,1)+pulse_times_aligned(:,2))*130);
        dFoverF_NREM=[];
        for k = 1:Npp:length(stim_NREM)-1
            eventTime = stim_NREM(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_NREM(k,:) = dFoverF_ds(windowStart:windowEnd) - mean(dFoverF_ds(windowStart+fs*(pre-3):windowStart+fs*pre));
            end
        end

        dFoverF_NREM_Session(i,:) = mean(dFoverF_NREM,1);
    end
    dFoverF_NREM_all{ii} = zscore(mean(dFoverF_NREM_Session,1));
end


