DE69320811T2 - MULTI-COMPONENT ANALYSIS FROM FT-IR SPECTRES - Google Patents

Description

The present invention relates to a method for analyzing an FT-IR spectrum of an unknown gas mixture in which the individual components and their partial pressures are unknown. However, a large set of archive spectra of molecular pure gases is known, which are measured at known pressures and with the same interferometer as the unknown sample.

In FT-IR spectroscopy for multicomponent analysis, the IR spectrum of an unknown gas mixture is available, in which the individual components and their partial pressures are unknown. What is known, however, is a large set of archive spectra of molecular pure gases, measured at known pressures and with the same interferometer as the unknown sample. Using these pure spectra, the partial pressures of the pure gases contained in the mixture should be calculated, with error limits. The errors in the values obtained in this way are due to the noise present in the spectra during measurement. Since calculating the partial pressures that best explain the spectrum of the mixture is relatively simple, it will only be briefly addressed here and the focus will be on calculating their error limits. It also deals with the optimal resolution that should be chosen in order to achieve the smallest possible error limits, as well as with the application of the non-negativity constraint for the partial pressures. In addition to gas spectra, all calculations also apply to the spectra of non-interacting liquids.

The measured spectrum of the mixture to be analyzed is called s and {Kj j = 1, ..., M} is the set of archive spectra that are used to explain s. It is assumed that all spectra are linearized , ie the negative logarithm of the transmittance is always given on the y-axis (perhaps apart from a constant coefficient). This means that the absorbance scale is used here. Then, according to Beer's law, s is a linear combination of the spectra Kj, or

Now coefficients xj must be determined that explain s as well as possible. Then the partial pressures of the pure gases contained in the mixture are determined by multiplying the measured pressures of the corresponding archive spectra by their coefficients in the linear combination.

Since only certain table points ∆1, ..., N measured at regular intervals are known from all spectra, s and Kj can be treated as vectors and equation 1 can be represented in vector form as follows:

This equation can be further simplified by collecting all Kj vectors together into a matrix K, which leads to the following matrix equation:

Kx = s, (3)

where

and x is the coefficient (column) vector containing the coefficients xj of the pure spectra.

Now all spectra that can be explained exactly by the used set of M pure spectra can be expressed in the form Kx. Even if the set {Kj j = 1, ..., M} would contain all components of the mixture, there is always noise in the spectra, so that s never lies exactly in the vector space created by the vectors Kj. So one has to make do with a linear combination Kx of the pure spectra that is as close as possible to the vector s. This means that the norm of the unexplained residual vector s-Kx has to be minimized. Thus, the present task can be represented in the form of an optimization task as follows:

s - Kx = min!.(4)

The optimization parameters are then the components of the coefficient vector x.

The optimal value of the coefficient vector depends on the vector norm used. This norm is now defined by the following condition:

This is the usual sum of squares vector norm with the exception of the constant coefficient 1/N. The purpose of this coefficient is to make the vector norms as independent as possible of the number of data N (and thus the resolution). It has no influence on the solution of optimization task 4. The solution of task 4 is now obtained by setting the gradient of the residual vector equal to zero or

x s - Kx = 0 (6)

A very simple calculation shows that this condition is fulfilled by the following coefficient vector:

xopt = (KTK)⊃min;¹KTs (7)

(of course, under the condition that N ≥ M.) This solution is known in the literature as the pseudo-inverse solution of the matrix equation Kx = s (see e.g. [1]). Due to the choice of the norm - see equation 5 - it can also be called a least squares solution. The matrix KTK is invertible if and only if the vectors Kj are linearly independent. If they are not, there is no single best solution. If the error bounds are not needed, the procedure described below is suggested for applying equation 7.

Now the inner product < a b> of two vectors a and b is defined as

The coefficient 1/N again ensures that the inner product is as independent as possible of the number of data. The usual relationship v = < v v½> between the norm and the inner product also applies. It is now clear that equation 7 can be written in the following form:

xopt = A⊃min;¹y, (9)

where

Aij = < Ki Kj> (14)

and

Yj = < Kj s> . (11)

The matrix A contains the inner products of the archive spectra used with each other and can be called the inner product matrix. If the error limits of the coefficients are not of interest, all possible inner products can be calculated in advance and collected in an archive. In a specific analysis, only those inner products that correspond to the set of archive spectra to be used in this analysis must then be selected in order to build the matrix A, calculate the vector y and solve the matrix equation Ax = y. If this is done using the Gaussian elimination method (see e.g. [2, p. 179]), A does not have to be inverted and the analysis is a very fast procedure.

All noise in the spectra creates errors in the optimal coefficient vector of equation 9. Here we only deal with the noise in the spectrum s of the gas mixture. The noise in the archive spectra Kj could also be taken into account approximately, but it is not taken into account here for two reasons. Firstly, the paper would then be too long and only of approximate accuracy. Secondly, in practical applications it is useful to register the archive spectra with much higher accuracy (i.e. by co-addition of a much larger number of individual measurement series) than the mixture spectrum obtained by a quick measurement. Thus, from a practical point of view, all noise is in the mixture spectra, and the archive spectra can be considered noise-free.

So far, the archive spectra Kj have been used as the basis vectors with which all measurements were explained in the form ΣMj=1 · jKj or Kx. However, when calculating the error limits, it is advantageous to use an orthogonal basis that generates the same vector space as the original archive spectra. These orthogonal basis vectors are constructed using the Gram-Schmidt orthogonalization method (see a textbook on linear algebra, e.g. [2, p. 138]). This means that the new basis vectors (K'j) are defined by the recursive group of equations

Thus, the old basis vectors are obtained from the new vectors by

Next, a transformation matrix Q is defined as

(Thus Q is an upper triangular matrix.) Now the transformation of the basis can be expressed very compactly as

K = K'Q. (15)

The transformation of coefficient vectors between these two bases is determined by the following equation:

x = Q⊃min;¹x'. (16)

Thus, the inverse of Q is required. After some complicated calculations, one arrives at the result that this is an upper triangular matrix, given by

The elements with j> i must be calculated in the order i = M-1, ..., 1.

Determining the best coefficient vector is very simple for the basis K'. Since this basis is orthogonal, < K'i K'j> = δij and the inner product matrix A' (see equation 10) is diagonal, we have

A' = diag(K'¹ ², K'² ², ... K'M ²). (18)

Thus, the best coefficient vector xopt is given by (see equation 9)

or

The vector x'opt gives the optimal coefficients of the orthogonal basis vectors K'j. The corresponding coefficient vector x of the real pure spectra Kj is then Q⊃min;¹x'opt. It is intuitively very clear that this is equal to the optimal coefficient vector xopt in the original basis. This result can also be formally proven as follows:

x'opt = (K'TK')&supmin;¹K'Ts = [(Q&supmin;¹)TKTKQ&supmin;¹]&supmin;¹(Q&supmin;¹)TKTs = Q(KTK)&supmin;¹QT(QT)&supmin;¹KTs = Q(KTK)⁻¹KTs = Qxopt.

This means that

xopt = Q⊃min;¹x'opt. (21)

Next, the effect of measurement noise on the optimal coefficient vector of equations 9 or 21 is considered. For this purpose, s is divided into two parts as follows:

s = s0 + se, (22)

where s⁰ is the correct noise-free mixture spectrum and se contains the noise. Since the dependence between s and xopt is linear, the coefficient vector xopt can also be divided into the correct coefficient vector x⁰opt and the error vector xeopt', for which the following applies:

x&sup0;opt = (KTK)&supmin;¹KTs&sup0; = Q&supmin;¹ (K'TK')&supmin;¹K'Ts&sup0; (23a)

and

xeopt = (KTK)-¹KTSe = Q-¹ (K'TK')-¹K'Tse (24a)

or equivalent

x&sup0;opt = A&supmin;¹y&sup0; = Q&supmin;¹A'&supmin;¹y'&sup0; (23b)

and

xeopt = A-¹ye = Q-¹A'-¹y'e (24b)

Consequently, the errors of the coefficients depend linearly on the measurements. The components of the noise vector se are normally distributed with a mean of zero. Even if the noise data of a single measurement are not normally distributed for some reason, according to the central limit theorem of probability, the sum of several noise data is always normally distributed, and in practice several individual series of measurements are always coadded to obtain a spectrum. Thus, if their standard deviation is denoted as σs, the following can be written:

it ÷ N(0, &sigma;s²). (25)

(This expression shows that it follows the normal distribution with a mean of zero and a variance σs². Note that N stands for the normal distribution, while N indicates the number of data.) Since the distribution of the noise data is known, the distributions of the Errors of the coefficients can be calculated using equation 24. This is done next.

First, a well-known result of probability theory is given without proof, which concerns the normal distribution.

It is assumed that zi ÷ N(ui, &sigma;i²) is independent and ai ε R. Then

This means that a linear combination of independent normally distributed random variables is also normally distributed. By inserting this result, the components of the vector y'e, which appears in equation 24b in the expression for xeopt, are given as:

Now another result of probability theory is needed, according to which two normally distributed random variables z₁ and z₂ are independent if and only if the correlation between them is zero. The correlation is defined by

where

cov(z 1 ,z 2 ) = E[(z 1 - E(z 1 ))(z 2 - E(z 2 ))]

and E is the expectation operator. Since the expectation values E (y'je) are zeros, two components y'je and y'ke are independent random variables if and only if E (y'je y'ke) is zero. Now

for E(sies1e) = δi1E [(sie)²], and all random variables sie are equally distributed. The independence of y'je and y'ke now results from the fact that the basis K' is orthogonal. It also follows from equations 19, 26 and 27 that

and that the components of x'eopt are independent. Now result 26 can be used again to arrive at the following equation:

If the non-orthogonal basis K is used, the equality of the left side of equation 24 results

where P is the pseudoinverse matrix

P = (KTK)⁻¹KT = 1/N A⁻¹KT (31)

Now that the distributions of the coefficient errors are known, error limits for the coefficients can also be given. Since these errors are normal distributed random variables, no upper limits can be specified for them. Instead, error limits can be calculated into which the correct coefficients fall with a given probability. It is now assumed that such an error limit j is desired, for which the coefficient xopt,j belongs with a probability p in the interval [x⁰opt,j - j, x⁰opt,j + j]. This means that

P(xeopt,j&epsi;[- j, j]) = p.

If the distribution function of the standard normal distribution (N(0,1)) is denoted as Φ, one obtains

which results in

or

The 50% error limit for the coefficient of the archive spectrum j is given, for example, by the following expression:

The variance of xeopt,j is given as follows from result 29 or 30:

It might at first seem that the original basis K is the most sensible choice to use in the analysis, since the coefficients and their variances in it are obtained from the very simple results 9 and 30. However, every time a new archive spectrum is introduced into the analysis, the inner product matrix A changes and its inverse, which is required at least in equation 31, must be recalculated, and there is no simple formula like equation 17 for inverting a general matrix. In addition, the matrix product A⊃min;¹KT must be recalculated to obtain P. The previous value of xopt also becomes useless. Using the orthogonal basis K', the coefficients and their variances are given by equations 20, 21 and 29. The following will discuss in more detail what happens when a new archive spectrum has to be added to the set of archive spectra used in the analysis when using the basis K'.

First, a new orthogonal spectrum K'M+1 must be calculated from equation group 12. When calculating the coefficients appearing in the expression K'M+1 in equation group 12, one also obtains the elements of column M + 1, which must be added to the transformation matrix Q. (The new row elements QM+1,j, j > M+1 are zeros according to equation 14.) This new column means that the inverse transformation matrix Q⊃min;¹ is also changed. However, it can be seen from equation 17 that any element of Q⊃min;¹ can only be influenced by the previously calculated elements of Q⊃min;¹ in the same column and by the elements of Q in the same and in previous columns. Thus, the previous elements of Q⊃min;¹ do not change and only M new elements Q⊃min;¹M,M+1, ..., Q⊃min;¹1,M+1 need to be calculated, which are given by equation 17 as

(Since Q⊃min;¹ is an upper triangular matrix, all elements of the new row Q⊃min;¹M+1 are either zero or one.) If the same measurement s is still analyzed as before the change of Q⊃min;¹, the following corrections must be made to the previous calculations:

1) A new component x'opt,M+1 must be added to the vector x'opt, which is obtained from equation 20.

2) According to equation 21: xopt,M+1 = Q⊃min;¹M+1x'opt = x'opt,M+1. The expression Q⊃min;¹jM+1x'opt,M+1 must be added to the other components xopt,j.

3) According to equation 29, Var(xeopt,M+1) = σs²/(N K'M+1 ²) applies. For the remaining variances Var(xopt,j) = Var(xeopt,j) the expression

be added.

Since the squares of the norms K'j ² are used repeatedly, it is useful to store them in a vector after they have been calculated in equation group 12. The use of double-precision real numbers is highly recommended in computer programs.

In the following, the invention is explained in more detail using examples and with reference to the accompanying drawings.

Fig. 1 shows a mixture spectrum s, the linear combination Kxopt of the archive spectra that best explains s, and the residual spectrum s-Kxopt that still remains unexplained (shown enlarged in the figure).

Fig. 2 shows an illustration of what happens when the archive used has errors; the same mixture spectrum was analyzed as in Fig. 1, but here the spectrum of 2-butanone was removed from the set of archive spectra used.

Fig. 3 shows a situation where the interferometer is correctly aligned, but the radiation source is circular rather than ideally point-shaped and the signal in each individual wavenumber α0 is evenly distributed over the interval from α0(1-αmax²/2) to α0.

Fig. 4 shows a situation where the observed lineshape e0 is always the convolution of the actual lineshape e, the boxcar broadening function wΩ described in Fig. 3, and the Fourier transform wL of the boxcar interferogram truncation function.

Fig. 5 shows a situation where for each component xj of the optimal non-negative coefficient vector xopt one of the following two conditions applies:

The preferred embodiment of the present invention is explained below. The method is based on an analysis of the spectrum of a simple gas mixture. The coefficients are calculated using equations 20 and 21 and their error limits using equations 29 and 32. 50% error limits are used here as they give a very good idea of the magnitude of the errors. Fig. 1 shows the mixture spectrum s to be analyzed using a set of 12 archive spectra. Table 1 shows the components of xopt with error limits. The first three components are the background spectra. Dias are all spectra or functions that are also present in the background measurement or that can arise from error data in the interferogram. They can have negative as well as positive coefficients. The mixture spectrum s and the archive spectra Kj are measured with identical devices, but artificial noise has been added to s to simulate a lower accuracy, which is assumed here. In addition, the linear combination Kxopt of the archive spectra that best explains s is shown. The residual spectrum s- Kxopt is also shown, which then remains unexplained. It can be seen that the residual spectrum consists of pure white noise. This indicates that the analysis was successful and that the 12 archive spectra used are sufficient to explain the measurement. The result of the analysis, or the optimal coefficient vector xopt, is given in Table 1. Optimal coefficient vector

The partial pressure of the archive gas j is now obtained by multiplying its measuring pressure by its coefficient xopt,j. The components that do not exist in the mixture have small positive or negative coefficients of the same order of magnitude as their error limits. If some background spectra are added to the analysis, the optimal removal of the background of the spectrum also happens automatically. In fact, the first three spectra (water, carbon dioxide and the constant function, which has a value of one everywhere) in the table are background spectra. Since their exact coefficients are not of interest here, but the coefficients are only to be eliminated, they do not even have to be pure, but can contain each other. The archive spectra can also contain some carbon dioxide and water without this affecting the value of their coefficients. The only effect is that the coefficients and error limits of the background spectra are not reliable. The coefficients of the background spectra can have negative as well as positive values. For example, it can actually happen that during the measurement of the background interferogram more carbon dioxide is present than during the measurement of the sample interferogram.

Fig. 2 shows what happens when the set of archive spectra used is not sufficient to explain the mixture spectrum. Table 2 shows the results of the analysis and the optimal coefficient vector xopt. Optimal coefficient vector

The same mixture spectrum is used as in Fig. 1, but the spectrum of 2-butanone, an important component of the gas mixture, is not included in the analysis. Now the minimized residual spectrum is not pure noise, and the coefficients of the other archive spectra have also changed. It can be seen that the minimized residual spectrum s-Kxopt is no longer white noise, but has a different structure and - what is very important - has structures on the same wave numbers on which the missing archive spectrum has spectral lines. In this way, it is possible to deduce from the residual spectrum which type of spectra should still be included in the analysis. Since the proportion of the missing spectrum must be explained as well as possible using the remaining spectra, their coefficients are also distorted and their error limits are no longer reliable. Consequently, new archive spectra must be included in the analysis until there is no longer any structure in the remaining spectrum.

If different pressure broadenings are present in s and the archive, it may be helpful to use multiple archive spectra for a single gas mixture. A better approach, however, would be to reduce the resolution so that all lines become sinusoidal (see below).

In Fig. 4, the function wL is a sine function, and the interval between its two consecutive zeros is equal to the data interval Δ in the spectral domain, which in turn is equal to 1/(2NΔx). The full width at half height (FWHH) of wL is approximately the same. The width W of wΩ is directly proportional to θ (location of the spectral line) and the area of the radiation source. Its height H is inversely proportional to θ and directly proportional to the surface brightness of the source. In an optimal situation, the full widths at half height (FWHH) of the three curves on the right-hand side of the equation are approximately equal. The total area of e0 is the product of the areas of e, wΩ and wL. The area of wL is always 1 and can be considered dimensionless.

After determining the wave number range to be dealt with, the sampling interval in the interferometer is also determined according to the Nyquist sampling theorem. However, the length of the recorded interferogram or the amplitude of the mirror movement must still be determined. which in turn is determined by the number of data used N. N denotes the number of data in the one-sided interferogram. The corresponding number used in the fast Fourier transformation algorithm is then 2N. In the following, we will examine how the number of data should be determined in order to minimize the error limits. As can be seen from equation 32, the error limits are directly proportional to the standard deviations of the coefficients, ie to the square roots of their variances. From equation 29 or 33 we get:

Thus, the error bounds are directly proportional to the standard deviation σs of the spectral noise and inversely proportional to the square root of the number of data N. According to the definitions of inner product and norm by equations 8 and 5, the remaining term in the square root is not an explicit function of N. It does, however, depend on the shapes of the archive spectra.

In the following we will examine what happens when the number of data N is reduced by a factor of 1/k. The negative effect that the coefficient 1/ N is increased by a factor of k1/2 is immediately apparent. However, the standard deviation σs of the noise of the mixture spectrum has also changed. This change is determined by Parseval's theorem

where ni and ns are the noise functions in the interferogram and spectrum, respectively. (These two random processes are a Fourier transform pair.) Since the noise ni is completely white noise, its "amplitude" is the same everywhere. If the length of the first integral is shortened to 1/k of its original value, the value of the integral is reduced by the same factor k-1. Consequently, the other integral must change by the same factor. Since the wave number range under investigation is not changed, the only possibility is that the "amplitude" of the noise ns or its standard deviation σs is reduced by the factor k1/2. This effect completely eliminates the 1/N dependence of equation 35.

As stated above, the coefficient σs/ N in equation 35 remains constant as the resolution is reduced. Thus, the only possible starting point for changes in the error limits is the expression

As already mentioned, the definitions of inner product and norm mean that this expression depends only on the shapes of the archive spectra when M is fixed. The number of data itself is not important. Any linear changes, where all archive spectra are multiplied by a constant coefficient C, change this expression in the square root by the constant C⁻¹. Now, in practice, the spectra are always calculated from the corresponding interferograms by applying the fast Fourier transform (FFT) algorithm. A fundamental property of this algorithm is that the data interval in the spectrum is 1/(2NΔx), where Δx is the sample interval of the interferogram. If the number of data is reduced by a factor of 1/k, the data interval of the spectral range will consequently increase by a factor of k. As long as the data interval (resolution/1.21) remains smaller than the full width at half height (FWHH) of the spectral lines, at least one data exists at each line and the shape of the spectral lines does not vary significantly. This means that using a resolution higher than the width of the spectral lines makes little sense. In the range of the interferogram, this means that the interferogram can be truncated without any problem as long as a significant portion of the signal is not clipped.

In the following, the function for truncating the interferogram is defined as a boxcar function with the value 1 between x = NΔx and otherwise the value 0. Since only a finite interval of the interferogram can be registered, the actual, infinitely long interferogram is always multiplied by this function. In the spectral range, this means that the spectra are convolved with the Fourier transformation wL of the truncating function, or

e⁻ = e⁻*wL (37)

where

wL( ) = 2NΔx sin(2NΔxπ ), (38)

e⁹ is the mixing spectrum or an archive spectrum and e− is the spectrum that would be obtained by reshaping the entire interferogram. The full width at half height (FWHH) of this sine function is approximately 1.21/(2NΔx), and this is the quantity used here as the resolution As long as N remains greater than (2Δx FWHH of spectral lines)-1, wL has a smaller width than the lines of e−1 and does not significantly affect their shapes. If the resolution is further reduced after this point, the spectral lines suddenly become broader and their shapes are determined primarily by wL rather than their true shapes. This means that the convolution of equation 37 then changes the spectra in a non-linear way, so that its effect is not a simple multiplication of the square root expression by a constant coefficient. If no apodization is applied, the lines start to become sinusoidal. (If apodization is applied, the noise data of the interferogram are no longer identically distributed and the error analysis is not valid.) Due to the broadening, the lines start to overlap, making it more difficult to distinguish them from each other. This in turn means that the sum term in equation 35 starts to grow. However, the growth rate depends on the number M of archive spectra used in the analysis and on the proximity of the lines to each other. For example, if the lines were originally grouped in sets of overlapping lines, the growth rate would not be as fast as if the lines were originally arranged at approximately equal intervals. Nevertheless, some rough results can be given. For example, if a set of at most 50 archive spectra is used, the growth coefficient of the square root term is usually between k1/3 and k1/2, depending on the arrangement of the lines. The square root term also depends on the number of lines in the spectra. This dependence follows the rough rule that the value of the square root is approximately inversely proportional to the square root of the average number of lines in a spectrum. Thus, it can be considered as a constant coefficient independent of N.

According to the above, the best choice for resolution would be the full width at half height (FWHH) of the spectral lines. Consequently, the recorded interferogram should extend from -1/(2 · FWHH) to 1/(2 · FWHH). However, this is only true if it is not possible to change the settings of the interferometer. If all parameters of the device can be freely determined, reducing the resolution offers two further advantages, which are discussed in more detail below.

As is well known, a non-point radiation source leads to the broadening of all spectral lines. The (opening of the) radiation source is usually round, and in this case every monochromatic spectral line is spread out into a boxcar line shape, as shown in Fig. 3. The width of the box is then directly proportional to the area of the radiation source. This means that the spectrum e− in Equation 37 is actually the true spectrum e convolved with the boxcar function wΩ arising from the non-zero area of the light source. Thus, Equation 37 can be rewritten as

e⁻ = e*wΩ*wL. (39)

Since the width of wΩ depends on the wavenumber 0 of the spectral line under consideration, a precise consideration would require the use of different wΩ for each line. Equation 39 is presented in Fig. 4. Since the full width at half height (FWHH) of the folded line is approximately equal to the sum of the FWHHs of the components of the convolution, the distortions wΩ and wL have a significant effect only when their widths are larger than that of the natural widths of the spectral lines. Thus, the signal can be easily increased by increasing the radius of the radiation source until the (maximum) width of the boxcar distortion equals the FWHH of the undistorted spectral lines. Similarly, the amount of computation can be reduced by reducing the resolution until the FWHH of the sine distortion equals that of the spectral lines. (This represents the optimal truncation of the interferogram.) In the case of a gaseous sample, however, the natural width of the lines may be so small that this situation cannot be achieved. In any case, it is still reasonable that the distortions wΩ and wL are of the same order of magnitude. Thus, this situation can be considered as a starting point. Now the number of data is reduced by a factor of k⊃min;¹. As already mentioned, this reduction leads to a broadening of the lines by a broadening of wL by a factor of k, which increases the square root expression of equation 35 by a factor of at most k1/2 if the number of archive spectra used is a few tens. On the other hand, if M is of the order of a few hundred, this factor may be of the order of k. Now, however, the area of the radiation source can be increased by a factor of k without causing a significant further increase in the widths of the lines. Since the area under the spectral lines must then increase by a factor of k due to the increase in the signal, the only possibility is that the heights of the lines are increased by the same factor k and the change in the spectra is approximately linear. Multiplying the spectra by the constant coefficient k reduces the square root expression by the coefficient k-1. This is more than sufficient to compensate for the growth of the square root expression in the non-linear interferogram. To prevent blunting surgery, however, in practice difficulties may arise when directing the magnified image of the radiation source onto the detector.

The second additional advantage offered by reducing the resolution is that now k interferograms can be recorded in the same time span as a single one was previously. Since the Fourier transform is a linear operation, co-adding these interferograms means that the corresponding spectra are also co-added. The error-free spectra e remain the same in each measurement, i.e. they are multiplied by k in the summation. This means a simple linear change in the spectra, which in turn means that the expression in the square root is multiplied by k-1. The noise of s, on the other hand, is different each time and is not summed linearly. From results 25 and 26 it can be seen that the summed noise has the distribution (N(0,kσs2)). Thus, the standard deviation σs of the noise is increased by a factor of k1/2. The overall effect is that the error limits are multiplied by a factor of k-1/2.

Finally, if all the above-mentioned effects are considered together, it becomes clear that the smallest possible resolution should be used if all parameters of the interferometer can be freely adjusted. However, the number of data N should be at least two or three times as large as the (maximum) number of archive spectra so that the structure of the residual spectrum s-Kx can be successfully investigated.

If the archive spectra were measured at a different resolution than the one at which the mixture spectrum is measured, the analysis may fail and large negative coefficients appear. A similar situation can arise when the line shapes in the archive spectra and in the mixture spectra are different due to non-linearities or different pressure broadenings. In this case, some improvement can be achieved by calculating the best non-negative solution instead of the best solution. A non-negative solution is understood to be a coefficient vector x which solves equation 4 under the condition that every component of x must be non-negative. This process introduces more information into the analysis, since a priori knowledge of the coefficients is applied. In the following, an algorithm for solving equation 4 under the condition of non-negativity constraint is derived.

The residual norm

d(x) = s-Kx

is denoted as d. Since a norm is always a non-negative set, the norm d has exactly the same minima as its square d ², so that its square can be minimized instead of the norm. Now d(x) ² is a convex function of x. This means that for every x₁, x₂, and λ, 0 < λ < 1, the following applies:

d(λx1 + (1 - λ)x2 ) 2 ? ? d(x 1 ) 2 + (1 - λ) d(x 2 ) 2 .

This is evident by applying the triangle inequality and by the fact that the geometric mean is always less than or equal to the arithmetic mean. The convexity implies that the square norm has only one minimum point, which makes minimization considerably easier.

In particular, if other coefficients are kept constant and only one coefficient xj is varied, d ² is a convex function of one variable. Thus, at the optimal point xopt there are two possibilities, either

(where all components except the j'th are bound to xopt) or, if the zero point of the derivative is not in the allowed range xj ≥ 0,

xopt,j = 0,

which means that xopt,j lies on the border between the allowed and forbidden regions. This can be proven as follows:

1) If the zero point of &part; d(xopt) ²/ &part;xj lies in the allowed range of x ≥ 0, then xopt,j must be uniquely equal to this zero point.

2) If the origin of &part; d(xopt) ²/ &part;xj lies in the forbidden region of xj < 0, the derivative due to the convexity of d(x) ² is positive if xj ≥ 0. Thus, if xopt,j > 0, decreasing xopt,j would decrease the value of d ² without leaving the allowed region. Thus, the only possibility is that xopt,j = 0.

Now the following optimal conditions can be established:

For the only solution point xopt of the minimization task 4 with non-negative components of x, one of the following conditions applies for each component xj:

or

This is demonstrated in Fig. 5. The so-called Kuhn-Tucker criteria (see a textbook on optimization, e.g. [3]), which are a necessary and sufficient condition for optimality in the case of a convex object function, would lead to pretty much the same condition.

The partial derivatives are very easy to calculate, and you get

x( d(x) ²) = 2 (Ax - y), (41)

where A is the inner product matrix defined in equation 10 and y is given by equation 11. The individual partial derivatives are the components of this gradient, or

If some background spectra are included in the set of archive spectra used, it is sensible not to subject their coefficients to the non-negativity constraint. The background spectra include all components that are also present during the measurement of the background (such as water and carbon dioxide), as well as a constant function and simple sine waves that can be generated by erroneous data in the interferogram. If this is also taken into account, the minimum point can be determined, for example, using the following algorithm:

1) A starting point x = x�0 is chosen, e.g. x�0 = (0,0, ..., 0)T or x�0 = A⊃min;¹y. j = 1 is inserted.

2) The following is calculated:

which according to equation 42 is the zero point of &part; d(x) ²/ &part;xj . If xj < 0 and if the coefficient xj is subject to the non-negativity constraint, xj = 0 is reinstated. j = j + 1 is inserted.

3) If j ≤ M (where M is the number of archive spectra used in the analysis), return to step 2. Otherwise, go to step 4).

4) The following is calculated:

G = 1/2 x( d(x) ²) = Ax - y.

If for each component Gj of G either Gj < ₁ or (Gj > 0 and xj < ₂), where ₁ and ₂ are suitable small real numbers, stop. Otherwise, insert j = 1 and return to step 2.

Reference works

1. C. L. Lawson and R. J. Hanson, Solving least squares problems (Prentice-Hall, New Jersey, 1974), p. 36.

2. B. Noble and J. W. Daniel, Applied linear algebra (Prentice-Hall, New Jersey, 1977), 2nd edition.

3. M. S. Bazaraa and C. M. Shetty, Nonlinear programming - theory and algorithms (John Wiley & Sons, New York, 1979), p. 137.

Claims (4)

1. A method for analyzing multicomponent FT-IR spectra of unknown gas mixtures by an analyzer, in which a set of archive FT-IR spectra of molecular pure gases measured at known pressures are stored as vectors k with N components, and in which inner products for all possible combinations of two vectors k are calculated and stored by the analyzer, comprising the following method steps: - Selecting M archive vectors k for the analysis of the unknown gas mixture, and thereby forming a basis matrix K, - Forming an orthogonal basis matrix K' by deforming the vectors k chosen for the matrix K into orthogonal vectors k', so that inner products of these vectors k' with each other but not with themselves result in zero, and calculating and storing the inner products of these orthogonal vectors k' with themselves, - Form and store a basis transformation matrix Q, such that K = K'Q, - Forming and storing an inverse matrix Q⊃min;¹ of the matrix Q, - Measuring an FT-IR spectrum s for the unknown gas mixture as a vector with N components, - Calculate an optimal coefficient vector x' for the spectrum s in the orthogonal basis by using the equation for its components x'j, - Calculate a coefficient vector x such that x = Q⊃min;¹ x', and - Calculating the partial pressures of the pure gases in the gas mixture to be analyzed by multiplying the measuring pressures of the corresponding archive spectra by their coefficients xj, and thus forming an analysis result for the spectrum s. 2. The method of claim 1, further comprising: - Compare the measured spectrum s with its analysis result and in case of a mismatch adding a new archive vector kM+1 to the base matrix K by - Calculate a new orthogonal vector k'M+1, - Calculate and store an inner product of the new orthogonal vector k'M+1 with itself, - Add a new column to the transformation matrix Q and store this new matrix, - Forming and storing a new inverse matrix Q⊃min;¹ based on the new matrix Q by adding a new column to the inverse matrix Q⊃min;¹, - Calculate a new component x'M+1 for the orthogonal coefficient vector x' of the spectrum s to be analyzed by the equation and adding them to the vector x', and - Calculate a new coefficient vector x such that x = Q⊃min;¹ x'. 3. The method of claim 1 or 2, further comprising: - Evaluation of a deviation &sigma;s for the noise of the spectrum s, - Calculate variances for the coefficients xj of the vector x by using the equation and - Calculate the error limits vj (p) for the coefficients xj by the equation where p is the probability that the error for the archive spectrum j lies between the limits (-vj, vj) and Φ⊃min;¹ is the inverse of the distribution function of the one-dimensional standardized normal distribution. 4. Method according to claims 2 and 3, further comprising: - Calculate the variances for the new coefficient xM+1 using the equation and - Calculate the values to be added to the variances of the coefficients xj, j = 1, ..., M by the equation - Add these values to the variances of the coefficients xj, j = 1, ..., and - Calculate new error limits vj (p) for the coefficients xj by the equation where p is the probability that the error for the archive spectrum j lies between the limits (-vj, vj) and Φ⊃min;¹ is the inverse of the distribution function of the one-dimensional standardized normal distribution.