This code determines the limit of detection from experimental data in a linear calibration curve. See the following article for details:

Hans-Peter Loock and Peter D. Wentzell, Detection limits of chemical sensors: Applications and misapplications, Sensors and Actuators B: Chemical, 173 (2012) 157–163, https://doi.org/10.1016/j.snb.2012.06.071

% Created by: Arthur Giron Santos 

 % Date: June, 2022

% Created by: Arthur Giron Santos 

 % Date: June, 2022

 % Version: 1.0

 % Objective: to calculate the Limit of Detection of chemical sensors using the following data provided:

 % x (concentration), y (signal), t (Student t-function) and k (number of replicate measurements).

 % cite as:

 % A.G. Santos, and H-P Loock (2022) Limit of Detection source code (Version 1.0) [Source code] <link_website>

 % based on: Hans-Peter Loock and Peter D. Wentzell, Detection limits of chemical sensors: Applications and misapplications,

 % Sensors and Actuators B: Chemical, 173 (2012) 157–163, https://doi.org/10.1016/j.snb.2012.06.071

 % In the code below we refer to equation numbers in Sens. Act. B: Chem, 173, (2012) 157–163.

 % Sample data => replace those row vectors with your data

     x = [1,4,9,16,25,36,49,64,81,100];

     % y = [10.967044171877, 13.7266562041864, 21.6774325819587, 28.0501802939311, 36.0314926863133, 51.9667357806104, 55.1532263138279, 74.6573287557684, 71.9956136960292, 129.299736817629];

     y= [0.948765056,2.868444736,9.156115413,10.93922882,29.1421891,35.21056678,47.19701892,78.89101104,71.86086959,88.32222164];

     t = 3; % Value of Student T-function

     k = 1; % number of replicates for each value

     [LOD_1 LOD_2 LOD_3] = calculateLOD(x,y,t,k);

     output = "LOD_1 = "+LOD_1+ " (dashed blue line; see equation (11)"

     output = "LOD_2 = "+LOD_2+ " (solid black line; Currie and Svehla (1994); see equation (15)"

     output = "LOD_3 = "+LOD_3+ " (dashed red line; should not be used)"

function [LOD_1 LOD_2 LOD_3] = calculateLOD(x,y,t,k)

     % --------------------------------------------------------------------------

     % Calculate LOD

     % --------------------------------------------------------------------------

     % Assuming that X and Y are ROW ARRAYS

         n = length(x);

     % perform linear regression

         x_aux = ones(size(x',1),1); % auxiliary variable

         x_aux = [x_aux x'];

         % [c,~,~,~,stats] = regress(y',x_aux); % calculates: a (slope), b (intersept), R^2 (correlation coefficient)

         c = polyfit(x,y,1);

         a = c(1); % slope (=sensitivity)

         b = c(2);

         % R = stats(1);

     % Calculate determinant in the denominator of equations (4) and (5)

         sum_x = sum(x); % auxiliary variable

         sum_x2 = sum(x.^2); % auxiliary variable

         D = sum_x2*n - sum_x*sum_x; % equation (6)

         d = zeros(1,length(x));

         for i = 1:length(x) % numerator in (8)

             d(i) = (y(i) - (a*x(i) + b))^2;

         end

         sigma_y = sqrt( sum(d)/(n-2) ); % equation (8)

         aux = (n/D - ((abs(a)/t)/sigma_y)^2);

     % given for completeness, but should not be used:

         LOD_3 = t*sigma_y / a;

     % LOD_2 from equation (11):

         LOD_2 = (2/D)/aux * ( sum_x - sqrt( sum_x^2 - D^2 * aux * (sum_x2/D + 1/k) ) );

     % LOD_3 determined from equation (15), which is derived from

     % Currie and Svehla, Nomenclature for the presentation of results of chemicalanalysis, Pure and Applied Chemistry 66 (1994) 595–608

         LOD_1 = 2*t*sigma_y / ( n * (t*sigma_y)^2 - D*a^2 ) * ( t*sigma_y*sum_x - sqrt(D*a^2 * sum_x2 + (a*D)^2) );

     % --------------------------------------------------------------------------

     % PLOT GRAPHIC

     % --------------------------------------------------------------------------

         x_fit = zeros(1,length(x));

         y_fit = zeros(1,length(x));

         sigma_x = zeros(1,length(x)); % for ERROR BARS in x

         MpC_x = zeros(1,length(x)); % mean_x + conf_x

         MmC_x = zeros(1,length(x)); % mean_x - conf_x

         MpS_x = zeros(1,length(x)); % mean_x + sigma_x

         MmS_x = zeros(1,length(x)); % mean_x - sigma_x

         y_fit = [y_fit, 0]; % "y_fit" needs an extra element for the plot

         y_fit(1) = b;

         for i = 1:length(x)

             x_fit(i) = (y(i) - b)/a; % x-values of the fit line

             sigma_x(i) = sigma_y/abs(a) * sqrt( 1/k + ( x_fit(i)^2 * n + sum_x2 - 2*x_fit(i)*sum_x )/D ); % equation (9)

             MpC_x(i) = x_fit(i) + t*sigma_x(i);

             MmC_x(i) = x_fit(i) - t*sigma_x(i);

             MpS_x(i) = x_fit(i) + sigma_x(i);

             MmS_x(i) = x_fit(i) - sigma_x(i);

             y_fit(i+1) = a*x(i) + b;

         end

     %auxiliary variables

         x0 = [0, x];

         aux = b+a*LOD_2/2;

         LOD3x = [0, LOD_3, LOD_3];

         LOD3y = [aux, aux, 0];

         LOD2x = [0, LOD_2, LOD_2];

         LOD2y = [aux, aux, 0];

         aux = b+a*LOD_1/2;

         LOD1x = [0, LOD_1, LOD_1];

         LOD1y = [aux, aux, 0];

     %plots

         f = figure;

         f.Position = [10 10 960 540];

         plot(x,y,"ko")

         ylim([0 inf])

         ax = gca;

             ax.YAxisLocation = 'origin';

         title('Limit of Detection')

         xlabel("Concentration")

         ylabel("Signal")

         txt={"LOD_1 = "+LOD_1+" (dashed blue line; eq'n [11]) ", "LOD_2 = "+LOD_2+ " (solid black line; eq'n [15])","LOD_3 = "+LOD_3+ " (dashed red line; don't use)"};

         dim = [.6 .2 .1 .1];

         annotation("textbox",dim,'String',txt,'FitBoxToText','on','BackgroundColor','white')

     hold on

         %errorbars

             yerr = zeros(1,length(x)) + sigma_y;

             xerr = sigma_x;

             errorbar(x,y,yerr,yerr,xerr,xerr,'.')

         %red line - linear fit

             plot(x0,y_fit,"r-")

         %blue lines - confidence intervals

             plot(MmC_x,y,"b-",MpC_x,y,"b-")

         %LOD points

             plot(LOD3x,LOD3y,"r--",LOD2x,LOD2y,"k-",LOD1x,LOD1y,"b--")

     hold off

 end