addpath('P:\3013074.03\scripts\data_scripts\EyeTracker')
fname = 'sub_22_ses1_Samples.txt';
fpath = 'P:\3013074.03\Test_data\physiological_data\eyetracker\raw\sub_22'

%1. Convert from .txt to .mat
samplerate = 50;
a_geteyetracker(fname,fpath);
load('P:\3013074.03\Test_data\physiological_data\eyetracker\raw\sub_01\eyetrackdat.mat');

%2a. interpolate holes using Linda's script
interpolated = interpoll(eyetrackdat.diameter{1,2},samplerate);

%2b. Set up data struct to use erno's interpolate script
a_setuppupdatstruct(1,1,2,fpath)
interpol(); %manualy select the resulting file and save results.



%3. plot the result
figure()
subplot(1,3,1)
plot(eyetrackdat.diameter{1,2});
title('Raw diameter');
subplot(1,3,2)
plot(interpolated);
title('Interpolated diameter Linda');
hold on
shocktimes = find(eyetrackdat.prescodes{1,2} == 100);
correctedtimes = [shocktimes(1)]
for i = 1:length(shocktimes)
    j = shocktimes(i);
    if((j -correctedtimes(end))>400)
        correctedtimes = [correctedtimes j]
    end
end
% for i = correctedtimes
%     plot([i,i],[10,45], 'r');
% end
subplot(1,3,3)
load('pupil_day1_2.mat');
plot(interpoldat.fullinterpoldat{1})
% hold on
% correctedtimes = [shocktimes(1)]
% for i = correctedtimes
%     plot([i,i],[10,45], 'r');
% end
title('Ernos interpolation');
% 
% for i = 1:size(shocktimes)
%     j = shocktimes(i);
%     figure()
%     subplot(1,2,1)
%     beginT = j - samplerate*1;
%     endT = j + samplerate*4;
%     plot(interpolated(beginT:endT));
%     prioraverage = mean(interpolated(beginT:j));
%     hold on
%     %plot([j,j],[10,45], 'r');
%     subplot(1,2,2)
%     plot((interpolated(beginT:endT)/prioraverage))
%     hold on
%     %plot([j,j],[10,45], 'r');
% end

figure()
title('Interpolated data - responses to shocks');
for i = 1:length(correctedtimes)
    j = correctedtimes(i);
    beginT = j - samplerate*1;
    endT = j + samplerate*4;
    prioraverage = mean(interpolated(beginT:j));
    subplot(4,5,i)
    plot((interpolated(beginT:endT)))
    
    %plot([j,j],[10,45], 'r');
end


figure()
title('Raw data - responses to shocks');
for i = 1:length(correctedtimes)
    j = correctedtimes(i);
    beginT = j - samplerate*1;
    endT = j + samplerate*4;
    prioraverage = mean(eyetrackdat.diameter{1,2}(beginT:j));
    subplot(4,5,i)
    plot((eyetrackdat.diameter{1,2}(beginT:endT)))
    
    %plot([j,j],[10,45], 'r');
end

figure()
title('Raw data - responses to shocks'); %Doesnt work yet becuase this contains the full time line
load('C:\Users\MAXHOU\Documents\MATLAB\Counterconditioning\raw\sub_01\pupil_day1_2.mat');

collect = [];
for i = 1:length(correctedtimes)
    j = correctedtimes(i);
    beginT = j - samplerate*1;
    endT = j + samplerate*4;
    prioraverage = mean(interpoldat.fullinterpoldat{1}(beginT:j));
    subplot(4,5,i)
    plot((interpoldat.fullinterpoldat{1}(beginT:endT)))
    ylim([20 45])
    %plot([j,j],[10,45], 'r');
    collect = [collect;[interpoldat.fullinterpoldat{1}(beginT:endT)/prioraverage]'];
end

%Blazej implemented an additonal step of filtering at 0.05 to 4 Hz bandpass
%and z-scoring this.

Fs          =   50; %Sampling frequency
f1          =   0.05; % Cut-off bandpass Hz
f2          =   5;% Cut-off bandpass Hz
Wn          =   [f1/(Fs/2) f2/(Fs/2)];
[b,a]       =   butter(2, Wn, 'bandpass');
PDmean      =   mean(interpoldat.fullinterpoldat{1});
%Filter the data
PD_filter   =   filtfilt(b,a,interpoldat.fullinterpoldat{1});
PD_filter   =   PD_filter + PDmean;


figure()
plot(PD_filter);
figure()
title('Filtered interpolated data'); %Doesnt work yet becuase this contains the full time line
collect = [];
for i = 1:length(correctedtimes)
    j = correctedtimes(i);
    beginT = j - samplerate*1;
    endT = j + samplerate*4;
    prioraverage = mean(PD_filter(beginT:j));
    subplot(4,5,i)
    plot((PD_filter(beginT:endT)))
    ylim([20 45])
    collect = [collect;[PD_filter(beginT:endT)/prioraverage]'];
end



Animaltimes = find(eyetrackdat.prescodes{1,2} == 4);
correctedtimesA = [Animaltimes(1)]
for i = 1:length(Animaltimes)
    j = Animaltimes(i);
    if((j -correctedtimesA(end))>400)
        correctedtimesA = [correctedtimesA j]
    end
end
Tooltimes = find(eyetrackdat.prescodes{1,2} == 5);
correctedtimesT = [Tooltimes(1)]
for i = 1:length(Tooltimes)
    j = Tooltimes(i);
    if((j -correctedtimesT(end))>400)
        correctedtimesT = [correctedtimesT j]
    end
end


figure()
title('Filtered interpolated data - TOOLS'); %Doesnt work yet becuase this contains the full time line
for i = 1:length(correctedtimesT)
    j = correctedtimesT(i);
    beginT = j - samplerate*1;
    endT = j + samplerate*4;
    prioraverage = mean(PD_filter(beginT:j));
    subplot(8,5,i)
    plot((PD_filter(beginT:endT)))
    ylim([20 35])
end
figure()
title('Filtered interpolated data - TOOLS NORMALIZED'); %Doesnt work yet becuase this contains the full time line
collect = [];
for i = 1:length(correctedtimesT)
    j = correctedtimesT(i);
    beginT = j - samplerate*1;
    endT = j + samplerate*4;
    prioraverage = mean(PD_filter(beginT:j));
    subplot(8,5,i)
    plot((PD_filter(beginT:endT)/prioraverage))
    ylim([0.5 1.5])
    collect = [collect;[PD_filter(beginT:endT)/prioraverage]'];
end


figure()
title('Filtered interpolated data - ANIMALS'); %Doesnt work yet becuase this contains the full time line
collectionT = [];
collect = [];
for i = 1:length(correctedtimesA)
    j = correctedtimesA(i);
    beginT = j - samplerate*1;
    endT = j + samplerate*4;
    prioraverage = mean(PD_filter(beginT:j));
    subplot(8,5,i)
    plot((PD_filter(beginT:endT)))
    collectionT = [collectionT (mean(PD_filter(beginT:endT))/prioraverage)]
    ylim([20 35])
    collect = [collect;[PD_filter(beginT:endT)/prioraverage]'];
end
figure()
title('Filtered interpolated data - ANIMALS NORMALIZED'); %Doesnt work yet becuase this contains the full time line
%CS+ for subject one
collectionA = [];
for i = 1:length(correctedtimesA)
    j = correctedtimesA(i);
    beginT = j - samplerate*1;
    endT = j + samplerate*4;
    prioraverage = mean(PD_filter(beginT:j));
    subplot(8,5,i)
    plot((PD_filter(beginT:endT)/prioraverage))
    collectionA = [collectionA (mean(PD_filter(beginT:endT))/prioraverage)]
    ylim([0.5 1.5])
    
end