Dear -
My user name is: pyff
I would like to submit to the GNU Octave Repository a piece of my work that
I'm wrote for the analysis of luminescence exponential decay signal.
In the present script, the function leasqr and expfit are used to analyze coma
separated data. Plots and results files are created and stored in a "Results"
folder created in the working directory.
Please let me know if you need more informations or if there is some problem
With my best regards
Filippo Piffaretti
_____________________________________
Filippo Piffaretti
PhD Student, EPFL-LPAS
Station 6
CH H5 534 (Bâtiment CH)
1015 Lausanne
SWITZERLAND
**************************************
Work: +4121/693 76 23
Labo: +4121/693 36 12
Work Fax: +4121/693 51 45
Mobile: +4176/399 37 59
E-mail: [email protected]
______________________________________
## Filippo Piffaretti, EPFL 2009
##
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public
## License as published by the Free Software Foundation;
## either version 2, or (at your option) any later version.
##
## Octave is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied
## warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
## PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public
## License along with Octave; see the file COPYING. If not,
## write to the Free Software Foundation, 59 Temple Place -
## Suite 330, Boston, MA 02111-1307, USA.
## usage: decay_lumi (param, indata, deg, pin)
##
## Small script that will help to analyse exponential decay
## phenomena. This script will coordinate diffenret type of
## approved program (Data retrivering, prony Method, partial
## differentiation, Levenberg-Marquart non-linear fitting
## algorithm.
## The results of the calculation will be displayed with plots
## and the results will be saved in different text documents in
## "Results" folder.
## Author: Filippo Piffaretti <[email protected]>
## Maintainer: Filippo Piffaretti <[email protected]>
## Created: November 2009
## Version: 0.6
## Keywords: Exponential fit, luminescence
function decay_lumi (param, indata, deg, pin)
ds = 0.001; %differential step
stol= 0.0001; niter=50;
minstep = [0; 0];
maxstep = [0.8; 0.8];
offset = 0;
options = [minstep, maxstep];
%##################################################
%### Argument processing ###
%##################################################
if (nargin < 3) % Check the number of arguments give by the user
%-----------------------------------------------
printf("\t Read the 'ReadMe.txt' file to set correctely the function parameters\n");
else
%-------------------------------------------------
%************************************************
%* Start to treat the parameters and the data *
%************************************************
for i = 1:nargin %Handling the execution mode chose
if (i==1) %---------------------------------
switch(param)
case{"Silent" "silent" "S" "s"}
param = 0; %silent
case{"Verbose" "verbose" "V" "v"}
param = 1; %verbose mode
case{"Print" "print" "P" "p"}
param = 2; %verbose mode
otherwise
printf("\t Incorrect setting of the first parameter %s ", param);
param = 0;
endswitch
endif
%---------------------------------------------
if (i==2)
if(isdir(indata) == true) %Check if "indata" is a directoy or a single file
%-----------------------------------------
files=readdir(indata); %Create files array
confirm_recursive_rmdir (0); %Do not ask for confirmation to cancel folder recursively
if (indata(length(indata)) == "/") %Test the folder parameter entred, is it ending with "/"?
rmdir (strcat (indata, "Results"),'s');
FileFolder = indata; %Set the file folder
else
rmdir (strcat (indata, "/Results"),'s');
FileFolder = strcat (indata, "/");
endif
mkdir (strcat (FileFolder, "Results")); %Create the directory for the results
else %If "indata" is a file
%---------------------
if(fopen(indata,'r') == -1) %Check if the input is an existing data
printf("\t !!ERROR: The file name entred does not exist!! \n");
else
if (rindex (indata, "/") != 0)
FileFolder = substr (indata, 1, rindex (indata, "/"));
files = cellstr(indata(rindex (indata, "/")+1:end));
else
FileFolder = "./";
files= cellstr(indata); %Put the file name to be treated in an array
endif
if (!isdir(strcat (FileFolder, "Results"))) %Create the folder only if it has not yet been created
mkdir (strcat (FileFolder, "Results"));
endif
endif
endif
endif
%---------------------------------------------
if (i == 4) % Transform the units of the entred parameters
for w = 1:(length(pin))
if (mod(w,2) == 0) % (Starting point entred in: Volt; Seconds) ...
pin(w)=1/pin(w); % to the parameters needed for the algorithm ...
endif % (Volt; 1/Second).
endfor
pin = [offset; pin]
if (param==1) %If verbose send to stdout
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% Print the informations to the SdtOut
printf("\t User estimated parameters: \n\n\t Offset %d Volt\n", pin(1));
for w = 2:(length(pin))
if (mod(w,2))
printf("\t Lifetime: %.2f µs\n\n", 1/(pin(w)*1e6));
else
printf("\t Pre-exp: %.2f Volt\n", pin(w));
endif
endfor
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
endif
endif
%--------------------------------------------------------
endfor
%********************************************************
%##########################################################
%### Loop to analyse the file or the folder files ###
%##########################################################
for file_n = 1:rows(files)
lastwarn ("ok"); %Reset last warning message to zero (Singular Matrix)
FileName = cellstr(split(files{file_n,1},".")); %File name of number "file_n" in the folder
if ((isdir(files{file_n,1}) == false) && strcmp(FileName{length(FileName),1},"txt"))
%Check if the FileName is not a directory and if it and with *.txt
%-----------------------------------------------------------------
if (rindex (FileName{1,1}, "/") != 0) %Define the file name that is going to be treated
FileName = FileName{1,1}(rindex(FileName{1,1}, "/")+1:end);
else
FileName = FileName{1,1};
endif
%Define the output filenames
savePlot = strcat (FileFolder, "Results/", FileName, "_Plot.ps");
saveRes = strcat (FileFolder, "Results/", FileName, "_Res.txt");
saveSum = strcat (FileFolder, "Results/", "summary.txt");
%Retrieve the data values form the data file
%-------------------------------------------
matrix = data_coma (strcat(FileFolder,files{file_n,1}));
h = (matrix(2,1) - matrix(1,1)); %Characterize the data
N = length(matrix);
t = (0:h:(N-1)*h)';
data = abs(matrix(:,2));
%********************************************
%*** Prony Method will Guess the ***
%*** the starting point needed ***
%*** for the L-M algorithm ***
%********************************************
if (nargin == 3)
[alpha,c,rms] = expfit(deg,0,h,data);
alpha =(abs(alpha(:)))'; % Evaluate the module of the complex number
c = (abs(c(:)))'; % Evaluate the module of the complex number
pin = [c;alpha];
pin = vec(pin); % Rearrenge the results of the Prony
pin = [offset; pin]; % method to feed the L-M algorithm
if (param==1) %If verbose send to stdout
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% Print the informations to the SdtOut
printf("\t Prony fitted parameters: Florescence fit %.2f Volt\n",c);
printf("\t Prony fitted parameters: Life time %.2f us\n",(1./alpha)*1e6);
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
endif
endif
%--------------------------------------------------------
%_________________________________________________________
%************************************
%*** Levenberg-Marquart ***
%*** itterative - algorithm ***
%************************************
switch (deg) %Define the number of terms needed for the mathematical
case {1} %model
dp = [ds; ds; ds]; %------------------------------------------------------
F = @exprf1;
case {2}
dp = [ds; ds; ds; ds; ds];
F = @exprf2;
case {3}
dp = [ds; ds; ds; ds; ds; ds; ds];
F = @exprf3;
case {4}
dp = [ds; ds; ds; ds; ds; ds; ...
ds; ds; ds];
F = @exprf4;
endswitch
%--------------------------------------------------------
wt1 = ones(length(data),1)./sqrt(var(data)); %The moust satisfaying scheme to have a weighted fit,
%wt1 = ones(length(data),1); %this will give more importance the the first value than
%wt1 = ones(length(data),1)./sqrt(data); %the last one, normally subjected higher noise
dFdp = @dfdp; %Calculate the estimated first derivative of the mathematical model
%------------------------------------------------------------------
[f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...
leasqr (t, data, pin, F, stol, niter, wt1, dp, dFdp);%, options); % Levenberg-Marquart algorithm
%_________________________________________________________ %------------------------------
% Send to stdout is matrix singular
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
file_warn = cellstr(split(lastwarn,",")); %Check for sigular matrix
%------------------------
if (strcmp(file_warn{1,1},"warning: inverse: matrix singular to machine precision"))
printf("!!! ERROR for %s, The results can not be interpreted, the fitting exponential order my not be adequate !!!\n",FileName);
endif
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
%***************************
%** Auto correlation plot **
%***************************
resid = (f1-data);
autocRes = autocor(resid);
ResMean=mean (resid);
ResSdt=sqrt(var(resid));
cl95 = 1.96/sqrt(N); %Confidential Limit P(95%)
cl99 = 2.576/sqrt(N); %Confidential Limit P(99%)
stdev=sqrt(diag(covp1)); %Standard Deviation of the parameters estimated
stdev=stdev*2.91;
p2 = ones (length(p1)/2,2);
stdev1= ones (length(p1)/2,2);
%--------------------------------------------------------
%*******************************************
%** Rearranging the solution vector **
%** to ease the handling of the results **
%*******************************************
offset = p1(1); % From a coloumn A1;B1;A2;B2; vector...to a more confortable A1,B1;A2,B2;
for i = 1:(length(p1)-1) % ----------------------------------------------------------------------
if (mod(i+1,2))
p2(int8(i/2),2)=p1(i+1,:);
stdev1(int8(i/2),2)=stdev(i+1,:);
pNL(int8(i/2),2)=pin(i+1,:);
else
p2(int8(i/2),1)=p1(i+1,:);
stdev1(int8(i/2),1)=stdev(i+1,:);
pNL(int8(i/2),1)=pin(i+1,:);
endif
endfor
[s,i] = sort (p2(:,2)); %Sort the results form the longest lifetime to the shortest one
p2 = (p2(i,:)); %--------------------------------------------------------------
% --------------------------------------------------------------
%***********************************
%** Calculate the rute chosed for **
%** the desexcitation **
%***********************************
rute = ((p2(:,1)'*diag(1./p2(:,2))))./((p2(:,1)'*(1./p2(:,2))));
%How many molecules choose the desexcitate through rute_x
% --------------------------------------------------------------
%*************************************
%** Write the results to differents **
%** output: stdout, files, plot **
%*************************************
if (param==1) %If verbose send to stdout
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% Print the informations to the SdtOut
printf("\n\t Lebenberg-Marquart fit\n");
printf("\t Offset %.2f Volt\n",offset);
printf("\t Florescence fit %.2f Volt\n",p2(:,1));
printf("\t Life time %.2f us\n",(1./p2(:,2))*1e6);
printf("\t Imortanza %.2f\n",rute(:));
printf("\t Number of iteration needed: %i; \n\t Max iteration %i \n", iter1, niter);
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
endif
% --------------------------------------------------------------
%*******************
%** Generate Plot **
%*******************
t=t*1000000; %Change the time scale for sec to us to have more readable graphs
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% Enable to show the graphs,or disable it
fig = gcf ();
if (param == 0 || param == 2)
set (fig, "visible", "off");
else
set (fig, "visible", "on");
endif
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
subplot (221); %First subplot to rappresent the mathematical fit to the data (log scale)
%-----------------------------------------------------------------------
semilogy (t, data, "-1; Decay data;", t, f1,"-2 ;Decay exponential Model;");
title ("Experimental data fitting",'FontSize', 13);
%axis([0,t(length(t)),0.2, 1])
xlabel ("Time, usec");
ylabel ("Log luminescece decay, a.u.");
%----------------------------------------------------
%subplot(222); %Second subplot to rappresent the mathematical fit to the data (normal scale)
%---------------------------------------------------------------------------
%plot(t, data, "-1; Decay data;", t, f1,"-2 ;Decay exponential Model;");
%axis([0,t(length(t)),0.2, 1])
%xlabel ("Time, usec");
%ylabel ("Luminescece decay, a.u.");
%----------------------------------------------------
subplot(223); %Third subplot to rappresent the residials of the fit
x=[0,t(length(t))]; %----------------------------------------------------
zero=[0,0];
plot (t,resid,"-3",x,zero,"-0 ;;");
title ("Residuals",'FontSize', 13);
xlabel ("Time, usec");
ylabel ("Fit residuals, a.u.");
axis ([0,t(length(t)),-.08,.08]);
%----------------------------------------------------
subplot(224); %Forth subplot to rappresent the autocorrelation of the residials
lag = 0:1:length(t)-1; %----------------------------------------------------------------
x=[0,lag(length(lag))+1];
l=[cl95,cl95]; %Setup the confidential bans limits
l1=-l;
b=[cl99,cl99];
b1=-b;
plot(lag,autocRes,"-3;Residuals Autocorrelation;",x,l,"-4;;",x,l1,"-4;;", x,b, "-5;;",x,b1,"-5;;");
title ("Residuals autocorrelation",'FontSize', 13);
xlabel ("Lags [a.u]");
axis ([0,x(2),-.2,1]); %Plot the confidential bands
text (x(2)*2/4, 0.22, "Confidential Bands:");
text (x(2)*6/7, 0.08, "99%");
text (x(2)*6/7,-0.13, "95%");
%----------------------------------------------------
%Write the numerical results of the ritting at the place of subplot 222
%----------------------------------------------------------------------
text (x(2)*1/2, 2.8, "Summary of the Results",'HorizontalAlignment', 'center',...
'FontName','Times-Roman', 'FontSize', 22, 'Interpreter', 'none'); %Write the title
text (x(2)*1/2, 2.55, strcat("Data file: ", FileName), ...
'HorizontalAlignment', 'center','FontName','Times-Roman', 'FontSize', 17, 'Interpreter', 'none');
offset_text = sprintf ("Offset: %0.2f Volt", offset);
text (x(2)*1/2, 2.30, offset_text, 'HorizontalAlignment', ...
'center','FontName','Times-Roman', 'FontSize', 17);
for i = 1:deg % Write the exponential parameters
fit_A = sprintf("%0.2f",p2(i,1));
fit_tau = sprintf("%0.2f",(1./p2(i,2))*1e6);
fit_deg = sprintf("%d",i);
pos_text = 2.30 - i/5;
text (x(2)*1/2, pos_text, strcat("A_",fit_deg,": ", fit_A, " Volt",
' \ldots {\fontsize{16}\tau}_', fit_deg, ": ", fit_tau, ' {\fontsize{12}\mu}s'),...
'HorizontalAlignment','center','FontName','Times-Roman', 'FontSize', 17);
endfor
%Check for sigular matrix
%------------------------
if (strcmp(file_warn{1,1},"warning: inverse: matrix singular to machine precision"))
text (x(2)*1/2, 1.6, strcat("! - Warning: Singular Matrix - !"),...
'HorizontalAlignment','center','FontName','Times-Roman', 'FontSize', 20);
endif
%----------------------------------------------------
print (savePlot, "-dpsc2", "-S1000,800", "-landscape"); %Save the plot to a file
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% Print the the plots to the default priter
if (param==2)
print ("-landscape");
endif
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
%**********************************
%** Write a complete text report **
%**********************************
% Write all the informations about the data fitting in the following text file
% Summary of the results of the specific data set
ResFit_id = fopen(saveRes, 'w',"ieee-le");
fprintf(ResFit_id, "The function has been fitted with an exponential function of degree \r\n\t%i \r\n", deg);
if (nargin == 2)
fprintf(ResFit_id, "\r\n The initial parameters feeded to NLLS were found by fittig the data by Prony Method:\r\n");
else
fprintf(ResFit_id, "\r\n The initial parameters has been given by the user:\r\n");
endif
fprintf(ResFit_id, "\t Florescence fit: %.2f \r\n",pNL(:,1));
fprintf(ResFit_id, "\t Life time: %.2f us \r\n",(1./pNL(:,2))*1e6);
fprintf(ResFit_id, "\r\n Number of iteration needed to reach a minima \r\n\t%i \r\n\r\n", iter1);
fprintf(ResFit_id, "\r\n The Expenential parameters found are Sum(A_i*(expB_i)) + Offset:\r\n");
fprintf(ResFit_id, "\t Offset %.2f \r\n",offset);
fprintf(ResFit_id, "\t Florescence fit: %.2f \r\n",p2(:,1));
fprintf(ResFit_id, "\t Life time: %.2f us \r\n",(1./p2(:,2))*1e6);
fprintf(ResFit_id, "\r\n The relevance of the desexitation rute 'x' is expressed by [(A_x*B_x)/(Sum(A_i*B_i)]:\r\n");
fprintf(ResFit_id, "\t Rute %.1f%%\r\n",rute(:)*100);
fprintf(ResFit_id, "\r\n Autocorrelation Plot, First Lag \r\n\t%0.4f [Confidence 95%%, %0.4f]\r\n", autocRes(2), cl95);
fprintf(ResFit_id, "\r\n Correlation coefficient R² \r\n\t%0.4f\r\n", r21);
fprintf(ResFit_id, "\r\n The proceded files and the saved ones are located here:\r\n\t%s \r\n", pwd());
fclose(ResFit_id);
% #########################################################################
% Write the fitting results in a single comma separeted summary *.txt file
% #########################################################################
if(fopen(saveSum,'r') == -1) %Check if the file has already been created
SumFit_id = fopen (saveSum, 'a', "ieee-le"); % If not -> write the file title and the columns identificator
fprintf(SumFit_id, "Results: Exponential fit\n");
for i = 1:deg % The columns identificator
fprintf(SumFit_id, "Pre-exponential coefficient %d [Volt], Lifetime %d [µs], ", i, i);
endfor
fprintf(SumFit_id, "Data file name, Fitting Offset [Volt], Degree of the exponential model \n\n");
else
SumFit_id = fopen (saveSum, 'a', "ieee-le");
endif
for i = 1:deg
fprintf(SumFit_id, "%.4f, %.4f, ", p2(i,1),(1./p2(i,2))*1e6);
endfor
fprintf(SumFit_id, "%s, %.3f, %d \r\n", FileName, offset, deg);
fclose(SumFit_id);
% #############################################################################################
endif
endfor
endif
endfunction
function y = exprf1(x,p)
y = p(1)+p(2)*exp(-p(3)*x);
endfunction
function y = exprf2(x,p)
y = p(1) + p(2)*exp(-p(3)*x) + p(4)*exp(-p(5)*x);
endfunction
function y = exprf3(x,p)
y = p(1) + p(2)*exp(-p(3)*x) + p(4)*exp(-p(5)*x) + ...
p(6)*exp(-p(7)*x);
endfunction
function y = exprf4(x,p)
y = p(1) + p(2)*exp(-p(3)*x) + p(4)*exp(-p(5)*x) + ...
p(6)*exp(-p(7)*x) + p(8)*exp(-p(9)*x);
endfunction
## Function: dacay_lumy.m
##
## usage: decay_lumi (param, indata, deg, [pin])
##
## Small script that will help to analyse exponential decay
## phenomena. This script will coordinate different type of
## approved program (Data retrieving, prony Method, partial
## differentiation, Levenberg-Marquart non-linear fitting
## algorithm.
## The results of the calculation will be displayed with plots
## and the results will be saved in different text documents in
## "Results" folder.
The function was written to ease the analyse the of exponential
luminescence decay data acquired with a home made time resolved
spectrophotometer.
You need to enter at least 3 parameter to the function.
param
Function parameter that will define the way the output will be print
- "Silent", "silent", "S", "s": Silent mode, no text and no plot,
the fitting results are just saved in the Results folder
- "Verbose", "verbose", "V", "v": Verbose mode, everything is shown
- "Print", "print", "P", "p": Printout mode, similar to the silent mode, but
the graphical results are sent to the default printer.
%%
"Print" or "S" or "v"
%%
--
indata
The experimental data, recorded with a Lecroy LT342, need to be
in a coma separated format to be analysed.
%%
Example of a coma separate data file "coma_separated.txt":
3.25e-006,-0.855649
3.35e-006,-0.775937
3.45e-006,-0.797849
3.55e-006,-0.740049
3.65e-006,-0.854062
3.75e-006,-0.761899
3.85e-006,-0.715024
3.95e-006,-0.700986
...
%%
--
deg
Exponential fitting degree, the data are fitted with a 1 to 4 terms
exponential mathematical model. By setting this parameter the
exponential model can be set. For this parameter, set the
smallest value possible to have stable residual autocorrelation
plot.
%%
1 or 2 or 3 or 4
%%
--
[pin]
Optional, first guess of the characterizing parameters of the
exponential function, needed for the Levenberg-Marquart algorithm.
The first guess will be used to find the minimum of the
residual function. The more the initial guess are close to the final
solution the more the iterative algorithm will be stable (less
sensible to possible local minima in the residual function).
The initial parameter has to be entered in the form:
[Volt_1; Lifetime_1 (in seconds); Volt_2; Lifetime_2; ...; Lifetime_2]
%%
[.1;5e-7;3e-2;1e-7]
%%
If the user do not provide any initial guess parameter the script
will use the Prony methods to estimate the fitting solution to
feed the iterative algorithm with acceptable initial solution.
Despite this possibility to find "automaticaly" the initial guess,
it's highly recommended that the user fix this parameter. It's also
very important the use the same initial guess for all the measurement
of the same series. The mathematical fit to a data set is
a very sensible process. The iterative Levenberg-Marquart
algorithm demonstrated very good stability but it is always
very wise to fix all possible parameter.
--
------------------------------------------------------------------------------
Join us December 9, 2009 for the Red Hat Virtual Experience,
a free event focused on virtualization and cloud computing.
Attend in-depth sessions from your desk. Your couch. Anywhere.
http://p.sf.net/sfu/redhat-sfdev2dev
_______________________________________________
Octave-dev mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
