898
Anal. Chem. 1989, 61, 898-904
Improved Discrimination of Noisy Spectra in Correlation-Based Spectral Analysis Douglas P. Hoffmann, Andrew Proctor, and David M. Hercules* Department of Chemistry, University of Pittsburgh, Pittsburgh, Pennsylvania 15260
Spectral comparlsons udng cross-correlation techniques generally rely on a single parameter, 7 0 , whlch is the magnitude of the correlation function at zero displacement. However In the presence of noise T~ decreases as the slgnai to noise (S/N) decreases. This means that noisy data may be unduly rejected in any spectral searching procedure. By autoconvolutingthe cross-correlatlon pattern, a new parameter, T ~ can~ be obtained ~ , that contains information about the symmetry of the correlation function. The symmetry of the correlatlon function represents additional evidence of spectral slmliarlty, slnce an autocorrelation function is totally symmetrlc. Use of the dual parameters, T~ and 7syM,gives much better dlscrlminatlon against nolsy data, effectively allowing lower T~ values to be accepted. Any standard or library spectrum can be callbrated over a wide S/N range to ~ ~ extended . gauge effective cutoff nmRs for T~ and ' T ~ This method Is examined by using Mossbauer data, which can often be e x t m l y noisy because of a combination of a weak source and samples with low concentrations.
1. INTRODUCTION Spectroscopic identification of an unknown sample is an important analytical application, but the information that can be extracted from a spectrum ultimately will be limited by spectral quality. There are many factors that impose limitations on the ability to collect high signal-to-noise (S/N) data. Such limitations include instrumental factors, sample purity, the time frame of the measurement, and sample concentration. One technique in which high S / N requires a relatively long collection period (hours), especially for low concentration samples, is Mijssbauer spectroscopy. Mossbauer spectroscopy uses a radioactive isotopic source (e.g. 57C0)to excite nuclear transitions in a corresponding isotope (e.g. 57Fe). For a given experimental setup the quality of Mossbauer data is related to the activity of the radioisotope source, the collection time, and the concentration of the sample. With respect to the latter point the isotope that is being probed may not be the major naturally occurring isotope for that element. In the case of iron, 57Fe has a natural abundance of -2.1%. Unless the sample is specifically enriched in "Fe (not practical or possible for most measurements), the collection time for spectra from low weight percent iron samples (4wt %) may be extremely long (>lo h). It is therefore necessary that a method of data analysis be applied that can supply accurate information about the sample with a minimum amount of interference from noise, which perturbs the true spectral shape. In this paper we propose an extension to correlation-based methods of spectral comparison described by other workers (1-15). In this extended method the effect of noise in the experimental data can be overcome by considering the overall shape of the correlation function derived from correlating the noisy experimental spectrum with a standard or library spectrum having very little noise. A measure of the symmetry of the correlation function by examination of its autoconvolution function is the basis for the proposed method. The
application of this method in spectroscopy is universal, but Mossbauer data will be used to illustrate the usefulness and limitations of the suggested approach. Calibration of any standard spectrum using noisy versions of itself is a means by which suitable cutoff limits can be set for the required parameters. In the following section we shall outline the mathematical basis and spectroscopic implementation of correlation and convolution analysis. An understanding of the similarities and differences between correlation and convolution is a prerequisite for full appreciation of the data analysis method proposed. This is most easily achieved by consideration of the mathematical manipulations that are carried out in Fourier space (16, 17).
2. EXPERIMENTAL SECTION Transmission Mossbauer spectra were obtained with a Nuclear Scientific and Engineering Corp. Mijssbauer effect spectrometer (Model AM-l), equipped with a 57C0source. The source was calibrated with sodium nitroprusside. The multichannel analyzer was interfaced to an Apple IIe microcomputer for data collection. The data could then be readily transferred to a VAX main frame cluster or IBM PC/AT and the appropriate data analysis performed. 3. THEORY Discrete Correlation and Convolution. For two discrete periodic functions a(k) and b(k) the correlation of a(k) by b(k) is defined as (16) N-1
ccorr(7)
= C a ( k ) b(7 + k ) k=O
(la)
and the convolution of a ( k ) by b ( k ) is defined as N- 1
c,,,(7)
= C a ( k ) b(7 - k) k=O
(1b)
where 7 has incremental values from 0 to N - 1,and N is the number of data points in one period of the discrete function. Equations l a and 1b are very similar; however the seemingly minor difference in the sign of k in function b is of great importance. The physical significance of the above equations cannot easily be grasped, and because a full appreciation of correlation and convolution is extremely important to the remaining discussion, it is worthwhile to describe the equations graphically. Figure 1illustrates the essential characteristics of eq l a and lb. For purposes of this illustration a ( k ) and b ( k ) (Figure la,b) are different asymmetric functions with no implied spectroscopic significance. For correlation, eq la, the function a(k) (Figure la) is multiplied by the function b(7 + k ) (Figure IC);b(7 + k ) is simply b(k) shifted along the k axis by an amount 7. For convolution, eq Ib, the function a ( k ) (Figure la) is multiplied by the function b(7 - k ) (Figure le); b(7 - k ) is b(-k) (Figure IC)shifted along the k axis by an amount 7. Note the difference between correlation and convolution is that for convolution the b(k) function is folded about the ordinate to yield b(-k) prior to the shift by 7. The total correlation or convolution of a ( k ) by b ( k ) is obtained by evaluating the sum of the products at all values of 7. Figure 2 illustrates several different correlation/convolution combinations of the functions a ( k ) and b ( k ) from Figure 1.
0003-2700/89/0361-089S$01.50/00 1989 American Chemical Society
ANALYTICAL CHEMISTRY, VOL. 61, NO. 8, APRIL 15, 1989 890
n
k=O
N-1
I b(T+k)
7 = 0
’
7 = (NB)
Flgure 3. Representation of discrete correlation of spectrum a(k)by
+
spectrum b ( r k ) illustrating the assumed periodic nature of the spectra. The spectrum b ( r + k ) Is shown for several different values of 7 .
+T
Flgwe 1. Two dissimilar asymmetric spectral profiles: (a) a ( k ) and (b) b ( k ) , (c) correlation function b(7 k ) shown with a posRive displacement of 7 , (d) b(-k), the same as b ( k ) but folded about the ordlnate, (e) convoluting function b(7 - k ) shown with a positive dis-
+
placement of
7.
b corr a
(4
a corr a
(f)
a conv b
(g)
b conv a
fi
(h)
-
0
r
a conv a
+
Figure 2. Correlation and convolution spectral profiles obtained for (a) a ( k )and (b) b ( k ) , (c) correlation of a ( k ) by b ( k ) , (d) correlation of b ( k ) by a (k),(e) autocorrelation of a ( k ) , (f) convolution of a ( k ) by b(k), (9)convolution of b ( k ) by a ( k ) ,and (h) autoconvolution of a ( k ) .
The full explanations for the similarities and differences can be appreciated by examining the Fourier space manipulations (vide infra) but it is important to point out some significant properties at this stage. Figure 2c,d shows that the order of the correlation [a(k)by b(k) or b(k) by a & ) ] of asymmetric functions affects the correlation function. For the convolution of [a(k)by b(k) or b(k) by a(k)]it can be seen from parts f and g of Figure 2, respectively, that the order of convolution is unimportant. For a(k)and b(k) the shapes of the correlation and convolution functions are different; however this is not a general result and arises because a(k) and b(k) are both
asymmetric and different. The autocorrelation or autoconvolution of a function is the correlation or convolution of that function with itself. Thus parts e and h of Figure 2 illustrate the autocorrelation and autoconvolution of function a(k), respectively. The important point to note is that the autocorrelation of a(k) is a totally symmetric function whose maximum value occurs at a displacement of 7 = 0. The value of any correlation function at zero displacement (7 = 0) is termed T~ In contrast to autocorrelation the autoconvolution of a&), Figure 2h, is asymmetric and furthermore has a zero or near zero value at 7 = 0. The value at this point in the autoconvolution function we shall refer to as ~ s y as ~ it , relates to the symmetry of the original function (vide infra). Figure 3 illustrates discrete correlation applied to spectra a(k) and b(k) (k = 0, N - 1). It is of tremendous importance to recognize that in carrying out such analysis the measured digital spectrum is considered as a single period in an infinite periodic function. Figure 3 shows the sliding of b ( k ) with respect to a(k) for several different displacements of 7. Instead of just showing a single period (from k = 0 to N - 1)of the various b(7 + k) positions, two additional periods are shown to give insight into how the function “wraps around” as it is displaced through its entire T range. These extra periods also illustrate that the two extreme ends of the spectrum (k = 0 and N - 1) should merge such that there are no sharp discontinuities which would introduce spurious end effects. The final important point of Figure 3 is that although 7 varies from 0 to N - 1, the effective 7 value will be both positive and negative. Beginning at 7 = 0 displacement (Figure 3), an increase in 7 up to 7 = N / 2 shifts b(k) 7 the right. Further increases in 7 > N / 2 continue to shift b(k) to the right, but this results in 7 moving out of the 12 = 0, N - 1 window. It is evident from the bottom diagram (Figure 3) that these positive shifts of 7 > N / 2 are equivalent to negative shifts of 7’ = ( 7 - N). Evaluation of the so or rsYM values can easily be determined by multiplying equivalent points in the two relevant spectra together at 7 = 0 (zero displacement) and summing all the individual products. However this does not produce the total correlation or convolution function, merely one point out of a possible N points. Production of the whole function clearly requires W multiplications. This computationally undesirable situation has many similarities to evaluation of the discrete Fourier transform (16, 18). The problem can be overcome in this latter case by application of the fast Fourier algorithm (19) which substantially reduces the number of required multiplications and corresponding computation time. This has significance here because Fourier transformation of eq l a and l b produce the
900
ANALYTICAL CHEMISTRY, VOL. 61, NO. 8, APRIL 15, 1989
following useful results (16):
Ccorr(n)= A h ) B*(n) (24 Cconv(n)= A b ) B(n) (2b) where n = 0 to N - 1and C,(n), C,,(n), A(n),and B(n) are the Fourier transforms (FT)of C ~ J T ) c,,(T), , a(k),and b(k), respectively. B*(n)is the complex conjugate of B(n). Thus the rather complicated sliding maneuvers carried out in the pseudo time domain can be carried out more easily in the pseudo frequency or Fourier domain simply by multiplication. The correlation or convolution function can then be retrieved by inverse Fourier transformation. Because of the FFT, the forward and inverse Fourier transformation route represents a significant savings in computation time. In addition, as with filtering (20), it is far easier to appreciate some of the properties of the correlation and convolution process by examining them in Fourier space. Fourier Space Multiplications. The Fourier transform X is a complex quantity equal to R, + iI,, where R, and I, are the real and imaginary parts of X and i is -ll/z. Expansion of eq 2a and 2b leads to
CCORR= @ARB + IAIB) i(IARB - IBRA) ~ICORR = RCORR + CCONV= (BARB- IAIB) + ~ ( I A F+B IBRA) + IICONV = RCONV
(34
(3b) From analysis of these equations the similarities and differences of correlation and convolution become apparent. (Note that both the correlation and convolution functions are the real part of the inverse transform.) Of particular interest in the discussion is the condition where A = B , i.e. autocorrelation and autoconvolution. Autocorrelation and Autoconvolution. For autocorrelation or autoconvolution of any function X eq 3a and 3b reduce to ACCoRR = R:
+ 1,' = 1x1'
(44
ACcoNv = (R: - I X 2 )+ i2R,I,
1x1
where is the magnitude of X. Clearly the autocorrelation transform is totally real and thus the autocorrelation function will always be an even or totally symmetric function (16). This is not true for autoconvolution in general since the transform is a complex quantity. However if x ( k ) is an even or symmetric function, then I, will be zero (17)and so the autoconvolution transform will be totally real resulting in an even autoconvolution function. In this case autocorrelation and autoconvolution are equivalent. Because of the reflection prior to displacement in convolution analysis, it seems clear that the autoconvolution operation will provide information about the symmetry of the original function x ( k ) . Normalization, r0,and rSyM.The value of any correlation function a t T = 0 has been termed T~ while the value of the autoconvolution a t T = 0 is referred to as T S Y M . In order to quantify these values we must first normalize the data. It can be shown by application of the inverse discrete F T to eq 4a that T~ for an autocorrelation is simply the weighted sum of the real coefficients, lXlz 1N-1 TO = (Rx2+Ix2) (5) Nn=o The 1/N term arises as a factor in the inverse FFT. Therefore if both R, and I , are first multiplied by F , where
then clearly To(auto)will equal unity. This will be the max-
imum value that T~ can have and this will only arise when the correlation is an autocorrelation. Applying the same normalization factor F, to the real and imaginary coefficients of the autoconvolution function (eq 4b), then the normalized autoconvolution, ACV(T), is (16) ACV( T ) =
Therefore a t
T
= 0 with substitution from eq 6
If the correlation function is symmetric (even), then its T S ~ M value = 1.0 since I, = 0. Conversely if it is an antisymmetric (odd) function, then R, = 0 and T S ~ M= -1.0. For any function that lies between these two extremes, T S ~ Mwill correspondingly have a value between 1.0 and -1.0. Phase Correlation. Kawata et al. (21) have suggested that better discrimination could be achieved during correlationbased analysis by utilizing the phase content of the correlation function alone rather than both the phase and the magnitude. Since a complex quantity X is also equal to exp(iO,), where 8, is arctan (IJR,) (16),eq 2a can be written
1x1
C = IAllBl exp[i(OA - OB)]
(8)
Division by the magnitude of C, IC1 = IA(IB1,will leave = exp[i(6A - OB)] = COS (0, - OB) + i sin (8, - 6,) (9) Inverse transformation of this expression will produce the phase correlation function having a value at T = 0 to which we shall refer as T~ Normalization by eq 6 (here F2 = 1)will again set the value of T~ for autocorrelation to be unity. For autocorrelation, examination of eq 9 shows that CPk is unity a t all points since BA = OB. The FT of this function will be a delta function or spike a t T = 0. For correlation in general, the correlation function will be a spike at T = 0 whose height ~p is less than or equal to unity. The major problem with this model is noise. The effect of noise can be understood by considering that division by IAIJBIin the Fourier domain is equivalent to multiplication by its reciprocal which is a function that dramatically increases with frequency and so enhances the high frequencies over the lower frequencies. This results in considerably more noise in the back transformed data. A comparison can be made with deconvolution. If a ( k ) is to be Fourier deconvoluted by b ( k ) then in Fourier space, we simply divide A by B Cphese
Thus deconvolution of a ( k ) by b ( k ) is equivalent to a weighted correlation of a ( k ) by b ( k ) . In a similar way the phase correlation is very much like deconvolution but here by Fourier weighting is IAllBl rather the IBI2. Thus phase correlation is a totally deconvoluted spectrum that contains no information about the peak shape. Extended Method Outline. Figure 4 outlines the extended method pathway used for the correlation/convolution analysis. Correlation of the measured spectrum and a standard spectrum is first carried out to yield the T~ (and TP) value. The correlation function itself is then autoconvoluted to yield the T S Y M value. The advantage of this extension is that T S Y M gives information about the symmetry of the correlation function. In general the more symmetric the corre-
ANALYTICAL CHEMISTRY, VOL. 61, NO. 8, APRIL 15, 1989
901
where RND is a random normal deviate which is defined as (22)
ri is a random number in the range 0-1. Choosing N = 1 2 reduces eq l l b to N
RND = (CrJ- 6 i-1
Flgure 4. Extended method pathway. The steps used in the comparison of an experimental spectrum with a set of library spectra or the calibration of a library spectrum using noise. See text for expianation.
Table I. Relationship between Spectral Symmetry and T~ and T ~ Y MValuesd
spectrum A
spectrum B
C"
To
identical' EVEN
identical' EVEN
ODD
ODD ODD
EVEN EVEN EVEN
1.0 51.0 51.0
ODD ODD
0.0 0.0
