
% Load your input data
Input = PhotometryData20220822T203705; %example data
InputBlue = Input(find(Input(:,3)==1),:);
InputGreen = Input(find(Input(:,3)==4),:);

InputBlueTime = Input(find(Input(:,3)==1),2);
InputGreenTime = Input(find(Input(:,3)==4),2);

% Plotting Green Channel Data
figure
plot(InputGreenTime, InputGreen(:,10))
hold on
plot(InputGreenTime, InputGreen(:,11))

% Loading stimulation data
StimtimeYFP = Stimulation20220822T203704;

% Isolate stimulation data for a specific condition
LEDTrunc73 = StimtimeYFP(find(StimtimeYFP(:,2)==73),:);
LEDTrunc73idx = intersect(find(StimtimeYFP(:,2)==73), find(StimtimeYFP(:,1)==3));
LEDTrunc73 = LEDTrunc73(find(LEDTrunc73(:,1)==3),:);

% Detect changes in stimulation state
count = 0;
TimeStimYFP = [];
for i = 2:length(LEDTrunc73)
    if LEDTrunc73(i-1,10) == 0 & LEDTrunc73(i,10) >= 1
        LEDTrunc73(i,10) = 4;
        count = count + 1;
        TimeStimYFP(count,1) = LEDTrunc73(i,3);
    end
end

% Find closest previous time points in InputGreen
maxvaluenegidx = zeros(length(TimeStimYFP(:,1)), 1);
for i = 1:length(TimeStimYFP(:,1))
    diffwithphoto = InputGreenTime - TimeStimYFP(i,1);
    maxvaluenegidx(i,1) = max(find(diffwithphoto(:,1) <= 0));
end

% Calculate baseline using median
baselineYFPVS = zeros(length(maxvaluenegidx(:,1)), 1);
baselineYFPDLS = zeros(length(maxvaluenegidx(:,1)), 1);
goback = 10;
for i = 1:length(maxvaluenegidx(:,1))
    if maxvaluenegidx(i,1) - goback > 0
        baselineYFPVS(i,1) = median(InputGreen(maxvaluenegidx(i,1)-goback:maxvaluenegidx(i,1)-1,10));
        baselineYFPDLS(i,1) = median(InputGreen(maxvaluenegidx(i,1)-goback:maxvaluenegidx(i,1)-1,11));
    end
end

% Delta F over F calculation
goforward = 250;
baselinelength = 50;
deltaFYFPVS = [];
deltaFYFPDLS = [];
for i = 1:length(maxvaluenegidx(:,1))
    if maxvaluenegidx(i,1)+goforward <= length(InputGreen(:,1))
        rawSignalVS = InputGreen(maxvaluenegidx(i,1)-baselinelength:maxvaluenegidx(i,1)+goforward,10);
        rawSignalDLS = InputGreen(maxvaluenegidx(i,1)-baselinelength:maxvaluenegidx(i,1)+goforward,11);
        
        % Applying custom moving average function
        signalYFPVS = simpleMovingAverage(rawSignalVS, 10);  % Window size of 10
        signalYFPDLS = simpleMovingAverage(rawSignalDLS, 10);  % Window size of 10

        deltaFYFPVS(i,:) = (signalYFPVS - baselineYFPVS(i,1)) / baselineYFPVS(i,1);
        deltaFYFPDLS(i,:) = (signalYFPDLS - baselineYFPDLS(i,1)) / baselineYFPDLS(i,1);
    end
end

% Visualization of Delta F/F
figure
imagesc(deltaFYFPDLS)
title('YFP2 VTA-->DLS')
colormap(gca, jet(256))
colorbar(gca)
caxis([-0.01 0.07])

figure
imagesc(deltaFYFPVS)
title('YFP2 VTA-->VS')
colormap(gca, jet(256))
colorbar(gca)
caxis([-0.01 0.05])

% Additional plotting
figure
plot(InputGreenTime, InputGreen(:,10))
hold on
plot(InputGreenTime, InputGreen(:,11))
for i = 1:length(TimeStimYFP(:,1))
    xline(InputGreen(maxvaluenegidx(i,1),2))
end

figure
plot(median(deltaFYFPDLS), 'r')
hold on
plot(median(deltaFYFPVS), 'b')
ylim([-0.015 0.1])
