clear all
close all
%bring in data and calibrate
%workpath
d=dir('*.dat');


alldata=[]
whole = {};
count_me = 5;
for  A=2: length(d)
%read in data
data=importintsphere([d(A).folder '\' d(A).name]);
    %extract data from raw UBAT file, BL is in column 4 and time or
    %recording number is in column 2
    BL  =double(data(:,4));
    RecNumb= double(data(:,2));
    %calibrate
    start=RecNumb-RecNumb(1);
    time_all=[1:1/60:start(end)];% represent timevector of data, at natural sampling rate
%1 start at 1, 60 samples per section

%take out first 10 columns as manual instructions (gain is three, double check)
    a=double(data(:,11:70)) .* double(data(:,3));
    %storing data 
    for xx = 1:size(data,1)
       if xx == 1
            fullresBL= a(xx,:)';
        else
            fullresBL= [fullresBL; a(xx,:)'];
       end

    end
    %time adjust to seconds
    newtime=[1:length(fullresBL)]*16.67;


%different stipulations that can be added:
[pks,locs,width,proms]=findpeaks(fullresBL,1:length(fullresBL), 'MinPeakProminence', 2, 'Annotate', 'extents')
%signal thresholding, find peak intensity and make a vector
     if max(pks)<1.5e7 % threshold
         continue
     else
[maxpeakintensity,xxx] = max(pks);
maxpeakloc = find(fullresBL==maxpeakintensity);
     end

plot(fullresBL,'k') 
hold on
plot(maxpeakloc(1),maxpeakintensity(1),'Marker','*','Color','g')
ten_thresh=maxpeakintensity(1)*0.01;
%yline(ten_thresh)

%actualtime=fullresBL(:,2)*16.67
plot(fullresBL(:,1),newtime,'k')

%e-fold time
efold_thresh=maxpeakintensity(1)*0.36788
[dontuse,ind]=min(abs(fullresBL(maxpeakloc(1):end)-efold_thresh))

if fullresBL(maxpeakloc(1)+ind)<efold_thresh
    time_efold=interp1([fullresBL(maxpeakloc(1)+(ind-1)),fullresBL(maxpeakloc(1)+ind)],[maxpeakloc(1)+(ind-1),maxpeakloc(1)+ind],efold_thresh,'linear');
elseif fullresBL(maxpeakloc(1)+ind)>efold_thresh
    time_efold=interp1([fullresBL(maxpeakloc(1)+(ind+1)),fullresBL(maxpeakloc(1)+ind)],[maxpeakloc(1)+(ind+1),maxpeakloc(1)+ind],efold_thresh,'linear');
end
efold_decay=(time_efold-maxpeakloc)*16.67;


% plot(fusi2(2),'k') 
% hold on
% plot(maxpeakloc(1),maxpeakintensity(1),'Marker','*','Color','g')
% ten_thresh=maxpeakintensity(1)*0.01;
% %yline(ten_thresh)
% efold_thresh=maxpeakintensity(1)*0.36788
% [dontuse,ind]=min(abs(fusi2(2)-efold_thresh))
% 
% if fusi2(2)(ind)<efold_thresh
%     time_efold=interp1([fusi2(2)(ind-1),fusi2(2)(ind)],[ind-1,ind],efold_thresh,'linear');
% elseif fusi2(2)(ind)>efold_thresh
%     time_efold=interp1([fusi2(2)(ind+1),fusi2(2)(ind)],[ind+1,ind],efold_thresh,'linear');
% end
% efold_decay=(time_efold-maxpeakloc)*16.67;


yline(efold_thresh)
%log plot
figure;semilogy(fullresBL,'.-b')

%get the slope of the decay
x1=1271
x2=1283
y1=9768000000
y2=74000000

m = (y2-y1)/(x2-x1) 

%interpolate linerarly between two points

x= [2121 2122]
v= [0 532800000]
xq= [2121.5]

figure
vq1=interp1(x,v,xq);
plot(fullresBL,'-b') 
hold on
plot(maxpeakloc(1),maxpeakintensity(1),'Marker','*','Color','g')
hold on
plot(x,v,'o',xq,vq1,':.')

risetime=(maxpeakloc-xq)/60

end

%integration
 newtime2=[1:length(peak)]./60'   
 trapz(newtime2,peak)

Q = cumtrapz(fullresBL)
intg=Q(end)

Qp = cumtrapz(peak)
intg=Qp(end)







% Assign the y-value threshold to a variable. Makes is clearer to see how and where it is used,
% as well as makes is easy to change if needed.
thres = 1;
% Extract x- and y-data for frequencies in the desired range of 1-6 Hz.
x = yourXData(yourXData>1 & yourXData<6); % x-values in the range
y = yourYData(yourXData>1 & yourXData<6)-thres; % y-values in range & with a threshold adjustment
% Set all values that are <0 (these were originally all values less than the threshold) to zero.
y(y<0) = 0;
% Use trapz function to calculate the area under the curve.
area = trapz(x, y )