%% Using Matlab to simulate data
clear
close all

% Lets generate a single FLIM decay.
time = .05:.05:10; % Time, in nanoseconds
tao = 2.5; % Time constant in nanoseconds
F0 = 1; % Initial fluorescence

% Fill in the FLIM equaition
FLIM = F0 * exp(-time./tao);

% We can also add some noise
noisePct = .50;

FLIMNoise = FLIM + (noisePct * F0 * rand(length(time),1)') - (noisePct/2);

figure
plot(time,FLIM,time,FLIMNoise)

% We can also look at the data on a log plot to see if it's a straight line
figure
plot(time,log(FLIM),time,log(FLIMNoise))
xlabel('Time (ns)')
ylabel('Log Intensity')

%%
% Now lets fit the data using the non-linear least square curve fit.
% First, lets define a general function

% First we define a function that we want to fit.  This is in the form of
% the solved variables, and then the independent variable, followed by the
% function.  Be sure to tell MATLAB this is a function with the @ symblol.
FLIMfun = @(taoF,time)exp(-time./taoF);

% Choose an inital guess to fit.
taoGuess = mean(FLIMNoise);

% The lsqcurvefit minimizes the residual error of the fit
taoFit = lsqcurvefit(FLIMfun,taoGuess,time,FLIMNoise)

% We can then plot out how well we did.
FLIMFit = FLIMfun(taoFit,time);

figure
plot(time,FLIMNoise,time,FLIMFit)
legend('Original','Fit')


%% We can see how well we fit as a function of the random noise
clear
close all

time = .05:.05:10; % Time, in nanoseconds
tao = 2.5; % Time constant in nanoseconds
F0 = 1; % Initial fluorescence

FLIM = F0 * exp(-time./tao);
FLIMfun = @(taoF,time)exp(-time./taoF);

% linspace will make an array of equally spaced variables between two
% limits
noisePct = linspace(.05,2,50);

for j = 1:50
    FLIMNoise = FLIM + (noisePct(j) * F0 * rand(length(time),1)') - (noisePct(j)/2);
    taoGuess = mean(FLIMNoise);
    taoFit = lsqcurvefit(FLIMfun,taoGuess,time,FLIMNoise);
    
    taoDiff(j) = abs(tao-taoFit)/tao;
end

figure
plot(noisePct,taoDiff)

% Look at the data at breakdown
noisePctBD = 1.2;
FLIMNoiseBD = FLIM + (noisePctBD * F0 * rand(length(time),1)') - (noisePctBD/2);
figure
plot(time,FLIM,time,FLIMNoiseBD)


%% We can also extend this to two variables.  First, generate our data
clear
close all

% Lets generate a single FLIM decay.
time = .05:.05:10; % Time, in nanoseconds
tao1 = 2.5; % Time constant in nanoseconds
F01 = .3; % Initial fluorescence
tao2 = 0.8;
F02 = .7;

% Fill in the FLIM equaition
FLIM = F01 * exp(-time./tao1) + F02 * exp(-time./tao2);

% We can also add some noise
noisePct = .05;

FLIMNoise = FLIM + (noisePct * (F01+F02) * rand(length(time),1)') - (noisePct/2);

figure
plot(time,FLIM,time,FLIMNoise)

% Again, lets try the log plot
figure
plot(time,log(FLIMNoise))


%% Now lets fit the noisy data
% Lets first try fitting to one constant.
FLIMfun1 = @(taoF,time)exp(-time./taoF);

% Choose an inital guess to fit.
taoGuess1 = mean(FLIMNoise);

% The lsqcurvefit minimizes the residual error of the fit
taoFit1 = lsqcurvefit(FLIMfun1,taoGuess1,time,FLIMNoise)


% We can then plot out how well we did.
FLIMFit = FLIMfun1(taoFit1,time);

figure
plot(time,FLIMNoise,time,FLIMFit)
legend('Original','Fit')

% Lets look at how we missed
resid = FLIMFit - FLIMNoise;

figure
plot(time,resid)

%% Now lets fit to two components
% This time our function needs 4 constants to fit.  The amplitudes and
% time constants of each component.
FLIMfun = @(fitC,time)(fitC(1)*exp(-time./fitC(2)) + fitC(3)*exp(-time./fitC(4)));

% We need to guess on all 4.  They go - Amplitude of component 1, tao of
% component 1, amplitude of component 2, tao of component2.
fitCGuess = [.4 2 .6 3];

% The lsqcurvefit minimizes the residual error of the fit
taoFit = lsqcurvefit(FLIMfun,fitCGuess,time,FLIMNoise)


% We can then plot out how well we did.
FLIMFit2 = FLIMfun(taoFit,time);

figure
plot(time,FLIMNoise,time,FLIMFit2)
legend('Original','Fit')

% And we can calculate out the residuals
resid2 = FLIMFit2 - FLIMNoise;

figure
plot(time,resid2)

