index previous next

Curve fitting B: Multicomponent Spectroscopy

[Classical Least Squares (CLS) calibration]   [Inverse Least Squares (ILS) calibration]   [Spreadheets]   [Matlab or Octave

The spectroscopic analysis of mixtures, when the spectra of the components overlap considerably, can be performed using special calibration methods based on a type of linear least-squares called multiple linear regression. This method is widely used in multi-wavelength instruments such as diode-array, Fourier transform, and automated scanning spectrometers. In this case I will derive the math, because it's instructive to do so:

Definitions:

n = number of distinct chemical species in the mixture
s = number of samples
s1,s2 = sample 1, sample 2, etc.
c = molar concentration
c1, c2 = species 1, species 2, etc.
w = number of wavelengths at which signal is measured
w1, w2 = wavelength 1, wavelength 2, etc.
Ɛ = analytical sensitivity (slope of a plot of A vs c)
A = analytical signal

Assumptions:

The analytical signal, A (such as absorbance in absorption spectroscopy, fluorescence intensity in fluorescence spectroscopy, and reflectance in reflectance spectroscopy) is directly proportional to concentration, c. The proportionality constant, which is the slope of a plot of A vs c, is Ɛ.

A = Ɛc

The total signal is the sum of the signals for each component in a mixture:

Atotal = Ac1 + Ac2 + ... for all n species.


Classical Least Squares (CLS) calibration. This method is applicable to the quantitative analysis of a mixture of species when the spectra of the individual species are known and where total signal of the mixture is the sum of the signals for each component in the mixture. Measurement of the spectra of known concentrations of the separate species allows their analytical sensitivity Ɛ at each wavelength to be determined. Then it follows that:

Aw1=Ɛc1,w1 Ɛc1 + Ɛc2,w1 Ɛc2 + Ɛc3,w1 cc3 + ... for all n species.

Aw2=Ɛc1,w2 cc1 + Ɛc2,w2 cc2 + Ɛc3,w2 cc3 + ...

and so on for all wavelengths. This set of linear equations is written compactly in matrix form:

A = ƐC

where A is the w-length vector of measured signals (i.e. the signal spectrum) of the mixture, Ɛ is the n x w rectangular matrix of the known Ɛ-values for each of the n species at each of the w wavelengths, and C is the n-length vector of concentrations of all the species.

If you have a sample solution containing unknown amounts of those species, then you measure its spectrum A and seek to calculate the concentration vector C. In order to solve the above matrix equation for C, the number of wavelengths w must be equal to or greater than the number of species n. If w = n, then we have a system of n equations in n unknowns which can be solved by pre-multiplying both sides of the equation by
Ɛ-1, the matrix inverse of Ɛ, and using the property that any matrix times its inverse is unity:

C = Ɛ-1A

Because real experimental spectra are subject to random noise (e.g. photon noise and detector noise), the solution will be more precise if signals at a larger number of wavelengths are used, i.e. if w > n. But then the equation can not be solved by matrix inversion, because the Ɛ matrix is a w x n matrix and a matrix inverse exists only for square matrices. A solution can be obtained in this case by pre-multiplying both sides of the equation by the expression (ƐTƐ)-1ƐT:

(ƐTƐ)-1ƐTA = (ƐTƐ)-1ƐTƐC = (ƐTƐ)-1(ƐTƐ)C

But the quantity (ƐTƐ)-1(ƐTƐ) is a matrix times its inverse and is therefore unity. Thus:

C = (ƐTƐ)-1ƐTA

In this expression, ƐTƐ is a square matrix of order n, the number of species. In most practical applications, n is typically only 2 to 5, this is not a very big matrix to invert, no matter how many wavelengths are used. In general, the more wavelengths are used the more effectively the random noise will be averaged out (although it won’t help to use wavelengths in spectral regions where none of the components produce analytical signals). The determination of the optimum wavelength region must usually be determined empirically. 

Two extensions of the CLS method are commonly made. First, in order to account for baseline shift caused by drift, background, and light scattering, a column of 1s is added to the Ɛ matrix. This has the effect of introducing into the solution an additional component with a flat spectrum; this is referred to as “background correction”. Second, in order to account for the fact that the precision of measurement may vary with wavelength, it is common to perform a weighted least squares solution that de-emphasizes wavelength regions where precision is poor:

C = (ƐT Ɛ-1Ɛ)-1 ƐT V-1 A

where V is an w x w diagonal matrix of variances at each wavelength. In absorption spectroscopy, where the precision of measurement is poor in spectral regions where the absorbance is very high (and the light level and signal-to-noise ratio therefore low), it is common to use the transmittance T or its square T2 as weighting factors.

The method is in principle applicable to any number of overlapping components. Its accuracy is limited by how accurately the spectra of the individual components are known, the amount of noise in the signal, the extent of overlap of the spectra, and the linearity of the analytical curves of each component (the extent to which the signal amplitudes are proportional to concentration).  In practice, the method does not work well with old-style instruments with manual wavelength control, because of insufficient wavelength reproducibility, but it works better with automated computer-controlled scanning instruments and is perfectly suited to diode-array and Fourier transform instruments. The method is widely applied in absorption spectrophotometry, where the well-known deviations from analytical curve linearity set a limit to the performance to this method (and which can be avoided, in a related method, by applying curve fitting to the transmission spectra). 

Inverse Least Squares (ILS) calibration. ILS is a method that can be used to measure the concentration of an analyte in samples in which the spectrum of the analyte in the sample is not known beforehand. Whereas the classical least squares method models the signal at each wavelength as the sum of the concentrations of the analyte times the analytical sensitivity, the inverse least squares methods use the reverse approach and models the analyte concentration c in each sample as the sum of the signals A at each wavelength times calibration coefficients m that express how the concentration of that species is related to the signal at each wavelength:

cs1 = mw1As1,w1 + mw2As1,w2+ mw3As1,w3 + ... for all w wavelengths.

cs2 = mw1As2,w1 + mw2As2,w2 + mw3As2,w3 + ...

and so on for all s samples. In matrix form

C = AM

where C is the s-length vector of concentrations of the analyte in the s samples, A is the w x s matrix of measured signals at the w wavelengths in the s samples, and M is the w-length vector of calibration coefficients.

Now, suppose that you have a set of standard samples that are typical of the type of sample that you wish to be able to measure and which contain a range of analyte concentrations that span the range of concentrations expected to be found in other samples of that type. This will serve as the calibration set. You measure the spectrum of each of the samples in this calibration set and put these data into a w x s matrix of measured signals A. You then measure the analyte concentrations in each of the samples by some reliable and independent analytical method and put those data into a s-length vector of concentrations C. Together these data allow you to calculate the calibration vector M by solving the above equation. If the number of samples in the calibration set is greater than the number of wavelengths, the least-squares solution is:

M = (ATA)-1ATC

(Note that ATA is a square matrix of size w, the number of wavelengths, which must be less than s). This calibration vector can be used to compute the analyte concentrations of other samples, which are similar to but not in the calibration set, from the measured spectra of the samples:

C = AM

Clearly this will work well only if the analytical samples are similar to the calibration set. The advantage of this method is that the spectrum of an unknown sample can be measured much more quickly and cheaply than the more laborious standard reference methods that are used to measure the calibration set, but as long as the unknowns are similar enough to the calibration set, the concentrations calculated by the above equation will be accurate enough for many purposes.
Spreadheets. Most modern spreadsheets have basic matrix manipulation capabilities and can be used for multicomponent calibration, for example Excel and OpenOffice Calc.  The spreadsheets RegressionDemo.xls and RegressionDemo.ods (for Excel and Calc, respectively) demonstrate the classical least squares procedure for a simulated spectrum of a 5-component mixture measured at 100 wavelengths. A screen shot is shown below. The matrix calculations described above that solves for the concentration of the components on the unknown mixture:
C = (ƐTƐ)-1ƐTA
are performed in these spreadsheets by the TRANSPOSE (matrix transpose), MMULT (matrix multiplication), and MINVERSE (matrix inverse) array functions, laid out step-by-step in rows 123 to 158 of this spreadsheet. These array operations may be combined into one cell equation:

C = MMULT(MMULT(MINVERSE(MMULT(TRANSPOSE(E);E));TRANSPOSE(E));A)

where C is the vector of the 5 concentrations of all the components in the mixture, E is the 5 x 100 rectangular matrix of the known sensitivities (e.g. absorptivities) for each of the 5 species at each of the 100 wavelengths, and A is the vector of measured signals at each of the 100 wavelengths (i.e. the signal spectrum) of the unknown mixture. (Note: spreadsheet array functions like this must be entered by typing Ctrl-Shift-Enter, not just Enter as usual, See http://www.openofficetips.com/blog/archives/2004/10/array_formulas.html).

.
OpenOffice Calc spreadsheet demonstrating
the CLS procedure for the measurement of a 5-component unknown mixture at 100 wavelengths

Alternatively, you can skip over all the details above and use the LINEST function, in Excel or OpenOffice Calc, which performs this type of calculation in a single function statement. This is illustrated in RegressionDemo.xls, in cell P23. A significant advantage of the LINEST function is that it can automatically compute the standard errors of the coefficients and the R2 value in the same operation; using Matlab or Octave, that would require some extra work. (LINEST is also an array function that must also be entered by typing Ctrl-Shift-Enter, not just Enter). Note that this is the same LINEST function that was used for polynomial least-squares; the difference is that in polynomial least-squares, the multiple columns of x values are computed, for example by taking the powers (squares, cubes, etc) of the x's, whereas in the multicomponent CLS method the multiple columns of x values are experimental spectra of the different standard solutions. The math is the same, but the origin of the x data is very different.

Matlab or Octave are really the natural computer approaches to multicomponent analysis because they handle all types of matrix math so easily, compactly, and quickly. In these languages, the notation is very compact: the transpose of matrix A is A', the inverse of A is inv(A), and matrix multiplication is designated by an asterisk (*). Thus the solution to the classical least squares method above is written in Matlab/Octave notation as

C = inv(E'*E)*E'*A

where E is the rectangular matrix of sensitivities at each wavelength for each component and A is the observed spectrum of the mixture. Note that the Matlab/Octave notation is not only shorter than the spreadsheet notation, it's also closer to the traditional mathematical notation. Even more compactly, you can write C = A * pinv(E')' or C = A/E, which yield essentially the same results in the same execution time. 
The script RegressionDemo.m (for Matlab or Octave) demonstrates the classical least squares procedure for a simulated absorption spectrum of a 5-component mixture at 100 wavelengths, illustrated above. Most of this script is just signal generation and plotting; the actual least-squares regression is performed on one line:
MeasuredAmp = ObservedSpectrum*A'*inv(A*A')
where different symbols are used for the variables: "A" is a matrix containing the spectrum of each of the components in each of its rows and "ObservedSprectum" is the observed spectrum of the unknown mixture. In this example the dots represent the observed spectrum of the mixture (with noise) and the five colored bands represent the five components in the mixture, whose spectra are known but whose concentrations in the mixture are unknown. The black line represents the "best fit" to the observed spectrum calculated by the program. In this example the concentrations of the five components are measured to an accuracy of about 1% relative (limited by the noise in the observed spectrum).  Comparing RegressionDemo.m to its spreadsheet equivalent, RegressionDemo.ods, both running the same  computer, you can see that the Matlab/Octave code computes and plots the results quicker that the spreadsheet, although both take no more than a fraction of a second for this example.

Extensions: The extension to "background correction" is easily accomplished in Matlab/Octave by adding a column of 1s to the A matrix containing the absorption spectrum of each of the components.:

background=ones(size(ObservedSpectrum));
A=[background A1 A2 A3];


where A1, A2, A3... are the absorption spectra of the individual components.  Performing a T-weighted regression is also readily performed:

  weight=T;
  MeasuredAmp=([weight weight] .* A)\(ObservedSpectrum .* weight);


where T is the transmission spectrum.  Here, the matrix division backslash "\" is used as a short-cut to the classical least-squares matrix solution (c.f.
http://www.mathworks.com/help/techdoc/ref/mldivide.html).

The cls.m function: Ordinarily, the calibration matrix M is assembled from the experimentally measured signals (e.g. spectra) of the individual components of the mixture, but it is also possible to fit a computer-generated model of basic peak shapes (e.g. Gaussians, Lorentzians, etc) to a signal to determine if that signal can be represented as the weighted sum of overlapping basic peak shapes. The function cls.m computes such a model consisting of the sum of any number of peaks of known shape, width, and position, but of unknown height, and fit it to noisy x,y data sets. The syntax is

heights=cls(x,y,NumPeaks,PeakShape,Positions,Widths,extra)

where x and y are the vectors of measured data (e.g. x might be wavelength and y might be the absorbance at each wavelength), 'NumPeaks' is the number of peaks, 'PeakShape' is the peak shape number (1=Gaussian, 2=Lorentzian, 3=logistic, 4=Pearson, 5=
exponentially broadened Gaussian; 6=equal-width Gaussians; 7=Equal-width Lorentzians; 8=exponentially broadened equal-width Gaussian, 9=exponential pulse, 10=sigmoid, 11=Fixed-width Gaussian, 12=Fixed-width Lorentzian; 13=Gaussian/Lorentzian blend; 14=BiGaussian, 15=BiLorentzian), 'Positions' is the vector of peak positions on the x axis (one entry per peak), 'Widths' is the vector of peak widths in x units (one entry per peak), and 'extra' is the additional shape parameter required by the exponentionally broadened, Pearson, Gaussain/Lorentzian blend, BiGaussain and BiLorentzian shapes.  Cls.m returns a vector of measured peak heights for each peak. clsdemo.m is a demonstration script that creates some noisy model data, fits it with cls.m, then computes the accuracy of the measured heights.

Weighted linear regression: The custom Matlab/Octave function "tfit.m" simulates the measurement of three absorbing components by weighted linear regression, demonstrates the effect of  the amount of noise in the signal, the extent of overlap of the spectra, and the linearity of the analytical curves of each component, and also compares the results to a more advanced method that applies curve fitting to the transmission spectra rather than to the absorbance spectra.


The Inverse Least Squares (ILS) technique is demonstrated in Matlab by this script and the graph above. The math, described above, is similar to the Classical Least Squares method, and can be done by any of the Matlab/Octave or spreadsheet methods described in this section. This example is a real data set derived from the near infrared (NIR) reflectance spectroscopy of agricultural wheat samples analyzed for protein content. In this example there are 50 calibration samples measured at 6 wavelengths. These calibration samples had already been analyzed for protein content by a reliable, but laborious and time-consuming, reference method.  The purpose of this calibration is to establish whether near-infrared reflectance spectroscopy, which can be measured on wheat paste preparations much more quickly than wet chemical methods, correlates to their protein content as determined by the reference method. These results indicate that it does, at least for this set of 50 wheat samples, and therefore is it likely that near-infrared spectroscopy should do a pretty good job of estimating the protein content of similar unknown samples.  The key is that the unknown samples must be similar to the calibration samples (except for the protein content). However, this is a very common analytical situation in commerce, where large numbers of samples of products of a similar predictable type must often be analyzed quickly and cheaply.  

Note: you can right-click on any of the m-file links above and select Save Link As... to download them to your computer for use within Matlab.
index previous next
Revised September, 2014. This page is maintained by Prof. Tom O'Haver , Department of Chemistry and Biochemistry, The University of Maryland at College Park. Comments, suggestions and questions should be directed to Prof. O'Haver at toh@umd.edu.
Unique visits since May 17, 2008: