Gibbs and Pichler (30) have made a very important connection between the Fourier power spectrum and the Hadamard (again called Walsh) power spectrum. With the discovery of their important link, the two domains can be used somewhat interchangeably and thus increase the usefulness of the Hadamard transform. The familiar Shift Theorem of Fourier theory has been proved for Hadamard theory (25), fast cyclic convolution and correlation algorithm have been developed (31), and several applications have been successfully completed (30) J. E. Gibbs and F. R . Pichler, ibid.. p 51. (31) D. A . Pitassi, ibid., p 130.
(32-34). These and other developments promise to make the Hadamard transform a useful tool for chemical data analysis. Received for review March 26, 1973. Accepted May 31, 1973. Paper presented at the 1973 Pittsburgh Conference on Analytical Chemistry and Spectroscopy, Cleveland, Ohio, March 1973. Work supported in part by the National Science Foundation. GP-36578X. (32) I . A. Davidson, ibid.. p 1 7 7 . (33) S. J. Campanella and G. S. Robinson, ibid.. p 199 (34) J. W. Carl and M. Kabrisky, ibid.. p 203.
Deconvolution Method for Identification of Peaks in Digitized Spectra Geert Brouwer and J. A. J. Jansen Philips
Research Laboratories, Eindhoven, The Netherlands
An intermediate step in the automatic evaluation of complex spectra is the determination of the positions and areas of the individual lines. The influence of the background is reduced by fitting the first derivative of the known profile to the first derivative of the spectrum. This method ultimately results in one set of deconvolution coefficients for the determination of the positions of the lines and another set for obtaining the areas. I n a specific case Gaussian profiles are assumed. To a certain degree, doublets are resolved. A 30,000-channel spectrum is processed in 1 minute. The method introduces a bandwidth limitation that has an influence on the signal-tonoise ratio and resolution. A comparison is made to other methods.
Spectra of varying complexity appear in several areas of analytical chemistry and computer programs are used for line identification. In this paper, attention is restricted to spectra where the line shape is reproducible, because it is determined by a stable instrument. This is the case with prism and grating spectrographs of moderate resolution, gamma-ray spectrographs and, to a lesser extent, mass spectrographs. Moreover, the line profiles can be approximated quite well by Gaussian curves with half-widths that vary only slowly with the position of the line in the spectrum. This knowledge is exploited in our peak-searching program. The data reduction program starts with a spectrum where the quantitative information obeys the superposition principle. Thus, in a photographic spectrum, the measured transmittances must first be converted to intensities. In the absence of background, the lines in a spectrum can be identified by a least-squared fit of the standard line profile, of which both position and amplitude are varied. The introduction of background leads to erroneous positions and areas. If the assumption is made that the slope of the background is constant in the neighborhood of the line that is processed, the influence of back-
ground can be eliminated in the following way: the first derivative of the standard profile is fitted to the first derivative of the measured spectrum, as is shown in the next section. The difference between both methods can be explained by their transfer characteristics: besides the common high-frequency limitation, the latter has an additional low-frequency cut-off that rejects low-frequency noise or background. Any data reduction program introduces some bandwidth limitation in the information flow and this affects both resolution and signal-to-noise ratio. In digital computer programs, the measurements are available in sampled form, and in this paper a deconvolution technique is developed for line identification. It yields two sets of deconvolution coefficients, one for the determination of positions, the other for the areas. The linear operations: differentiation and least-squares fitting need not be performed sequentially but can be compacted into one operation ( I ) . This enhancement of speed is desirable because the data reduction is performed simultaneously with the reading of the information to avoid the necessity of a vast computer memory space. Overlapping peaks are resolved by the deconvolution method presented in this paper, as far as this is permitted by the resolution of the spectrograph and the data reduction program. The sampling of the spectrum may have an influence on the final results. Special attention is paid to the width of the plate reader slit that is used in connection with photographic spectra. A special section is devoted to a comparison of the various data reduction methods that have been published.
INSTRUMENTAL The plate reader employed in our laboratory was built by Witmcmr and Van Gool ( 2 ) . It scans photographic plates u p to 25 cm in length and measures the transmittance automatically a t a rate of 0.5 cm/sec. A built-in measuring system generates trigger puls(1) A. Savitzky. and M. J. E. Goiay, Anal. Chem.. 36,1627 (1964). (2) A . W . Witmer and G . H . Van Gool, Colloquium Spectroscopicum lnternationaie XVI, Heidelberg, preprints, 1 9 7 1 , p 254.
A N A L Y T I C A L C H E M I S T R Y , V O L . 45, NO. 13, NOVEMBER 1973
2239
The best position of the line is given by the zero of the location correlation function, /3( t l ) , defined by Equation 2, and the best amplitude is obtained from the value of the amplitude correlation function cu(t1). The right hand expressions are derived by partial integration and require no differentiations of F( t ) any more. I t can easily be proved that Equations 2 and 3 are insensitive to linear background. Assume the latter to be of the form:
t
Location coefficients
.-
Amp/; tude meffic!ents
en
F(t)
-
=
f
+ gt
(4)
where f and g are constants. In the absence of a peak, Equations 2 and 3 should become zero when Equation 4 is inserted. This leads to four types of integrals:
Figure 1. Example of a spectrum. Time or qualitative coordinate is plotted horizontally, intensity vertically Sampling interval is 8 p m on the photographic plate. In the lower part, deconvolution coefficients for location and amplitude determination are plotted as vectors ri+i
Figure 2. Result of location correlation applied to the previous
spectrum. Zero's of transitions from negative to positive indicate peak positions. es at 8 p m intervals. These are used to activate the analog-to-digita1 converter to take a sample and convert .it to a digital value. The computer is a Philips P 9205 with a memory space of 16 kilowords of 16 bits each. It has a cycle time of 1.6 psec and is fitted with a fast arithmetic. The spectrographic data were obtained with a Hilger E 742. The range of wavelengths was from 2400 to 3400 A.
LINE FITTING AND BACKGROUND REDUCTION Recognition of a given line profile is commonly done by application of the least squares method. When the measured function of time, t, is F ( t ) one can look for the presence of a symmetrical, but otherwise arbitrary peak of standard shape, 4 ( t ) . For reasons already explained, we shall search here for $ i l ) ( t ) in F [ l ) ( t )where the exponent between parentheses means the order of differentiation with respect to the argument. Both the position tl and amplitude CY of the standard can be varied until the best fit is found, expressed by a minimum of the integrated quadratic deviations:
+
J ; { F ' l ( t ) - cu@W tl)12dt
(1)
By differentiation with respect to the parameters t l and condition is obtained by the choice
CY,this
+
p ( t , ) = J-lF(lYt)C$(zJ(ttJdt =
+
1 r F ( t ) 4 ! 3 ) ( t t,)dt
2240
= 0
(2)
In Equations 5 and 8, the integrands are odd functions and integration yields the result zero. Partial integration is employed in (6). Provided that d ( t ) converges more rapidly than t - 2 , the constant term disappears a t both extremities and the same applies to d ( t ) . After partial integration, Equation 6 becomes similar to Equation 7 . So far, the line profile remained undefined, but from now on it is assumed that a Gaussian shape is present and its halfwidth is known in advance. A Gaussian curve easily complies with the convergence requirements. One can make up the ratio of the Gaussian function and t - 2 and, to determine its value a t infinity, differentiation may be applied to the numerator and denominator separately, after inversion of the exponents and the ratio. In Figure 1, a spectrum is given, consisting of a number of Gaussian peaks and a linear background. Preliminary, the spectrum will still be treated as if it were a continuous function, although in reality it is available in tabulated form only. When the correlation procedure (Equation 2) is applied to F ( t ) , without the restriction that the result should be zero, Figure 2 is obtained. The correlation integral is a function of the displacement tl of the standard line profile with respect to F ( t ) . The scale of tl is chosen identical with t in F ( t ) so that F ( t ) and its correlation result can be plotted in the same diagram. Thus t E tl and p( t l ) 3 /3( t ) .The function p( t ) is derived from F( t ) by the correlation operation (Equation 2). When the function P ( t ) passes from negative to positive, a peak is detected. When the position of a peak is known, its amplitude can be calculated by application of Equation 3. This is done in Figure 3 where the amplitudes of the various peaks are given by vectors. The areas of isolated peaks are correct. If the distance between two adjacent peaks is too small, the apparent amplitudes will require corrections as is shown in the section on error corrections. Interference between neighboring lines may introduce fake peaks with negative amplitude, such as a t the second zero in Figure 2. The correlation of a single Gaussian line with its rzth derivative leads to quasi n-fold differentiation, by integration. A Gaussian profile is given by
A N A L Y T I C A L C H E M I S T R Y , V O L . 45, NO. 13, N O V E M B E R 1973
Table I. Deconvolution Coefficientsa
d
2.60
+ 2.817
o m
+ 0.586
al +
-1.272 a 3 -0.428 a4 -0.059
a5 +0.006
W Result of amplitude correlation applied to the spectrum. At the zero's of the previous correlation, the uncorrected amplitudes are found Figure 3.
cr, -0.005 a7 + 0.003 Q -0.002
ag
io.001 a10 0.000
d
aj where d is the width a t half intensity and
a2
the variance.
a2 CZj
It is convenient to introduce dimensionless quantities: a 5 Q a7
r = t/a r1 = t , / a
as ag a,
Then the correlations 2 and 3 lead to n t?
-
-
+ ~ ~ ) @ ' *-' ( 77Jdr
o5 =-0.086 =
0.354@'*)(r1fi)(14)
These results are not straightforward and are explained in the section entitled Mathematics. The mutual displacement of @(T 7 1 ) and 4'3)(7- 71) is 271. The argument on the right-hand side includes the factor 2 ~ 1 / & ? that gives the extent of broadening, x n , obtained by correlation with respect to exact differentiation. Although Equations 2 and 3 give us the basic approach with which to solve the spectral data reduction, it must be borne in mind that F ( t ) is available as a sampled function only and that a numerical procedure must be found by which to perform the location and amplitude correlations 2 and 3.
+
3.00
-0.396 1.258 -0.917 0.377
-
- 0. 1m -0.018
- 0.002 0.OOO 0. OOO
3.m
3.80
- 1.214 - 0.795 - 0.181 - 0.016
- 1.122 - 1.001 - 0.952 - 1.079 - 0.282 - 0 . w 2 - 0.039 -0.075 - 0.004 - 0.008
- 0.002
+ 0.001
0.000 0.000
4.40 t2.829 + 1.745 -0.240 -1.232 1.020 - 0.480 - 0.149 0.032 - 0.005 -0.001
-
-
4.60
+ 2.828 t1.828 -0.089 -1.184 1.106
-
-0.536 -0.210 -0.053 -0.010 -0.001 0. om
0.000
@1=-1.193
e6=+a028
0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000
- 0.001
~ 2 0
+ 1.653
3.40
+2.827 + 2.828 + 2.829 + 2.829
-0.003 + 0.002 -0.001 + 0.001 0.000
+ 2.829
3.20
+ 2.824
+ 0.796 + 0.9w + 1.13 + 1.301
- 1.268 - 0.616 - 0.106 - 0.002
4.00 +,2.829 f1.545 0.53 -1.260 -0.798 0.281 -0.062 -0.009 -0.001 0.000 o.Oo0
Po = 0.000
Jr@(r
2.80
+2.800
@=
p7=
+ 0.850 0.000
- 1.173 - 0.534 - 0.128 - 0.015 - 0.001
- 0.710 - 1.232 - 0.669 - 0.197 - 0.035 - 0.004
am0 0.m 0.000
0.000 0.000 0.
-
-
- 0.001
wo
5.00
5.23
+2.817 f1.978 +0.186 -1.033 1.230 -0.790 -0.366 -0.115 0.033 -0.0% 0.001
+2.779 f2.075 +0.2w -0.913 -1.287 -0.863 -0.469 -0.149 -0.057 -0.006 -0.004
&BO
+2.826 f 1.903 +0.055 1.118 -1.175 0.693 -0.262 -0.081 -0.019 -0.003
t 1.433
- 0.861
-
-
-
c)4= f0.176 p g = 0.000
(?3=-0.308 0.000
p8=
a The amplitude coefficients en are given as a function of the half-width, expressed in sampling interval units The location coefficients are independent of the variance t
Magnitude
/
x
\
DECONVOLUTION COEFFICIENTS Let us consider next how the deconvolution coefficients can be calculated that are required for the numerical operations on the sampled functions. The spectrum F ( t ) is known a t discrete intervals only but intermediate values can be reconstructed by utilizing polynomials. A fourth order polynomial can be made to fit five consecutive sampled values. The coefficients of such a polynomial are linear combinations of these sampled values. When the functions 4 are defined as Gaussians, the correlations 2 and 3 can be executed by dividing the integration interval into a number of separate regions, in each of which a single polynomial is valid. The result consists of a linear combination of sampled values. This means that the correlations can be performed with two sets of deconvolution coefficients, one for the location and one for the amplitude correlation. The exact results of the correlations are given by Equations 13 and 14. Now one can derive those deconvolution coefficients that satisfy these equations. The location coefficients are anti-symmetric and the amplitude coefficients form a symmetric set. Let us assume that 11 independent coefficients have to be found. Then Equations 13 and 14 must be satisfied a t 11 different values of the displacement 71. In each case, this yields 11 equations and permits the calculation of all deconvolution coefficients. The convergence of the system of equations is but poor and the modification of Gauss' method has been
Figur 4. Deconvolution coefficients for the amplitud as a function of the halfwidth expressed in number of sampling intervals
used that was introduced by Crout ( 3 ) . Once the coefficients have been established, the correlation operations are quite simple. T o determine the location of a peak in a spectrum such as the one in Figure 1, the sampled values of F ( t ) must be multiplied by the corresponding location deconvolution coefficients Pn. The sum can be plotted a t the position of the center of the set, that is represented in vector form. Then the set is displaced one interval to the right and the process is repeated. When the correlation result in Figure 2 passes zero during a transition from nega(3) P. D.Crout, Trans. A I € € . 60, 1235 (1941).
A N A L Y T I C A L C H E M I S T R Y , V O L . 45, NO. 13, N O V E M B E R 1973
2241
As soon as the intensity of a line exceeds a certain value, it is classed as saturated and only the centroid is determined. If the line shapes deviate from Gaussians, the data reduction procedure can be modified.
I
Holfwidth4.5
-+Re/.
Figure 5.
freq.
Frequency distributions of an emission spectrum
Crosses are Fourier components present in the spectrum on entrance into the computer program. Curve 1 gives the theoretical distribution for the case that only Gaussian profiles are present. The amplitude and location correiations give a frequency pass band represented by curves 2 and 3. AI the output of the program, the resulting distributions are indicated by circies and triangles. Curve 4 gives the pass band for increased resolution and the squares the output frequency distribution
tive to positive, a peak is detected. To find the corresponding amplitude, the set of amplitude deconvolution coefficients a , has to be applied in a similar way and the sum thus obtained gives the area of the peak, subject to corrections given in the section on error corrections. The result of a continuous correlation is given in Figure 3. The positions of the peaks and their magnitudes are indicated by arrows. A fake line is detected but can be rejected on behalf of its negative amplitude. The coefficients are given in Table I and Figure 4 as a function of the width a t half intensity d, expressed in sampling intervals. If the number of points per peak is such that the corresponding variance falls outside this range of values, a different sampling interval can be chosen. The deconvolution coefficients for the amplitude are smooth functions of the variance and can be represented by low order polynomials. Thus, a slowly varying variance as a function of time can be tackled as well by interpolation of the deconvolution coefficients. The location correlation coefficients are not given as a function of d for reasons to be explained in the section on resolution. The number of deconvolution coefficients depends on the variance and on the required precision. Thus a 1%determination of the peak area with a relative half-width of 3.00 can be performed with 5 independent coefficients. The numerical procedure is quite fast. When no line is present, the program has to perform the location correlation only and this requires 21 multiplications a t most or roughly 200 psec. Only if a peak is detected, the apparent amplitude must be calculated. It takes 33 multiplications to obtain the interpolated coefficients. Another 21 multiplications are needed for the actual amplitude determination, totaling about 700 psec. At a speed of 0.5 cm/sec and a sampling interval of 8 pm, the time per measuring point is 1600 ysec. The amplitude correlation only gives values a t grid points 8 p m apart. When the centroid of the peak does not coincide with the grid point, a correction has to be made. Two adjacent peaks can interfere with each other and again a correction has to be applied. This is treated in the section entitled Error Corrections. 2242
FREQUENCY SPECTRUM AND BANDWIDTH LIMITATION The continuous correlation according to Equations 2 and 3 yields a quasi differentiation of the Gaussian peaks. The distortion consists of a pure broadening by a factor 1.41. The correlation technique provides smoothing of the data when the spectrum contains a good deal of noise, An abrupt increase in intensity can mean two things: the presence of a strong peak or noise. When no large peak is recorded in the vicinity, the discontinuity must be noise. A persistently constant slope of the intensity points to low frequency noise or background. To describe the smoothing aspects more quantitatively, a technique borrowed from the field of data transmission is used. The spectrum is divided into small regions of the order 3d and a Fourier analysis is made. Sinusoidal waves of various frequencies and amplitudes are fitted to the measured parts of the spectrum by the least squares method or by correlation. The amplitudes of the various frequency components are averaged over the complete spectrum, ignoring the phase of the sinusoidal waves. The distribution of amplitudes us. frequencies is shown in Figure 5 . The crosses apply to a spectrum with comparatively weak lines and are based on 5000 measuring points. Had the spectrum consisted of a superposition solely of Gaussians, the resulting distribution would also have been Gaussian and be represented by curve 1. Serious differences between the measured values and the theory are found in the low frequency region. The background obscures the information on lines in this part. To get rid of this background, the information must be drawn from those parts of the frequency spectrum that have a good signal-to-noise ratio. This can be achieved by introducing a frequency pass band or filter. Curve 2 illustrates the filtering action of the program. It gives the relative reduction of amplitude of the Fourier components in the spectrum after passing through the amplitude correlation of the program. A new distribution results, indicated by circles. The most important frequency components are all situated around 0.07. As very high and very low frequency components are absent, the inverse Fourier process can lead only to smooth curves with a mean value equal to zero. A more concrete mathematical form can be obtained by replacing the function F ( t ) in Equation 3 by a cosine wave of frequency W. It is symmetric with respect to the origin. Asymmetric sine waves do not contribute to the Fourier spectrum. Amplitude correlation 2 is applied to the cosine wave. Taking 40)
-
r2 = wa
cos 727$'*'(7)d7
-
(15) (16)
This result gives the amplitude of the cosine wave that, after differentiation, gives the best fit to the first derivative of a Gaussian, situated a t the origin. The resulting function is strongly frequency-dependent and the quadratic term serves to reduce the influence of the lower frequencies while the Gaussian reduces the higher frequencies. The result obtained above is not so straightforward and is explained in the mathematics section. What has been said about the amplitude correlation can be repeated for the location correlation. The frequency
A N A L Y T I C A L C H E M I S T R Y , VOL. 45, NO. 13, N O V E M B E R 1973
Origimirpectrum
intensity
'
'\ // /'
A,
'.
Peak 1
Peak2
-Time
Figure 7. Instrumental resolution defined by a 50% valley
Carrebted with second derimtivc
Figure 8. Limit of resolution by valley Figure 6. Noise reduction by correlation Circles give the sampling values of the synthetic spectrum with a high noise level. The dots give the result of amplitude correlation. The heavy line gives the contribution of the hidden Gaussian line
pass-band and distribution of the output signal are given by curve 3 and the triangular points. The noise reduction of the correlation technique is illustrated in Figure 6. In the upper part, a synthetic spectrum is indicated consisting of a linear background, a single Gaussian line profile, and random noise with a white frequency spectrum. With the aid of correlation 3, performed with the corresponding deconvolution coefficients, an inverted quasi second derivative of a Gaussian is obtained, indicated by dots. Without noise, the drawn curve would have been obtained. In spite of the high noise level, the area of the hidden Gaussian can be obtained with a spread of roughly 10%. If the original spectrum is differentiated twice, the contribution of the peak is less than 10% of the noise and the former can no longer be recognized, unlike the case of the correlation result. A frequency analysis of a spectrum can also be made without dissecting it into small regions but this brings about adverse effects. If the spectrum consisted of a repetitive pattern, a large frequency component would be found a t the repetition frequency. Such effects of long range correlations obscure the actually wanted information on line shape; that requires a correlation over the width of the line only. The deconvolution coefficients cover a limited range at a time and do not involve long range correlations.
RESOLUTION The resolution of a spectrum is primarily determined by the spectrometer. Whether the information produced by the apparatus is fully exploited, depends on the subsequent treatment of the data. This can be illustrated by a number of cases given in Figures 7 to 10. The Gaussians of equal magnitude and equal halfwidth d are gradually merging. In Figure 7 , the distance is 1.41 d and this corresponds to instrumental resolution defined by a 50% valley. Visual inspection will still reveal two separate lines until the valley disappears a t a distance of 0.85 d. This is by no means the limit of resolution. If the intensity spectrum of Figure 8 is differentiated twice, the heavily drawn curve of Figure 9 is obtained in which two peaks can clearly be distinguished ( 4 ) . Finally, when the distance is decreased (4) J. R. Morrey, Ana/. Chem.. 40, 905 (1968)
v u
Figure 9. Zero valley case resolved by second derivative
M k l
h k 2
-
Vine
Figure 10. Limit of resolution by second derivative
to 0.64 d as in Figure 10, the limit of resolution for this method is reached. In the absence of noise, one could continue looking a t increasingly higher derivatives to improve resolution. The central lobe of the higher derivates continues to slim down but new side lobes are generated and the total width increases. This causes additional interferences that obscure the picture. In practice, noise imposes even stricter limits. As the order of differentiation increases, the high frequency noise is enhanced and becomes prohibitive unless bandwidth limitations are introduced. Correlation techniques such as Equation 2 do this, but result in a broadening by a factor 1.41; the resolution is limited to an inter peak distance of 0.90 d and no improvement is obtained with respect to the visual method. To increase the resolution, the following trick is applied. Instead of correlating the measured spectrum F ( t ) in Equation 2 with the third derivative of a Gaussian of prescribed variance, a smaller value of d is used; the width a t half intensity is given by 2.3556 b. By analogy with Equation 13, the following result is proved in the section on mathematics:
A N A L Y T I C A L C H E M I S T R Y , V O L . 45,
NO. 13,
N O V E M B E R 1973
2243
3
f
One
it+
Two peaks
c
t -
-tr-tp I
0
1
2
,
'
I
' 3 ' 4
1
I
5
6
-w90 1
7
I
8
- .- - - - - - - - - - - - - - - -
-W%
Figure 12. Error correction for overlapping lines In the upper part, two overlapping lines are drawn. The lower curve gives the contribution to the amplitude correlation of an interfering iine at a distance I t and with an area of unity. The horizontal axis gives the reduced distance
Figure 11. Experimental resolution and precision Along the horizontal axis, the distance between two gradually merging peaks is plotted. In the upper part, the areas of the peaks, determined by the program are plotted; in the lower part, the relative errors in the positions. The point pairs give separable peaks. The crosses apply to unresolved doublets
resolution, as it leads to a poor signal-to-noise ratio. The program will find more peaks but with less confidence. The subsequent amplitude correlation does not suffer from this drawback, and erroneous peak positions will be corrected by attributing small amplitudes to them. The location correlation result obtained from Equation 17 passes from negative to positive a t a position that is independent of a and b. Thus it is not really necessary to adapt the location correlation coefficients to the actual variance of the peak.
ERROR CORRECTIONS The variance of the correlation result is greater than both that of the function used for the correlation and that of the function which is correlated. When b = 0.5 a, the broadening is only 12%. Thus a minimum distance for resolution of 0.72 d is found. Experimentally a resolution of 0.78 d is achieved. This is shown in Figure 11 where there are two gradually merging lines with Gaussian profiles. The halfwidth is chosen as 2.8 and the area as 150. In the upper part of the diagram, the values of the peak areas, found by the computer, are given as a function of the separation of the peaks. Up to a point corresponding to a 50% valley, the deviations from the value 150 stay well within the 10% error area, bounded by streaked lines. As the valley disappears, the errors start scattering. The case drawn in the insert is indicated by the heavy arrow and it still gives reasonable results. There is a minimum separation of 2.2 units beyond which the separation criterion breaks down. Then the two peaks are treated as a single one and the errors become large. Finally, when both peaks coincide, the value 300 is correct again. In the absence of sufficient resolution, nothing better can be done in the intermediate range. The errors in the determination of the position are shown in the lower part of the diagram. Actually, the dislocation divided by the variance of the line is plotted. As long as two separate peaks can be found, the error stays within the 10% limits. Afterward the centroid of the conglomerate is traced. The influence of the improved resolution on pass band is illustrated in Figure 5 , curve 4. It is twice as broad as curve 3, and the influence of higher frequencies is apparent from the square output points. Those frequencies in the input spectrum that do not contribute to the assumed Gaussian profile have about as much weight as those that do contribute. This illustrates the hazards of improving 2244
The deconvolution coefficients yield results a t points of the measuring grid only. The zero value of the location correlation can be determined with sufficient precision by linear interpolation. The corresponding amplitude correlation can be applied only to the nearest grid point. But for a normalization to unity, the error made follows Equation 14 and is plotted in Figure 12, lower part. The horizontal axis gives the displacement At of the correlating function with respect to the actual peak center, divided by the local half-width. The curve is independent of half-width. If the distance from peak 1 to nearest grid point 3 in reduced quantities is A t / d , the error can be read from this correction curve. The result is always a reduction, cu(1,3), and to arrive a t the actual value of the peak area, the correlation result must be divided by a (1,3). The amplitudes of overlapping lines can be corrected with the aid of the same curve. In Figure 12, an example is given. The overlapping peaks are situated a t 1 and 2 and the nearest grid points are a t 3 and 4. The areas of the lines are k(1) and l ( 2 ) . With this knowledge, the apparent amplitudesp(3) and q ( 4 )can be calculated:
+ cu(2.3)1(2) + cu(l.4)k(l)
p ( 3 ) = ai4,3)h(l) q ( 4 ) = a(2,4)1(2)
(18) (19)
Solving these linear equations leads to the required inverse process: k(l) =
a ( 2 . 4 ) p ( 3 )- a(2,3)q(4) at1.3)a(?,4)- a(1.4)0((2,3)
120)
It is possible to extend the correction procedure to three or more peaks interfering with each other. The required
A N A L Y T I C A L C H E M I S T R Y , VOL. 45, NO. 13, N O V E M B E R 1973
'1
Re/. ampi.
&
Figure 13. Influence of the width of the plate reader slit on the amplitude correlation W h e n s l d I S increased beyond t h e value 2 t h e presence of a double peak suggests itself
resolution and precision limit the process to three adjacent peaks. When, instead of a correlation with the second derivative, a higher order of a Gaussian had been used, the correction curve would have been wider and the interferences would have multiplied. As the correction curve can easily be stored in the computer memory as a table or a polynomial, no iterations are necessary and the corrections can be applied quite rapidly.
WIDTH OF PLATE READER SLIT So far it has been assumed without further comment that with a certain width, 5, of the plate reader slit, a Gaussian ensues when the transmittances are properly converted to intensities. Increasing the width leads to a distortion of the Gaussian peak and reduces the higher frequencies of the Fourier spectrum. This can be shown by mathematical treatment. Assume that a Gaussian profile is found in the limiting case of nearly zero slit width. When s is increased, an average is measured between the limits t - to and t + t o where s = 2to. Taking T = t / a and T O = t o l a , the new intensity function is
1-+-
~ ( T . 7 0=
@(TICIT =
+
erf(.r
To)
- erf(.r
- To)
I,
(22) in which the error function is defined by
d erf(r)/d.r
(23)
= $(T)
This new intensity function has to be correlated o((T1,,T1)
- Jm -
(7
T Z
{erf(.r
+ + To
TI)
erf(7 sinh no-
-
TO
-
2
70
cosh
+
-
T ~ ) \ ~ ( ~ ) ( Tl)dT T
TTO
-
This is explained in the section on mathematics. In the 0 and thus T O 0, expressions 24 and limiting case s 14 become identical, but for a difference in amplitude. The influence of the slit width on the distortion of the correlation result is given in Figure 13. When s / d is increased to 3, a maximum flanked by two minima is found, suggesting a double peak and the results become erroneous. I t is safe to chose s equal to the sampling interval. In the calculation of the frequency pass band of the amplitude correlation, only cosine functions are required. Defining 7 2 = o u , the integration by a finite slit gives a cos
T ~ TdT
= 2 a ~ ~ -cos l r 2 T sin T ~ T( 2~5 )
Figure 14. Reduction of the frequency pass band of the amplitude correlation by the width of the plate reader slit W h e n s / d = 0, the result
I S identical with Figure 5 T h e h i g h frequency contributions are already severely curtailed w hen s l d = 1
and when this function passes the amplitude correlation, according to the mathematics section, the result is:
-
-JAX.)TI-I cos r27 sin T2T04'2'(T)dT
T?@(T?)
sin
7270
-
(26) In the limiting case T O 0, this result has the same shape as Equation 16. The influence of the finite slit width is contained in the factor 72-1 sin 7270 in Equation 25. The frequency pass band of the location correlation is based on sine functions: -J'?TL~
sin
T2r
-
sin T ~ T ~ ~ ' ~ ' ( TT2 2) 4~( TT2 ) sin
7270
(27) As can be seen in Figure 14, the wider the slit, the better the high frequency noise suppression. As the resolution is also affected, a good compromise must be found. Although the amplitude correlation depends on the value of s, the error made could be corrected by correlating with a distorted Gaussian function, obtained by the same integration procedure. Then new deconvolution coefficients have to be calculated along the same lines as before.
ERRORS RESULTING FROM INSTRUMENTAL INSTABILITIES When the variance of the lines as a function of the qualitative coordinate, time, has been established, this vital information must be exploited. Methods introducing a fitting procedure of the variance do not fully exploit the available knowledge. Emission spectrographs and kick sorters in y-ray spectroscopes have sufficient stability, but a mass spectrometer cannot always be expected to have the stability required over a period of hours. The errors resulting from a poor guess a t the variance to be expected can be derived from the result of the correlation
("8)
where a2 is the expected variance and b2 the actually present variance. When b = (1 - ?)a, y is the fractional
A N A L Y T I C A L C H E M I S T R Y , VOL. 45, NO. 13, N O V E M B E R 1973
2245
MATHEMATICS
Error in
amplitude
L
\
For the reader interested especially in the mathematics, a number of hints are offered that may facilitate the calculations. In the correlation integrals, the product of two mutually displaced Gaussians or derivatives of the latter is frequently encountered. By straightforward calculation and starting from the definition of a Gaussian, one proves the statement
+M%
-50%
-50%
Figure 15. Errors caused by erroneous prediction of variance 1
2
intensity
Making use of this result one obtains
(33) Figure 16. Precision of
-t position determination
Two peaks a and b are so close that resolution fails. When the ratio of the areas is changed from 3 : l to 1 . 3 , the envelope shows a shift of the top of the order of resolution
in which we exploit the property that the area of a Gaussian is not affected by a shift of position. By repeated differentiation with respect to the parameter tl and application of theorem 35 one derives
error in halfwidth. Comparing Equation 28 with the result for a = b, one arrives a t the fractional error in the area determination, E,
E
= (1
+ y + j/?y')-'/*
(29)
a function plotted in Figure 15.
LIMIT OF ACCURACY OF LINE POSITIONS The precision of determination of the position of a line on a photographic plate is Zt2 pm. The half-width of the lines is of the order of 30 pm. Lines can be resolved when the distance between equally large peaks is about 20 pm. There seems to be a discrepancy between the values 2 and 20. When looking for possibly interfering element lines in qualitative analysis, one has to consult tables such as ( 5 ) . Taking for granted a determination of the position with a precision of h 2 pm, corresponding to a wavelength variation of k0.008 A in a specific emission spectrograph, the number of interferences would amount to 1/3, as a mean. Assuming an uncertainty of the order of Zt20 pm or Zt0.08 A, an average of three possible interferences is found. Thus the first case is the most profitable. Yet, this case must be discarded as being overoptimistic. Let us assume the line in question consists of a doublet that has not been resolved because the distance between the components is 20 pm or less. This case is given in Figure 16. When the ratio of the line intensities is inverted, the position of the top of the conglomerate is seen to shift over a distance of the order of the resolution. Thus the uncertainty of the position is of this same order of magnitude. The situation would change if one knew that the line to be examined was free from interference. Then one could borrow information from instruments with a higher resolution, such as exact wavelength. A series of lines all lying a t exactly the wavelengths given by high resolution instruments would suggest freedom from interference of all these lines and would present a basis for the r t 2 pm accuracy. ( 5 ) W. F. Meggers. C. H. Corliss, and F. Scribner. Tables of Spectral Line Intensities, Nat. Bur. Stand. (U.S.1Monogr.. 32, Part 1 1 (1961): U.S. Government Printing Office, Washington 25 D.C.
2246
Partial integration proves the theorem
For the determination of the limits of resolution, the sum of two displaced Gaussians and their derivatives is needed.
S(7) =
(T .
+
+
T ~ ) +(T
-
-
S'"'(0)
7J
- @(TI
cosh
(T~T)+(T~)
(36) (37)
.'n'(71)
The calculation of the distortion caused by a finite plate reader slit leads to the difference of two Gaussians.
D(T)=
(T .
+
71)
-
-
$(T
71)
-~XT)
sinh
(38) The Fourier analysis uses a set of formulas derived from the Poisson integral by inserting a complex variable 1 = S*;o(r!dr = =
~ I @+( T+ T~
i7&d7
(39)
When the real parts are equated, the result is
SI@(.;) cosh
T ~ Tcos T?T
d7 =
{ C # J ( T ? ) / @ ( T Ocos )~ T
~
T
(40) Again, differentiation with respect to T O or T Z introduces new formulas required to solve Fourier problems. In general the nth deviations of 4 ( ~ can ) be expressed as a polynomial in T , of order n, times @ ( T ) . In the next section, it is stated that a Gaussian is a solution of a diffusion process. When the concentration of counts is given by the function c ( x , t ) where x denotes the channel location
A N A L Y T I C A L C H E M I S T R Y , V O L . 45, N O . 13, N O V E M B E R 1973
~
c ( x ,t )
= t-1/2 exp(-x2/4Dt)
W
(41)
bl
it satisfies the diffusion equation with D as the diffusion constant.
Dd2c/dx2 = b c / d t
A
(42)
The total concentration given by Equation 41 is constant. If the time decreases, the magnitude of the Gaussian increases as the variance decreases. All counts are swept back to the center of the peak. This is time-reversed diffusion. Normal diffusion results when time progresses in a positive direction.
COMPARISON WITH OTHER METHODS Inouye, Harper, and Rasmussen (6) have described a data reduction program that relies essentially on Fourier transformation, filtering in the w space, and inverse Fourier transformation. The filtering results in a noise reduction and distortion of the lines. The authors are a t liberty to choose any type of filter. A filter with a very narrow band, however, resembles a resonant circuit with a low damping. It introduces line satellites. Although a special and fast Fourier transform program is used, about 20 msec are required per measuring point; a 25-cm long spectrum sampled a t 8 p m takes 10 min to decode. This dramatic increase in computer time does not pay off in terms of increased signal-to-noise ratio if the resolution is comparable to the program in our paper and no side lobes are introduced. The peak identification program published by Mariscotti (7) is, a t first glance, completely different from this work. He employs a smoothing operation that in principle, but not in practice, averages the measured spectrum, or its derivatives, in a number of successive steps z . The averaging is performed over a number of adjacent elements UI. Experimentally the values u’ and z are optimized. Apparently no line shape is assumed. This is not true, however, as can be shown by a certain reasoning, not given in Mariscotti’s paper. At the same time, the averaging procedure can be made more translucent. Let us assume that a single y emission line is studied with a multichannel counter and that 3125 counts are integrated in the central channel and, by some miracle, none in the other channels. Next, the counts are allowed to jump 0, 1, or 2 channels in either direction and with equal probability. This process can be repeated a number of times, say 5 . The changes in distributions are given in Figure 17. The final distribution is also displayed in graph a . Its shape resembles a Gaussian. Not by coincidence, however, for the process described above is a diffusion process of counts serving as particles. Such a diffusion inevitably leads to a Gaussian profile. Had one chosen u: = 3 instead of 5, then the mobility of the counts would have been lower and a greater number of steps would have been necessary to arrive a t the same variance of the distribution. When a second derivative of the curve u: = 5 , z = 5 in Figure 17b is obtained by finite differences, the distribution is the same as the deconvolution coefficients used in Mariscotti’s paper. This means that he, too, essentially employs the technique of least squares to recognize line shapes corresponding to the first difference of the distribution u: = 3, z = 5 to the first differences of the measured spectrum. In fact he has to correlate these functions
c) 3125 625 625 625 625 500 625 500 375 250 125 450 475 450 3375 250 150 75 m u s NO 340 a60 in 100 a)-c365 381 365 320 255 185 121 16 -16 -4 -65 -m -tx -51 b)+-a -32 -29 -20 -5 6 13 c)-3 3 9 15 11 7 3
a5 50
20 s 70 35 15 5 -35 -20 -10 -0 10 6 16 15 -1 -5 -4 -3
I -1
3 1 -2 -7
Mariscotti’s averaging procedure and generation of line profiles in y-ray counting techniques
Figure 17.
When 3125 counts in a multichannel counter are allowed to diffuse by making five successive jumps over not more than 2 channels, a distribution results that is illustrated in a ) . Differentiating twice and thrice leads to b ) and c ) . The center of the distributions i s indicated by the heart symboi
or the measured spectrum with the second differences of the standard shape. This should be compared to Equation 3. The averaging process generates a line profile fitting best with that measured if u! and z are chosen accordingly. The amplitude correlation is almost identical with ours. The position determination is not performed in the same manner though it could have been done so with success. The location correlation is an essential part of the least squares method. These coefficients are displayed in Figure 17c. Helmer, Heath, Putnam, and Gipson (8) have written a program for y-ray counting that uses Gaussian fit with adaptable variance. The fitting procedure becomes nonlinear and the resulting iterative process would take about an hour, if a spectrum containing 360 lines were to be processed. Such a program is well suited to finding the function variance us. energy. Once it has been established, the adaptation of the variance means a loss of information, because previously obtained knowledge is discarded. The increased computer time does not balance the gain of information. Faced with automatic evaluation of gas chromatograms, Westerberg (9) made use of digital filtering techniques that can be simulated by computer programs. After discarding crude techniques to resolve overlapping peaks, he resorted to curve fitting. A drawback of gas chromatography is the irreproducibility of the line profile. A number of modified Gaussians were considered. Mills ( I O ) advocates a variant to Mariscotti’s method of averaging in which the number of successive steps, z, is not given in advance, but decided by the operator after having reviewed the process. A lower number of z leads to correlating with a Gaussian of smaller variance and to increased resolution, as is explained in our paper. Such a process is time consuming and places the burden of decision on the operator. In a stable gamma-ray counter, this is quite unnecessary, once the right value of z has been chosen. In this case too, the position determination should follow the same lines. Received for review March 8, 1973. Accepted May 29, 1973. (8) R G Helmer, R L Heath M Putnam, and D H Gipson, Nuci In-
( 6 ) T. Inouye. T. Harper, and N. C. Rasmussen, NucI. instrum. Methods. 67, 125 (1969) (7) M. A. Mariscotti, Nuci. Instrum. Methods. 50, 309 (1967).
strum Methods 57. 46 (1967) ( 9 ) A W Westerberg, Ana/ Chem 41, 1595 (1969) (10) S J Mills Nuci instrum Methods 81, 217 (1970)
ANALYTICAL CHEMISTRY, VOL. 45, NO. 13, NOVEMBER 1 9 7 3
2247
