- Email: [email protected]
Numerical methods for the correction of apparent band shapes due to finite slit width
Spectrochimica Acta, 1962,Vol. 18,pp. 1629to 1639.PergamonPressLtd. Printedin NorthernIreland
Numerical methods for the correction of apparent band shapes due to finite slit width A. L. KHIDIR
and J. C. DECIUS*
Department of Chemistry, Oregon State University, Corvallis, Oregon (Received 7 &Iuy 1962) Abstract-The problem of transforming apparent spectroscopic band shapes to their true values is considered. The integral equations which express the influence of a spectrometer with finite slit width are given a matrix representation in which the problem is reduced to the inversion of Although an exact solution of the the matrix which represents the spectrometer transformation. inversion problem is readily found for a wide class of slit functions, it is shown to be impractical in the sense that a small experimental error in the observed apparent intensities leads to a much larger error in the true intensities calculated with the aid of the inverse spectrometer matrix. Nevertheless, it is demonstrated that a satisfactory approximate solution can be found by an iterative calculation which converges sufficiently rapidly provided the slit width is not too much greater than the band width, and that the band does not contain extensive fine structure. The method makes no assumptions concerning the true shape of the band and may, therefore, be used in the case of asymmetrical bands. Neither is it restricted to any particular shape of the slit function (e.g. triangular, trapezoidal, gaussian). Several numerical examples are used to illustrate the flexibility of the method and the speed of convergence to be expected in practice.
THE correction of errors in spectrometric data due to the finite resolving power of real spectrometers is of special importance in connection with infrared absorption measurements of the so-called absolute (integrated) intensities. The techniques introduced by WILSON and WELLS [l] are widely used in the investigation of gases. In the case of condensed phase samples, it is possible to use extrapolation procedures or to employ tabulated corrections to the maximum absorbance, the band width, or the integrated absorption coefficient as given by RAMSAY [2]. The latter procedure suffers from the disadvantage that it is premised on the two assumptions (i) that the true band shape is Lorentzian and (ii) that the spectrometer transmission function is triangular, while the former procedure makes it difficult to detect concentration effects on the integrated intensity arising from solute-solvent interactions. It is, therefore, of interest to consider whether more flexible methods can be devised to make such corrections. A very general solution of this problem has been given by HARDY and YOUNG [3], who, by the use of the theory of Fourier transforms, obtain expressions for the true intensity at any point in the spectrum as a series which involves the apparent intensity and its derivatives of all orders at that point; the properties of the spectrometer * Alfred P. Sloan Foundation
Fellow.
[l] E. B. WILSON, Jr. and A. J. WELLS, J. Chem. Phys. 14,578 (1946). [2] D. A. RAMSAY, J. Am. Chem. Sot. 74, 72 (1952). [3] A. C. HARDY and F. M. YOUNG, J. Opt. Sot. Am. 39, 265 (1949). 1629
A. L. KRIDIR and J. C. DECIUS
1630
are expressed through a set of coefficients for the derivatives of the apparent absorption curve. Despite the generality of HARDY and YOUNG’S solution it does not seem to have been widely used in practice. The same fundamental difficulties which impede the application of their technique will limit the use of the method which is described in the present paper, which may, however, have the advantage of exposing these difficulties in a more explicit fashion. General discussions of these difficulties have been given by RAUTIAN (4). The mathematical technique to be described below fulfils the aim of flexibility, and the tedium of the computation required can be mitigated by the application of There are, nevertheless, some serious difficulties an automatic digital computer. involved, even when automatic computers are applied to the problem. Although an exact formal solution can be found for a wide class of spectrometer transforms, this solution proves to be without direct value, since it is excessively sensitive to error in the experimental data. The exact solution is nevertheless of value in suggesting convergent iterative methods and the conditions under which they may be usefully applied. The relation of these methods to the basic problem of resolution is discussed, and some plausible examples are worked out numerically. The discussion is restricted to the case of bands which are practically isolated in the spectrum. GENERAL SOLUTIONOF THE SPECTROMETERTRANSFORMATION Matrix approximation The influence of the spectrometer on the true incident energy spectrum, expressed in the form of a transmitted energy spectrum, t(d) as a function nominal spectrometer frequency v’ where t(v’) = p3(v’,v)
i (v) dv
i(v), is of the
(1)
in which the spectrometer transformation is expressed by S(v’, v). A formal solution of our problem is obtained if an inverse transformation function, X-l(v”, v’), can be found which satisfies js-l(v”, v’)S(v’, Y)dV’ = B(vN,v) (2) in which B(v”, v) is a Dirac function, since then i(v”) = js-l(v”,
v’)t(v’) dv
(3)
In practice of course t(d) is not given analytically, but is known graphically and can be tabulated from the observed data for come convenient set of discrete values of v’ . This suggests that (l), (2), and (3) should be replaced by the system of matrix-vector equations : t = Si (I’)
S-IS = E i = S-It [4] S. G. RAUTIAN,Sov. Phys. Uspekhi 66, 245 (1958).
(2’) (3’)
Methods for the correction of apparent band shapes
1631
in which the components of the vectors i and t are simply the values of the incident and transmitted spectra at a set of discrete frequencies: %?I= ;(Y,);
(4)
t, = t(d,)
In such a representation a spectrometer of unlimited resolving power would have a matrix S which would be a unit matrix. A spectrometer of limited resolving power would be represented by a matrix with off-diagonal terms. Cyclic form of the spectrometer matrix, S It is convenient in several respects to make the assumption that S is cyclic in the sense that each row is a cyclic permutation of the first row, i.e. that S has the form
S=
‘abed
....
babe
.......
cbab
........
.
.
.
.
bed....
..dcb
-
dc d (5)
. dcba
For one thing, such a form for S ensures that if i, is constant, ization of S is achieved by the requirement &
= 1
so is t,.
Normal-
(6)
In practical examples to be given below, a few of the elements a, b, c, d, etc. near the diagonal will be taken as different from zero, with many vanishing elements farther from the diagonal. If we suppose in particular that there are four nonvanishing elements, a, b, c, d, in the notation of (5), with N, say of the order of 20 to 50, we must justify the inclusion of the non-vanishing elements d, c, b in the upper right-hand corner of S. The first row of S will ordinarily be applied to calculate the intensity transmitted by the spectrometer at a frequency vl’, i.e. at the low-frequency edge of the band. The elements in the upper right-hand corner of S represents a fictitious transmission of energy arising from incident frequency components $N_2, iN_l, i, at the high-frequency edge of the band. But provided the interval of frequency and N are chosen properly, the energy in the region of iN_z, i,.,_i, i,,, will just compensate for the energy which is lost in consequence of the finite-frequency range. In other words, at spectrometer setting vl’, the contribution to the transmitted intensity by real frequencies less than v1 = vl’ 1s ’ compensated, provided i(v_J = i(vN), i(v_J = i(~~__~),etc. This condition will be satisfactorily fulfilled provided the incident spectral intensity does not change rapidly over the width of the absorption band. Later we shall introduce an iterative method of solution which does not require that S be cyclic, but it has been found convenient to retain the cyclic form both to demonstrate the convergence of the iterative method and to conserve memory space in calculations performed with an automatic digital computer.
A. L. KHIDIR and J. C. DECIUS
1632
Inversion of the spectrometer matrix We are now ready to reap the main advantages consequent from the assumed form of S. Our fundamental problem is simply the calculation of the inverse of the S matrix. This will be done by appealing to the diagonal form of 5. Any cyclic matrix of the assumed form may be diagonalized by a similarity transformation with a unitary matrix, U,
D = USU-l
= USUt
where the dagger indicates transpose conjugate. given by
(‘1
In particular, the elements of U are
in which N is the number of rows and colums in S, and k and m will be assigned the ranges 0 to N-l and 1 to N respectively. The elements of D are D,,,
= d,,, D(k) = a + 2b cos (27&/N) + 2c cos (2~2Tc/N)
(9)
+ 2d cos (2~*31c/N) + etc. Then since
(D-l),,
= &,/D(k)
= 6,,.D-l(k)
and
(10)
S-l = Uf D-l U
(11)
the elements of S-l are easily calculated and have the form I s&l = ; 1 + 2 cn_d”eD-l(k) cos 2nk(p - q)/N k=l C where N is odd, or A,-, = ;
1 + 2
(N--2)/2 2 D-l(k) cos 2vk(p k=l
1
q)/N + (-l)“-*D-l
(12)
(N/2)
1
(13)
where N is even. For large values of N the calculation required in (12) or (13) may be carried out very simply on a digital computer ; we have, in fact, programmed such a computation for the Alwac-III of the Department of Mathematics at Oregon State University.* Resolving power an& singularity of the S matrix One may well raise the question whether the general solution of the transform problem is capable of regaining the information presumed lost due to the finite resolving power of the spectrometer. Put differently, if the inversion of S expressed by (12) or (13) remains valid as N --f co, which is to say that Av = v,+~ - v, -+ 0, then in principle the true energy incident upon the spectrometer, i(v), could be calculated from the energy transmitted by the spectrometer, t(v’), with unlimited resolution, provided, of course, that one knows the spectrometer transform function, S(v’, v). * However,
see Table 1.
Methods for the correction of apparent band shapes
1633
But suppose that S is singular, i.e. that for some 4 = S&/N, D(4) = D(k) = a + 2b cos + + 2c cos 2+ + etc. = 0. Then, of course, 8-l will not exist, and our method will fail. It is easy to see how this will in fact occur for some reasonable spectrometer functions. Consider, for example, the triangular silt function. Designate by s the half width of the slit in units of Av and for simplicity take s as an integer. The constants a, b, c, cl, etc. are then given by a = l/s = s/s2 b = (s -
l)/s2
c =
(s -
S)/s2
d = (s -
3)/s2
(14)
etc. Then S will be singular if s20(+)
= s + 2(s -
1) cos 95 +
2(s -
2) cos 2+ + etc. = 0
(15)
This will in fact be the case if k = N/s
(16)
i.e. if N is divisible by s. This will surely occur in the limit as N + co. Even in the case of a finite approximation where N and s are chosen so that N is not divisible by s, one may expect to encounter near singularities when both N and s are reasonably large, and therefore great caution must be exercised in performing the calculation of (12) or (13) to avoid undue loss of arithmetic significance. As an example, we have calculated S-1 for the triangular slit with N = 41 and s = 3. The eigenvalues of S do not exceed unity, but for k = 14, D(k) = 0*00084365 so that D-l(k) is known to only about four significant places even though the calculation is done originally with 7 or 8 places. The values of D(k) and of SF, are given in Table 1. This, however, is not the most serious difficulty. In practice, the exact* solution turns out to be essentially useless for the following reasons. The elements of S-1 are all large compared with unity, at least for practical values of N such as N > 20. Under these circumstances, the errors in the calculated values of i become large compared with the errors in the observed values of the t. Specifically, if the errors in the t are approximately O-1 per cent the errors in the calculated values of i will be approximately lo-15 per cent. Put somewhat differently, the calculated values of the i are linear combinations of all the t, of a form which involves the approximate cancellation of many large numbers of opposite sign. Thus an unbearable burden is put upon the experimentalist, since no spectroscopist-at any rate none working in the infrared region with conventional equipment-will relish the task of recording an absorption spectrum with the O-01 per cent intensity precision which would be required in the above example to yield a final result good to 1 per cent. We may now consider a different approach, which still utilizes our knowledge of the diagonalization of S. Consider the difference i -
t = (E - S)i
(17)
* Here and elsewhere we mean by “exact solution,” the exact inversion of the assumed S matrix, which is of course itself only an approximation,
of finite resolving power, to (1) and (3).
1634
A. L. KHIDIR and J. C. DEGIUS
and, assuming the difference to be small, suppose that, starting with the initial approximation i. = t on the right-hand side of (17), the iterative process defined by
i, - t = (E - S)i,_,
(18)
Table 1. Eigenvalues of S and elements of S-l for triangular slit with N = 41 and a = 3 Eigenvalues
Inverse matrix, S-l
k
D(k) = D( -ii)
P-9
SZ
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1~00000000 0.98443472 0.93882698 0.86634570 0.77197224 0.66207677 054388195 0.42486048 0.31212103 0.21183997 0.12878662 0.06598623 0.02454835 0.00367495 0.00084365 0.01214789 0.03275841 0.05746238 0.08122593 0.09972774 0~10981087
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
s3* -43 -34 74 -43 -25 65 -43 -16 56 -43 7 47 -43 2 38 -43 11 29 -43 20
* Candour compels us to admit that the SG values were calculated numerically using (12) ; the fact that S,-,l are integers was noted empirically and verified by checking that SPIS is indeed the unit matrix.
will converge. imation is
In fact the difference
between
the n-th and (n -
l)-th
approx-
i, - in_, = (E - S)“t so that if t is expanded in the eigenvectors u tk) of E - 8, which are identical with those of S and whose components are defined by (8), one finds i, -
i,_l = (E - S)nz tkutk) = 2 t,jl,~(~) k
where 1, is an eigenvalue of E -
(19)
S and is therefore given by Ak = 1 -
D(k)
(20)
(see equation 9). Thus provided D(k) f 0, the difference between successive approximations to i approaches zero as n goes to infinity. Clearly the difficulty in this
Methods for the correction of apparent band shapes
1635
process in the slowness of convergence associated with very small values of D(k). Moreover, if it proves necessary to carry the iteration in (18) to a large value of n, it can easily be shown that the result will be equivalent to the use of S-l which we have already shown to be impractical. However, if t, in (19) is negligible for values of k such that 2, is very close to unity, then the iteration may converge with practical accuracy for sufficiently small values of n thereby yielding a useful result. We have found this to be true for several reasonable examples. These examples have two characteristics : (1) the band is reasonably smooth and (2) the slit is not more than a few times the true band width. Despite the success of our method in these cases, to be described in greater detail below, one can anticipate difficulties if one pushes the procedure too far, for instance in attempting to resolve extensive fine structure, since in such cases the t, for higher-frequency Fourier components would be larger. ITERATWE INVERSION OF SOME EXAMPLES Although zian, i.e.
we have chosen as some sample bands i(y) which are either Lorenti(y) = exp --E(V) cc(Y) =
(21)
a (Y -
(22)
YO)2+ b2
or Gaussian, U(Y) = a exp [-(Y
-
aQ2/b2]
(23)
it must be emphasized that our method is in no sense restricted to these cases and indeed was developed specifically to avoid this restriction. Also, we have used as a spectrometer function the two forms of (i) a triangular and (ii) a Gaussian transform, but this restriction is also not essential. In this connection it is interesting to note that KOSTKOWSKIand BASS (5) in the course of a related investigation have demonstrated experimentally the adequacy of a Gaussian spectrometer function in the near infrared. The program for the iterative calculation is based essentially upon (18), except that in order to accelerate convergence, those components of i, which have already been calculated during the n-th cycle are used in place of the corresponding components of i,_, on the right hand side of (18). Each cycle is initiated at a wing of the band, and the program has not taken advantage of the symmetry of the bands chosen for examples, since we desired explicitly to treat asymmetric bands. Experience has shown that a desk calculator would be adequate for the slightly distorted functions t which arise when the silt width does not exceed one-half of the apparent band width. Triangulur slit ; Lorentzian band Table 2 illustrates the results of the numerical calculations along with Fig. 1. The true band shape was based on the following parameters: a = 4, b = 2 in (21). The values of i, and t, are given for p and q running from 1 to 32 ; p = q = 16 corresponds to the band centre, which amounts to 36.79 per cent transmission. [5] H. J. KOSTEOWSEIand A. M. BASS, J. Opt. Xoc. Am. 46, 1060 (1956).
A. L. KHIDIR and J. C. DECIUS
1636
The t, were calculated using a triangular S function with s = 3. Thus the apparent band shape itself involves the approximation due to taking discrete values of the data, and our numerical test of the inversion in which we seek to back calculate i is a Table2. Lorentzian band; transform table a = 4,6 = 2,Av = 1,s = 3 i
t
i-t
1”26
L - t2f3
0.9828 0.9803 0.9770 0.9733 0.9683
0.9771 0.9733 0.9685
0.9819 0.9797 0.9765 0.9725 0.9674
+0.0005 $0.0006 +0~0008 +0~0011
0.9623 0.9540 0.9428 0.9273 0.9048
0.9608 0.9519 0.9397 0.9224 0.8969
+0.0015 10.0021 +0.0031 +0.0049 +0.0079
0.9626 0.9542
0.8712 0.8187 0.7351 0.6065 0.4493
0.8580 0.7977 0.7084 0.5972 0.4978
+0.0132 $0.0210 +0.0267 +0.0093 -0.0485
0.8690
0.8185 0.7383 0.6039 0.4476
+0.0022 +0.0002 -0.0032 +0.0026 +0.0017
0.3679 0.4493
0.6065 0.7351 0.8187
0.4570 0.4978 0.5972 0.7084 0.7977
-0.0891 -0.0485 +0.0093 +0.0267 +0.0210
0.3722 0.4476 0.6039 0.7383 0.8185
-0.0043 _tO~OO17 +0.0026 -0.0032 +0~0002
0.8712 0.9048 0.9273 0.9428 0.9540
0.8580 0.8969 0.9224 0.9397 0.9519
+0.0132 +0.0079 to.0049 +0.0031 +0.0021
0.8690
0.9061 0.9278 0.9418 0.9542
0.9623 0.9685 0.9733
0.9771 0.9802
0.9608 0.9674 0.9725 0.9765 0.9797
+0.0015 +0.0011 +0.0008 +0.0006 +0~0005
0.9626 0.9683 0.9733 0.9770 0.9803
0.9827 0.9847
0.9819 0.9827
+0.0008 +0.0020
0.9828 0.9845
0.9827 0.9802
+0~0008
0.9418 0.9278 0.9061
-0~0001 +0.0002
little more charitable than we expect to encounter using real experimental data. As shown in Table 2, the maximum difference between i and t is -0.0891 and occurs at the absorption maximum. In this case, four iterations reduce the maximum error to -0*0105, and approximately twenty minutes of computation reduces the largest error to the order of 0.001 which is within typical experimental accuracy. Fig. 2 shows the results of a similar calculation in which t was obtained numerically at
Methods for the correction of apparent band shapes
1637
,-
0.06
-2-w 0.04 I-
0.02
C, .
-0.02
I-
-0.04
Lorentz 1
-0.06
A
bond
slit;
A’; = 0.75
i,,;
k 0.3679
t opprox. -0.06
\r I
-0.10 -15
-IO
I -5
I
I
0
5
I
I IO
15
Fig. 1. Finite slit width error initially (i - t) and after twenty-six of iterative correction (i - t,,) (t approximate)
cycles
0.04 L--h
0.02
0
i-f
w
-0.06
5
IO
IC
Fig. 2. Finite slit width error initially (i - t) and after forty-nine of iterative correction (i - td9) (t “exact”)
cycles
A. L. KHIDIR and J. C. DECIUS
1638
high resolution (s = 30) ; the back calculation was done at s = 3, and it is seen that a maximum initial error is reduced from -0.097 to -0*016 after 49 iterations. A third example given in Fig. 3 illustrates a case in which s = 3 with the parameter a chanced so that the minimum transmission is 0.1; after seventeen iterations, the largest error is reduced from -0.074 to +0*0018.
I
-15
I
-10
I
-5
I 0 v- Y0
I 5
I IO
Fig. 3. Finite slit width error initially (i - t) and after seventeen cycles of iterative correction (i - tI,)
There are some tricks which one can play to accelerate convergence, but we prefer not to discuss them at this time. It is very interesting to examine the periodicity in the residual error. For the case s = 3, we expect slowest convergence for the Fourier component Ic w N/s = Il. The plots of the errors, shown in Figs. 1, 2 and 3 confirm this expectation, and a similar effect was observed for the case s = 6 where k RS 5. The success of the approximation calculation therefore is clearly due to the smallness of the corresponding Fourier components in the t function. Gaussian slit ; Lorentzian band One example of this type will now be described. The (unnormalized) slit function is taken as exp [ -2*77(a~’ - v)~/s~] to retain s as the definition of the width of the slit function at half height. An initial maximum error of -0.0923 is reduced to 0.0029 rather slowly (after 200 iterations) ; here again the ratio of slit width (s = 3) to band width (Av,,, = 4) is O-75.
Methods
for the correction
of apparent
band shapes
1639
An experimental example A final illustration is presented based on data recorded with a Beckman IR-7 prism-grating spectrometer for the 771 cm-l band of m-chlorotoluene in CS, solution. The absorption was recorded at two nominal slit widths (i) at s = 4.5 cm-l, for which Avl,aa = 7.5 cm-l and (ii) at s = 0.54 cm-l, for which Avflz = 5.1 cm-i. The latter recording was assumed to be a reasonably accurate representation of the true transmission curve. Using the same procedure as was employed above, the --o-o-_ 60
$ a E E F E
*,.-.-O-CO-H
-
60-
40
.
I
.
I 770
780 Wovenumber,
I 760
I 750
cd
Fig. 4. Infrared absorption of m-chlorotoluene in CS, --observed spectrum, nominal slit 4.5 cm-1 ---observed spectrum, nominal slit 0.5 cm-l * . . calculated spectrum based on wide slit scan
data obtained from the curve traced at s = 4.5 cm-l were subjected to iteration (N = 32) on the assumption of a triangular slit shape, one condition of the experiment being that the calculator was not shown the “true,” i.e. high-resolution curve until he had concluded that the iterative process had converged. The integrated absorption coefficient, whose calculation was incorporated in the program for the digital computer, was used as a criterion for convergence. Although a relatively high error (3 per cent) persists at the band center (see Fig. 4), there is enough improvement to warrant the assumption that part of the final error may be attributed to the oversimplification in taking a triangular slit shape, and perhaps in adopting the nominal value of the slit width. CONCLUSION Provided the ratio of slit width to true band width does not appreciably exceed the maximum ratio reported in RAMSAY’S [Z] tables, it appears feasible to obtain improved estimates of the true shapes of absorption bands without making assumptions concerning their functional shape. Thus only the spectrometer function remains as an unknown, which presumably can be reasonably well defined by the analysis of a number of examples. The difficulties inherent in this sort of calculation have been explicitely stated in terms of the near singularity of the spectrometer function. 11
Numerical methods for the correction of apparent band shapes due to finite slit width A. L. KHIDIR
and J. C. DECIUS*
Department of Chemistry, Oregon State University, Corvallis, Oregon (Received 7 &Iuy 1962) Abstract-The problem of transforming apparent spectroscopic band shapes to their true values is considered. The integral equations which express the influence of a spectrometer with finite slit width are given a matrix representation in which the problem is reduced to the inversion of Although an exact solution of the the matrix which represents the spectrometer transformation. inversion problem is readily found for a wide class of slit functions, it is shown to be impractical in the sense that a small experimental error in the observed apparent intensities leads to a much larger error in the true intensities calculated with the aid of the inverse spectrometer matrix. Nevertheless, it is demonstrated that a satisfactory approximate solution can be found by an iterative calculation which converges sufficiently rapidly provided the slit width is not too much greater than the band width, and that the band does not contain extensive fine structure. The method makes no assumptions concerning the true shape of the band and may, therefore, be used in the case of asymmetrical bands. Neither is it restricted to any particular shape of the slit function (e.g. triangular, trapezoidal, gaussian). Several numerical examples are used to illustrate the flexibility of the method and the speed of convergence to be expected in practice.
THE correction of errors in spectrometric data due to the finite resolving power of real spectrometers is of special importance in connection with infrared absorption measurements of the so-called absolute (integrated) intensities. The techniques introduced by WILSON and WELLS [l] are widely used in the investigation of gases. In the case of condensed phase samples, it is possible to use extrapolation procedures or to employ tabulated corrections to the maximum absorbance, the band width, or the integrated absorption coefficient as given by RAMSAY [2]. The latter procedure suffers from the disadvantage that it is premised on the two assumptions (i) that the true band shape is Lorentzian and (ii) that the spectrometer transmission function is triangular, while the former procedure makes it difficult to detect concentration effects on the integrated intensity arising from solute-solvent interactions. It is, therefore, of interest to consider whether more flexible methods can be devised to make such corrections. A very general solution of this problem has been given by HARDY and YOUNG [3], who, by the use of the theory of Fourier transforms, obtain expressions for the true intensity at any point in the spectrum as a series which involves the apparent intensity and its derivatives of all orders at that point; the properties of the spectrometer * Alfred P. Sloan Foundation
Fellow.
[l] E. B. WILSON, Jr. and A. J. WELLS, J. Chem. Phys. 14,578 (1946). [2] D. A. RAMSAY, J. Am. Chem. Sot. 74, 72 (1952). [3] A. C. HARDY and F. M. YOUNG, J. Opt. Sot. Am. 39, 265 (1949). 1629
A. L. KRIDIR and J. C. DECIUS
1630
are expressed through a set of coefficients for the derivatives of the apparent absorption curve. Despite the generality of HARDY and YOUNG’S solution it does not seem to have been widely used in practice. The same fundamental difficulties which impede the application of their technique will limit the use of the method which is described in the present paper, which may, however, have the advantage of exposing these difficulties in a more explicit fashion. General discussions of these difficulties have been given by RAUTIAN (4). The mathematical technique to be described below fulfils the aim of flexibility, and the tedium of the computation required can be mitigated by the application of There are, nevertheless, some serious difficulties an automatic digital computer. involved, even when automatic computers are applied to the problem. Although an exact formal solution can be found for a wide class of spectrometer transforms, this solution proves to be without direct value, since it is excessively sensitive to error in the experimental data. The exact solution is nevertheless of value in suggesting convergent iterative methods and the conditions under which they may be usefully applied. The relation of these methods to the basic problem of resolution is discussed, and some plausible examples are worked out numerically. The discussion is restricted to the case of bands which are practically isolated in the spectrum. GENERAL SOLUTIONOF THE SPECTROMETERTRANSFORMATION Matrix approximation The influence of the spectrometer on the true incident energy spectrum, expressed in the form of a transmitted energy spectrum, t(d) as a function nominal spectrometer frequency v’ where t(v’) = p3(v’,v)
i (v) dv
i(v), is of the
(1)
in which the spectrometer transformation is expressed by S(v’, v). A formal solution of our problem is obtained if an inverse transformation function, X-l(v”, v’), can be found which satisfies js-l(v”, v’)S(v’, Y)dV’ = B(vN,v) (2) in which B(v”, v) is a Dirac function, since then i(v”) = js-l(v”,
v’)t(v’) dv
(3)
In practice of course t(d) is not given analytically, but is known graphically and can be tabulated from the observed data for come convenient set of discrete values of v’ . This suggests that (l), (2), and (3) should be replaced by the system of matrix-vector equations : t = Si (I’)
S-IS = E i = S-It [4] S. G. RAUTIAN,Sov. Phys. Uspekhi 66, 245 (1958).
(2’) (3’)
Methods for the correction of apparent band shapes
1631
in which the components of the vectors i and t are simply the values of the incident and transmitted spectra at a set of discrete frequencies: %?I= ;(Y,);
(4)
t, = t(d,)
In such a representation a spectrometer of unlimited resolving power would have a matrix S which would be a unit matrix. A spectrometer of limited resolving power would be represented by a matrix with off-diagonal terms. Cyclic form of the spectrometer matrix, S It is convenient in several respects to make the assumption that S is cyclic in the sense that each row is a cyclic permutation of the first row, i.e. that S has the form
S=
‘abed
....
babe
.......
cbab
........
.
.
.
.
bed....
..dcb
-
dc d (5)
. dcba
For one thing, such a form for S ensures that if i, is constant, ization of S is achieved by the requirement &
= 1
so is t,.
Normal-
(6)
In practical examples to be given below, a few of the elements a, b, c, d, etc. near the diagonal will be taken as different from zero, with many vanishing elements farther from the diagonal. If we suppose in particular that there are four nonvanishing elements, a, b, c, d, in the notation of (5), with N, say of the order of 20 to 50, we must justify the inclusion of the non-vanishing elements d, c, b in the upper right-hand corner of S. The first row of S will ordinarily be applied to calculate the intensity transmitted by the spectrometer at a frequency vl’, i.e. at the low-frequency edge of the band. The elements in the upper right-hand corner of S represents a fictitious transmission of energy arising from incident frequency components $N_2, iN_l, i, at the high-frequency edge of the band. But provided the interval of frequency and N are chosen properly, the energy in the region of iN_z, i,.,_i, i,,, will just compensate for the energy which is lost in consequence of the finite-frequency range. In other words, at spectrometer setting vl’, the contribution to the transmitted intensity by real frequencies less than v1 = vl’ 1s ’ compensated, provided i(v_J = i(vN), i(v_J = i(~~__~),etc. This condition will be satisfactorily fulfilled provided the incident spectral intensity does not change rapidly over the width of the absorption band. Later we shall introduce an iterative method of solution which does not require that S be cyclic, but it has been found convenient to retain the cyclic form both to demonstrate the convergence of the iterative method and to conserve memory space in calculations performed with an automatic digital computer.
A. L. KHIDIR and J. C. DECIUS
1632
Inversion of the spectrometer matrix We are now ready to reap the main advantages consequent from the assumed form of S. Our fundamental problem is simply the calculation of the inverse of the S matrix. This will be done by appealing to the diagonal form of 5. Any cyclic matrix of the assumed form may be diagonalized by a similarity transformation with a unitary matrix, U,
D = USU-l
= USUt
where the dagger indicates transpose conjugate. given by
(‘1
In particular, the elements of U are
in which N is the number of rows and colums in S, and k and m will be assigned the ranges 0 to N-l and 1 to N respectively. The elements of D are D,,,
= d,,, D(k) = a + 2b cos (27&/N) + 2c cos (2~2Tc/N)
(9)
+ 2d cos (2~*31c/N) + etc. Then since
(D-l),,
= &,/D(k)
= 6,,.D-l(k)
and
(10)
S-l = Uf D-l U
(11)
the elements of S-l are easily calculated and have the form I s&l = ; 1 + 2 cn_d”eD-l(k) cos 2nk(p - q)/N k=l C where N is odd, or A,-, = ;
1 + 2
(N--2)/2 2 D-l(k) cos 2vk(p k=l
1
q)/N + (-l)“-*D-l
(12)
(N/2)
1
(13)
where N is even. For large values of N the calculation required in (12) or (13) may be carried out very simply on a digital computer ; we have, in fact, programmed such a computation for the Alwac-III of the Department of Mathematics at Oregon State University.* Resolving power an& singularity of the S matrix One may well raise the question whether the general solution of the transform problem is capable of regaining the information presumed lost due to the finite resolving power of the spectrometer. Put differently, if the inversion of S expressed by (12) or (13) remains valid as N --f co, which is to say that Av = v,+~ - v, -+ 0, then in principle the true energy incident upon the spectrometer, i(v), could be calculated from the energy transmitted by the spectrometer, t(v’), with unlimited resolution, provided, of course, that one knows the spectrometer transform function, S(v’, v). * However,
see Table 1.
Methods for the correction of apparent band shapes
1633
But suppose that S is singular, i.e. that for some 4 = S&/N, D(4) = D(k) = a + 2b cos + + 2c cos 2+ + etc. = 0. Then, of course, 8-l will not exist, and our method will fail. It is easy to see how this will in fact occur for some reasonable spectrometer functions. Consider, for example, the triangular silt function. Designate by s the half width of the slit in units of Av and for simplicity take s as an integer. The constants a, b, c, cl, etc. are then given by a = l/s = s/s2 b = (s -
l)/s2
c =
(s -
S)/s2
d = (s -
3)/s2
(14)
etc. Then S will be singular if s20(+)
= s + 2(s -
1) cos 95 +
2(s -
2) cos 2+ + etc. = 0
(15)
This will in fact be the case if k = N/s
(16)
i.e. if N is divisible by s. This will surely occur in the limit as N + co. Even in the case of a finite approximation where N and s are chosen so that N is not divisible by s, one may expect to encounter near singularities when both N and s are reasonably large, and therefore great caution must be exercised in performing the calculation of (12) or (13) to avoid undue loss of arithmetic significance. As an example, we have calculated S-1 for the triangular slit with N = 41 and s = 3. The eigenvalues of S do not exceed unity, but for k = 14, D(k) = 0*00084365 so that D-l(k) is known to only about four significant places even though the calculation is done originally with 7 or 8 places. The values of D(k) and of SF, are given in Table 1. This, however, is not the most serious difficulty. In practice, the exact* solution turns out to be essentially useless for the following reasons. The elements of S-1 are all large compared with unity, at least for practical values of N such as N > 20. Under these circumstances, the errors in the calculated values of i become large compared with the errors in the observed values of the t. Specifically, if the errors in the t are approximately O-1 per cent the errors in the calculated values of i will be approximately lo-15 per cent. Put somewhat differently, the calculated values of the i are linear combinations of all the t, of a form which involves the approximate cancellation of many large numbers of opposite sign. Thus an unbearable burden is put upon the experimentalist, since no spectroscopist-at any rate none working in the infrared region with conventional equipment-will relish the task of recording an absorption spectrum with the O-01 per cent intensity precision which would be required in the above example to yield a final result good to 1 per cent. We may now consider a different approach, which still utilizes our knowledge of the diagonalization of S. Consider the difference i -
t = (E - S)i
(17)
* Here and elsewhere we mean by “exact solution,” the exact inversion of the assumed S matrix, which is of course itself only an approximation,
of finite resolving power, to (1) and (3).
1634
A. L. KHIDIR and J. C. DEGIUS
and, assuming the difference to be small, suppose that, starting with the initial approximation i. = t on the right-hand side of (17), the iterative process defined by
i, - t = (E - S)i,_,
(18)
Table 1. Eigenvalues of S and elements of S-l for triangular slit with N = 41 and a = 3 Eigenvalues
Inverse matrix, S-l
k
D(k) = D( -ii)
P-9
SZ
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1~00000000 0.98443472 0.93882698 0.86634570 0.77197224 0.66207677 054388195 0.42486048 0.31212103 0.21183997 0.12878662 0.06598623 0.02454835 0.00367495 0.00084365 0.01214789 0.03275841 0.05746238 0.08122593 0.09972774 0~10981087
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
s3* -43 -34 74 -43 -25 65 -43 -16 56 -43 7 47 -43 2 38 -43 11 29 -43 20
* Candour compels us to admit that the SG values were calculated numerically using (12) ; the fact that S,-,l are integers was noted empirically and verified by checking that SPIS is indeed the unit matrix.
will converge. imation is
In fact the difference
between
the n-th and (n -
l)-th
approx-
i, - in_, = (E - S)“t so that if t is expanded in the eigenvectors u tk) of E - 8, which are identical with those of S and whose components are defined by (8), one finds i, -
i,_l = (E - S)nz tkutk) = 2 t,jl,~(~) k
where 1, is an eigenvalue of E -
(19)
S and is therefore given by Ak = 1 -
D(k)
(20)
(see equation 9). Thus provided D(k) f 0, the difference between successive approximations to i approaches zero as n goes to infinity. Clearly the difficulty in this
Methods for the correction of apparent band shapes
1635
process in the slowness of convergence associated with very small values of D(k). Moreover, if it proves necessary to carry the iteration in (18) to a large value of n, it can easily be shown that the result will be equivalent to the use of S-l which we have already shown to be impractical. However, if t, in (19) is negligible for values of k such that 2, is very close to unity, then the iteration may converge with practical accuracy for sufficiently small values of n thereby yielding a useful result. We have found this to be true for several reasonable examples. These examples have two characteristics : (1) the band is reasonably smooth and (2) the slit is not more than a few times the true band width. Despite the success of our method in these cases, to be described in greater detail below, one can anticipate difficulties if one pushes the procedure too far, for instance in attempting to resolve extensive fine structure, since in such cases the t, for higher-frequency Fourier components would be larger. ITERATWE INVERSION OF SOME EXAMPLES Although zian, i.e.
we have chosen as some sample bands i(y) which are either Lorenti(y) = exp --E(V) cc(Y) =
(21)
a (Y -
(22)
YO)2+ b2
or Gaussian, U(Y) = a exp [-(Y
-
aQ2/b2]
(23)
it must be emphasized that our method is in no sense restricted to these cases and indeed was developed specifically to avoid this restriction. Also, we have used as a spectrometer function the two forms of (i) a triangular and (ii) a Gaussian transform, but this restriction is also not essential. In this connection it is interesting to note that KOSTKOWSKIand BASS (5) in the course of a related investigation have demonstrated experimentally the adequacy of a Gaussian spectrometer function in the near infrared. The program for the iterative calculation is based essentially upon (18), except that in order to accelerate convergence, those components of i, which have already been calculated during the n-th cycle are used in place of the corresponding components of i,_, on the right hand side of (18). Each cycle is initiated at a wing of the band, and the program has not taken advantage of the symmetry of the bands chosen for examples, since we desired explicitly to treat asymmetric bands. Experience has shown that a desk calculator would be adequate for the slightly distorted functions t which arise when the silt width does not exceed one-half of the apparent band width. Triangulur slit ; Lorentzian band Table 2 illustrates the results of the numerical calculations along with Fig. 1. The true band shape was based on the following parameters: a = 4, b = 2 in (21). The values of i, and t, are given for p and q running from 1 to 32 ; p = q = 16 corresponds to the band centre, which amounts to 36.79 per cent transmission. [5] H. J. KOSTEOWSEIand A. M. BASS, J. Opt. Xoc. Am. 46, 1060 (1956).
A. L. KHIDIR and J. C. DECIUS
1636
The t, were calculated using a triangular S function with s = 3. Thus the apparent band shape itself involves the approximation due to taking discrete values of the data, and our numerical test of the inversion in which we seek to back calculate i is a Table2. Lorentzian band; transform table a = 4,6 = 2,Av = 1,s = 3 i
t
i-t
1”26
L - t2f3
0.9828 0.9803 0.9770 0.9733 0.9683
0.9771 0.9733 0.9685
0.9819 0.9797 0.9765 0.9725 0.9674
+0.0005 $0.0006 +0~0008 +0~0011
0.9623 0.9540 0.9428 0.9273 0.9048
0.9608 0.9519 0.9397 0.9224 0.8969
+0.0015 10.0021 +0.0031 +0.0049 +0.0079
0.9626 0.9542
0.8712 0.8187 0.7351 0.6065 0.4493
0.8580 0.7977 0.7084 0.5972 0.4978
+0.0132 $0.0210 +0.0267 +0.0093 -0.0485
0.8690
0.8185 0.7383 0.6039 0.4476
+0.0022 +0.0002 -0.0032 +0.0026 +0.0017
0.3679 0.4493
0.6065 0.7351 0.8187
0.4570 0.4978 0.5972 0.7084 0.7977
-0.0891 -0.0485 +0.0093 +0.0267 +0.0210
0.3722 0.4476 0.6039 0.7383 0.8185
-0.0043 _tO~OO17 +0.0026 -0.0032 +0~0002
0.8712 0.9048 0.9273 0.9428 0.9540
0.8580 0.8969 0.9224 0.9397 0.9519
+0.0132 +0.0079 to.0049 +0.0031 +0.0021
0.8690
0.9061 0.9278 0.9418 0.9542
0.9623 0.9685 0.9733
0.9771 0.9802
0.9608 0.9674 0.9725 0.9765 0.9797
+0.0015 +0.0011 +0.0008 +0.0006 +0~0005
0.9626 0.9683 0.9733 0.9770 0.9803
0.9827 0.9847
0.9819 0.9827
+0.0008 +0.0020
0.9828 0.9845
0.9827 0.9802
+0~0008
0.9418 0.9278 0.9061
-0~0001 +0.0002
little more charitable than we expect to encounter using real experimental data. As shown in Table 2, the maximum difference between i and t is -0.0891 and occurs at the absorption maximum. In this case, four iterations reduce the maximum error to -0*0105, and approximately twenty minutes of computation reduces the largest error to the order of 0.001 which is within typical experimental accuracy. Fig. 2 shows the results of a similar calculation in which t was obtained numerically at
Methods for the correction of apparent band shapes
1637
,-
0.06
-2-w 0.04 I-
0.02
C, .
-0.02
I-
-0.04
Lorentz 1
-0.06
A
bond
slit;
A’; = 0.75
i,,;
k 0.3679
t opprox. -0.06
\r I
-0.10 -15
-IO
I -5
I
I
0
5
I
I IO
15
Fig. 1. Finite slit width error initially (i - t) and after twenty-six of iterative correction (i - t,,) (t approximate)
cycles
0.04 L--h
0.02
0
i-f
w
-0.06
5
IO
IC
Fig. 2. Finite slit width error initially (i - t) and after forty-nine of iterative correction (i - td9) (t “exact”)
cycles
A. L. KHIDIR and J. C. DECIUS
1638
high resolution (s = 30) ; the back calculation was done at s = 3, and it is seen that a maximum initial error is reduced from -0.097 to -0*016 after 49 iterations. A third example given in Fig. 3 illustrates a case in which s = 3 with the parameter a chanced so that the minimum transmission is 0.1; after seventeen iterations, the largest error is reduced from -0.074 to +0*0018.
I
-15
I
-10
I
-5
I 0 v- Y0
I 5
I IO
Fig. 3. Finite slit width error initially (i - t) and after seventeen cycles of iterative correction (i - tI,)
There are some tricks which one can play to accelerate convergence, but we prefer not to discuss them at this time. It is very interesting to examine the periodicity in the residual error. For the case s = 3, we expect slowest convergence for the Fourier component Ic w N/s = Il. The plots of the errors, shown in Figs. 1, 2 and 3 confirm this expectation, and a similar effect was observed for the case s = 6 where k RS 5. The success of the approximation calculation therefore is clearly due to the smallness of the corresponding Fourier components in the t function. Gaussian slit ; Lorentzian band One example of this type will now be described. The (unnormalized) slit function is taken as exp [ -2*77(a~’ - v)~/s~] to retain s as the definition of the width of the slit function at half height. An initial maximum error of -0.0923 is reduced to 0.0029 rather slowly (after 200 iterations) ; here again the ratio of slit width (s = 3) to band width (Av,,, = 4) is O-75.
Methods
for the correction
of apparent
band shapes
1639
An experimental example A final illustration is presented based on data recorded with a Beckman IR-7 prism-grating spectrometer for the 771 cm-l band of m-chlorotoluene in CS, solution. The absorption was recorded at two nominal slit widths (i) at s = 4.5 cm-l, for which Avl,aa = 7.5 cm-l and (ii) at s = 0.54 cm-l, for which Avflz = 5.1 cm-i. The latter recording was assumed to be a reasonably accurate representation of the true transmission curve. Using the same procedure as was employed above, the --o-o-_ 60
$ a E E F E
*,.-.-O-CO-H
-
60-
40
.
I
.
I 770
780 Wovenumber,
I 760
I 750
cd
Fig. 4. Infrared absorption of m-chlorotoluene in CS, --observed spectrum, nominal slit 4.5 cm-1 ---observed spectrum, nominal slit 0.5 cm-l * . . calculated spectrum based on wide slit scan
data obtained from the curve traced at s = 4.5 cm-l were subjected to iteration (N = 32) on the assumption of a triangular slit shape, one condition of the experiment being that the calculator was not shown the “true,” i.e. high-resolution curve until he had concluded that the iterative process had converged. The integrated absorption coefficient, whose calculation was incorporated in the program for the digital computer, was used as a criterion for convergence. Although a relatively high error (3 per cent) persists at the band center (see Fig. 4), there is enough improvement to warrant the assumption that part of the final error may be attributed to the oversimplification in taking a triangular slit shape, and perhaps in adopting the nominal value of the slit width. CONCLUSION Provided the ratio of slit width to true band width does not appreciably exceed the maximum ratio reported in RAMSAY’S [Z] tables, it appears feasible to obtain improved estimates of the true shapes of absorption bands without making assumptions concerning their functional shape. Thus only the spectrometer function remains as an unknown, which presumably can be reasonably well defined by the analysis of a number of examples. The difficulties inherent in this sort of calculation have been explicitely stated in terms of the near singularity of the spectrometer function. 11