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.
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. e = 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 e.
A = ec
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.
Measurement of the spectra of known concentrations of the separate
species allows their analytical sensitivity e at each wavelength to be
determined. Then it follows that:
Aw1=ec1,w1
cc1
+ ec2,w1
cc2
+ ec3,w1
cc3
+ ... for all n
species.
Aw2=ec1,w2
cc1
+ ec2,w2
cc2
+ ec3,w2
cc3+ ...
and so on for all
wavelengths. This set of linear equations is written compactly in
matrix form:
A
= eC
where
A is the w-length vector of measured signals (i.e. the signal spectrum)
of the mixture, e is the n x w rectangular matrix of the known e-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 e-1, the matrix inverse of e, and using the property that any matrix times its inverse is unity:
C
= e-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
thee 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 (eTe)-1eT:
(eTe)-1eTA
= (eTe)-1eTeC
= (eTe)-1(eTe)C
But the quantity
(eTe)-1(eTe)
is a matrix times its inverse and is therefore unity. Thus:
C=
(eTe)-1eTA
In this expression, eTe 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 e 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=
(eTV-1e)-1
eTV-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 linearily 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 wavelengh 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.
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.xlsandRegressionDemo.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=
(eTe)-1eTA
are performed in thses 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:
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
Matlab
or Octaveare really the natural computer approachs 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 *. 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. Alternatively, 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 spreasheet 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 eaily accomplished in Matlab/Octave by adding a column of
1s to the A matrix containing the absorption spectrum of each of the
components.:
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
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=exponentionally broadened
Gaussian; 6=equal-width Gaussians; 7=Equal-width Lorentzians;
8=exponentionally 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 linearily
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. This 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 much more
quickly on wheat paste preparations, 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 a similar predictable
type must 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.Revised March, 2013. 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:
.
