12670
J . Phys. Chem. 1993,97, 12670-12673
Wavelet Fast Fourier Transform (WFFT) Analysis of a Millivolt Signal for a Transient Oscillating Chemical Reaction Delmar N. S. Permand Canadian Government Laboratory Visiting Fellow, Emergencies Science Division, Environment Canada, River Road Environmental Technology Centre, 3439 River Road, Ottawa, Ontario, Canada K I A OH3
Heshel Teitelbaum Department of Chemistry, The University of Ottawa, Ottawa, Ontario, Canada KIN 6N5 Received: October 15, 1993"
A procedure is described for the W F F T (wavelet fast Fourier transform) analysis of a real millivolt signal that has been collected with a moderately fast data acquisition system, 40 million samples per second (msps). The output signal is noisy, and the frequency components are resolved by using a wavelet analysis followed by an FFT of selected parts of the wavelet output. After processing, the signal is then reconstructed (from selected components) resulting in a waveform that is almost identical to the original but with reduced noise, thus simulating signal averaging. The method is demonstrated on a signal from a single-shot experiment on a gas-phase oscillating chemical reaction occurring on the microsecond time scale.
Wavelet analysis has been used to analyze geological signals, speech and vision, and also irregular heartbeats.' It has also been used to detect a weak forcing term in an otherwise chaotic signal.* The present study is directed toward the application of wavelet analysis to real measured signals complete with noise as in ref 3. Mallat and Hwang4 have done work toward recovering a computer-generated signal to which they had added noise, and their work is used to explain some of the results obtained in this study. The Fourier transform (FT) and the fast Fourier transform (FFT) are known to be good methods for analyzing repetitive signal^,^,^ but these methods have some limitations when it comes to analyzing transient signals. There are ways to overcome some of these limitations.5v6 This work will not pursue those techniques, however, but will instead show, for the first time, how the wavelet analysis in conjunction with a FFT can analyze a real millivolt transient signal obtained from monitoring an oscillating chemical reaction. The chemical reaction and apparatus, on the one hand, and the laser schlieren detector, on the other hand, have been described in separate previous p~blications.~-9 The data collection system, for this study, and the modified schlieren detection system are described in detail in ref 3. A short summary is presented below: CC13F, which used to be a common refrigerant, was shockcompressed to produce a supersaturated gas (S < 2). Under these conditions the initial distribution of clusters, (CC13F),,, was perturbed, and a new distribution was formed via reactions of the form M
where * denotes an unstable species which can be stabilized upon collision with inert molecules, M, such as argon (Ar) or CC13F itself. The rate of homogeneous nucleation is highly dependent on the degree of supersaturation, and therefore under our conditions the formation of liquid is very slow. However, with appropriately sensitive techniques it is possible to monitor the reequilibration of clusters of size n = 1-5 on the millisecond time scale.I0 The rate of reaction 1 is governed by differential equations containing terms proportional to [(CC13F)n][CCI3F]. Under our conditions (for a mixture of 25% CC13Fin 75% Ar, the postshock e Abstract published
in Advance ACS Absrracts, November 15, 1993.
0022-365419312097- 12670$04.00/0
conditions were 185.9 Torr of CC13F at a temperature of 336.0 K), a significant number of dimers are present at equilibrium. Hence, the rate equations are highly nonlinear, and it is not surprising to observe a rich variety of behavior ranging from monotonic, to oscillating, to chaoticvariation of the rate of change of temperature, depending on such parameters as the degree of dilution in Ar. Shock-compressed CCl3F traveled down a stainless steel rectangular cross section tube and arrived at a series of stations traversed by narrow He-Ne laser beams of low power ( 5 mW). The chemical reaction's thermal history was spread out spatially, thus exposing a density gradient and hence a refractive index gradient to the laser beams, one after the other. The time intervals for passage between stations of known separation gave us the shock velocity which calibrated the reaction temperature and pressure. The intercepted laser beams were deflected in proportion to the local refractive index gradient, and this time-dependent deflection off the center of a quadrant-connected photodiode generated a time-varying electrical signal with a response of 0.01 ps. Because the signal is proportional to the density gradient, it is a reflection of the chemical reaction's rate of heat release and is most sensitive at the early stages of the reaction. One of the three schlieren station's signals could be examined for the detailed behavior and is presented in this work. The millivolt analog signal was captured using a modified apparatus that is described in ref 3. The heart of this apparatus is an 8-bit AID (analog-to-digital) CompuScope data collection card (manufactured by Gage Applied Sciences, Montreal, Quebec, Canada) installed in a personal computer. The millivolt analog signal was sampled at 40 million samples per second (msps), and the data that are shown in this work are for 32 768 data points, Le., for 8 19.4 ps. The signal that is plotted in this work is a subset of these collected data points. Every eighth data point was kept, and the wavelet analysis was performed on these 4096 data points. This procedure of keeping only every eighth data point could be considered a smoothing of the original data. In fact, a wavelet analysis was also performed on a data set that was obtained by keeping every second data point of the original 32 768. That wavelet analysis did not work as well as the wavelet analysis of the 4096 data set, which are the only plots that are presented here. This is an indication that the wavelet analysis is adversely affected by a certain level of noise. The usual signal from the schlieren detection system is a sharp 0 1993 American Chemical Society
The Journal of Physical Chemistry, Vol. 97, No. 49, 1993 12671
Letters a
si-02
(e)
FFT V 3 9iiO3
C GEiOO
4
Of-04
8 O b 0 4
3 5E-02-
0 0i400
4
OE-04 8 O E - 0 4
Time
Time
Hz
.
9E-0
2 5;-02
2 5E-02
3 9Et01
0 GEtGO
4
OE-04
Time
8 OE-04
o ortoo
4 OE-04
a
OE-04
Time
Figure 1. Plot of the millivolt signal from the schlieren detection system
-1
OEfO3 0 OEt03
6 4Et05
1
3Et36
of the shock tube apparatus: (a) the "smoothed" data, (b) dc(-25 X IO-j)-shifted voltage, (c) waveform reconstructed from the wavelet analysis, and (d) waveform reconstructed from edited details D,, Dz,D3, and Dq and the unedited details of the analysis.
Hz Figure 2. Fast Fourier transform of the complete millivolt signal (32 768 data points) lasting 819.2 ps: (a) the full scale plot for 0-1.3 MHz; (b) a vertical scale expansion by a factor of 100.
negativegoing peakdue to the passageof theshock front, followed by a much larger positive peak which decays back to zero.8 For CC13F, in contrast, the observed signal has been a quasiperiodic one. The usual features are present as well, but there are sometimes small-scale oscillations and other times very largescale oscillations superimposed on the decaying part of the curve.' The intriguing series of spikes observed in this study and shown in Figure l a are only one manifestation of the rich variety of behavior which can beobserved. Thecomplex behavior is believed to be due to the formation and breakup of n-mers of CC13F. This has been discussed in a preliminary manner by Cheng et a].,' and a complete explanation is being presented elsewhere.lO Here we concentrate on describing our novel procedure for extracting the frequency components of the oscillations. The signal that is used in this work was from an experiment that is presented in detail in ref 3. The velocity of the shock wave that produced the signal, used here, was only slightly larger than the velocity of an experiment in which a mixture of 25% CC13Fin 75% Ar produced no oscillations (postshock conditions were 148.7 Torr of CC13F and a temperature of 317.9 K) but instead merely showed the usual slow decay back to zero. Thus, the temperature jump that the system was involved in was only slightly larger ( A T = 18.1 K). Nevertheless, this slight increase had a dramatic effect. Since it is the oscillations that are the interesting part of the signal the portion of the signal that was caused by passage of the shock front was not subjected to the analysis. The wavelet analysis presented in this work is accomplished with the Mallat pyramidal algorithm,llJ2 which is a multiresolution method. In other words, the signal is represented at different resolutions, one being coarser than the next, with the difference between the information present at one resolution and the next resolution being termed the detail, Dj, which Mallat symbolized by D,. Thus, the detail is a data set that contains the information missing in the signal representation at two different resolutions. The signal can be viewed as being a collection of data, f, that are approximated vectorially. f is equivalent to Vo = VI @ Wl, with @ meaning orthogonal sum. This could be continued and represented with 5-1 = V,@ W,. The vector, 6,is equivalent to Afand the vector, Wj, is equivalent to Df Note that VI has half of the data points that Vo has, and the detail Dj is the difference between the signal at the two resolutions, j and j - 1. This decomposition of the signal into lower and lower resolution representation is stopped a t j = 8 for this work. The reconstruction of the original signal is achieved by a reversal of the above described process.
A two-parameter, (m,n), set of functions called the scale function wavelet set, 4mn,$mn, can be defined and used as basis functions for V, and W,, respectively. For example, &,,(r) = 2-m/24(2"t - n) and h n ( t ) = 2-m/2$(2-mt - n) where the m parameter dilates the functions and the n parameter translates the functions. The choice of wavelet, $ ( t ) , is large, but the orthogonal basis sets, as defined by Da~bechies,~3-15 result in functions that are fractal in nature for the lower numbered wavelets, Le., db04, where db stands for Daubechies. As the number of the wavelet increases, the functions are less fractal and are considered to be smoother.15 Here the Daubechies 08 (db08) wavelet is used, and the coefficients for these wavelets are tabulated in ref 13 and 15. The wavelet analysis has been described by Mallat12 as a process of partitioning the signal into different octaves or frequency bands, and this process is illustrated in this work. In Figure 1, the millivolt signal is plotted a t various stages of the analytic process. Each of the plots in Figure 1 has a reference line drawn in at 0.0 V. Figure l a is the millivolt signal that has been reduced to 4096 of the original 32 768 data points. In Figure l b the millivolt signal has been dc shifted by -25 X 10-3 V. This shifting was done to decrease the amplitude of the zero-frequency component of the FFT of the millivolt signal. This is a common FFT practice as described by Brigham.6 Figure lc,d will be discussed below. The FFT of the millivolt signal, of Figure 1b, is shown in Figure 2. The maximum amplitude is less than 7900 (a) while the zerofrequency component amplitude was approximately 4. Figure 2b is a vertical scale expansion of Figure 2a by a factor of 100. There are several frequencies present, and the frequencies that our new technique resolves are present here also but with smaller amplitudes. Thus, they are buried in the noise of the millivolt signal. The wavelet analysis was then performed on the dc-shifted 4096 data point signal of Figure 1b. Since each of the resulting detail data files has a different number of data points, each file is filled with zeros added at the end to make up a total of 4096 points, and then an FFT is calculated. To resolve separate frequencies, a procedure is used, that was developed in ref 3, of setting a threshold value of the amplitude of the detail data. If the detail amplitude is larger than the threshold value, the data point value is set to zero. The threshold value is decided upon by observing the amplitude of the detail plot. Usually the waveform in the detail will have a central set of values with only a few values exceeding this amount. This center value is then
Letters
12672 The Journal of Physical Chemistry, Vol. 97,No. 49, 1993 2 OE-02
5 DE-03
D1 0 OEtOO
0 OEtOO
1
1 j
1
-2 OE-021 0 OEtOO
2 F5E F31
03
4 OE-04
0 OEtOO
8 Of-04
Time
4
OE-04
2 5;-01, 04
0 OEtOO
0 OEtOO
-2 O
E 0 2 1 0 O E I O O 4 OE-04 8 OE-04
Time
-5
OE-03-
o
FFT D2
octo0
1
4 O E - O ~
,
’
5 F I E - 0F2
m
T
D
3
L
i
2 5E-02
‘
OE-02 0 OEtOO
- 1
8 OE-04
D
2E-01
Time
02
r
2E-011
6 4Et05
OE-OS 0 OEtOO
-1 1 3Et06
Hz
6Et05
3 2Et05
Hz
5
(b)1
1
r 1
~
‘E-02
FFT
D4
2 5E-02
I
8 OE-O~
Time
Figure 3. Plots of detail files: (a) D I , (b) D2,(c) D3,and (d) D4.
usedas the thresholdvalue in the FFTprogram. Several threshold values near this selected value are chosen to see whether the FFT output changes. If there are only minor changes for an order of magnitude change of the threshold value, and if the amplitude of the 0-Hz component is relatively small (a tool used in FFT analysis), then that threshold value is kept. The edited detail data file is also saved in this process to be used in the second reconstruction of the original signal. The procedure uses eight detail files to reconstruct the original signal. The reconstructed signal is shown in Figure ICand is, for all intents and purposes, identical to Figure 1b. (In fact, the data files agree to the three significant figures supplied by the original 8-bit (A/D).)In the second reconstruction the first four edited detail files and the four remaining unedited detail files are used to produce the waveform in Figure Id. Notice that Figure Id has a significant signal-to-noise improvement over Figure 1b and reveals an underlying oscillation which is otherwise not evident in the original Figure la,b because of the high amount of noise present. Since the second reconstructed waveform is different from the original waveform, in theory, it should be possible to repeat the process and pull out further information. In fact, this is not the case at all. A subsequent wavelet-FFT of the second reconstructed waveform results in the already known frequencies and no new information. Although it does provide an increase in the signal-noise-ratio,unfortunately the decrease of the noise is not as dramatic as that shown from a comparison of Figure Ic,d. Thus, there is a sense of diminishing returns. In Figure 3 the detail files are presented with vertical scales which are the same for DI and D2 and the same for 0 3 and 0 4 . The detail data points are plotted, and a line is drawn between them to guide the eye. There is too much information present to be able to detect a single frequency in the waveform, and that is one reason for performing an FFT on the waveform. The threshold values used for the FFT were a3.0X 10-3 for D1,f5.0 X le3 for D2,f 3 . 0 X IC3for D3, and f2.0 X IC3for Dd. In Figure 4 the FFT of the details, as described above, are plotted. In Figure 4a frequencies of 1.18 MHz and 500 KHz are observed while in Figure 4b 250 and 63.2 kHz are the observed frequencies. In Figure 4c,d the frequencies are 251, 94.2 and 61 -4kHz. These are dramatically more pronounced than in the simple FFT’s of Figure 2. The fact that the edited detail files reconstruct the original signal with a removal of some noise is an indication that the frequencies detected are an important component of the original signal and not the noise itself, whether those frequencies are chemically significant or not (see below). A wavelet analysis of a 16 382 data point tile was also calculated, but only a few of the higher frequencies were observed and then only when the threshold values were set extremely low. Thus the
Hz
HZ
Figure 4. Plots of the fast Fourier transform of the detail data files with threshold settings slightly different for each detail (sec text): (a) FFT of DI,(b) FFT of D2,(c) FFT of D3,and (d) FFT of D4.
“smoothed” data file, that the 4096 data point file represents, had less noise than the 16 384 data point file, and the wavelet analysis was more effective on the “smoothed” data. Such preliminary smoothing is apparently advantageous. The success of the WFFT analysis of this work in “denoising” the signal can be explained by referring to Mallat and Hwange4 They have developed a special algorithm to push the “denoising” process even further. However, their basic explanationwill suffice for this work. Mallat and Hwang4point out that when the wavelet analysis is performed on the derivative of a signal, the higher details actually contain information pertaining to the location of the edges of the original waveform. In our work the higher numbered details are not edited. The second and third reconstruction processes, of this work, edit only the lower numberdetails. It is evident, from Figure 4, that the higher frequenciesare present in the lower numbered details, and decreasing the amplitude of some of the points in the detail eliminates high-frequency noise. Since the higher numbered details are unedited, the relatively low frequencies needed to build up the original signal are unaffected. Thus, as pointed out earlier, this process has broken the original signal into frequency bands, some of which have been edited to remove noise, while the remainder have been left unaltered. In order to identify the chemical oscillations with a set of frequencies, a WFFT was performed on the base line. A highfrequency spectrum very similar to the one in Figure 4, but with the low frequencies of Figure 2 absent, was observed. (In fact, the He-Ne laser produces a 0.1% rms background signal with frequencies in the range of 10 kHz to 1 MHz.) The chemically interesting information is present in the low-frequency portion. (A FFT of the denoised signal Figure Id showed the same lowfrequency components as in Figure 2.) The WFFT, as used in this study, works best for high frequencies; for low frequencies it does not present any particular advantage over the normal FFT (of e.g. Figure 2). Nevertheless, the WFFT technique has succeeded in processing a noisy data set and has effectively done the work of time-averaging repeated experiments. Time-averaging works for experiments whose initial conditions can be reproduced exactly. If, as in this case, initial conditions are not reproducible or the signal is a fast transient with a long time interval (- 1 h) between experiments, time averaging is not practical, and therefore it is much more advantageous to pull out the existing information present in the single experiments. There are other algorithms, such as developed by Delprat et al.,l6 for removing frequencies from noise, then ratioing out these frequencies, and finally revealing new frequencies. A combination of approaches may be useful in the future.
Letters This seems to be the first application of a wavelet analysis to chemical kinetics and one of the first applications of a combined wavelet-FFT (WFFT) analysis to the processing of signals in science. The authors would like to close this study by stating that, although the wavelet code used in this work was obtained from an anonymous ftp n0de,17 other code exists.I7
Acknowledgment. The authors gratefully acknowledge a grant of computer time from the University of Ottawa. We also acknowledge the National Science and Engineering Research Council (Canada) for partial funding of this work. We thank Jacques Froment and Stephane Mallat of the Courant Institute of Mathematics, New York University, for use of their wavelet code. D.N.S.P. thanks Zhikai Cheng for hours of useful discussion on wavelets, chaos, and homogeneous nucleation. D.N.S.P. also thanks the Faculty of Graduate Studies at University of Ottawa for a partial scholarship.
References and Notes (1) Combes, J. M., Grwmann, A., Tchamitchian, Ph., Eds. Wauelets; Springer-Verlag: New York, 1989.
The Journal of Physical Chemistry, Vol. 97, No. 49, 1993 12673 Permann, D.; Hamilton, I. Phys. Rev. Lett. 1992, 69, 2607. Permann, D. Ph.D. Thesis, university of Ottawa, 1993, submitted. Mallat, S.; Hwang, W. L. IEEE Trans. I f o . Theory 1992,38,617. Cartwright, M. Fourier Methods; Ellis Horwood: Chicheater, 1990. Brigham, E. 0. The Fast Fourier Transform and its Applications; Prentice-Hall: Englewood Cliffs, NJ,1988. (7) Cheng, Z.; Hamilton, I. P.; Teitelbaum, H.J. Phys. Chem. 1991,95, (2) (3) (4) (5) (6)
4931. (8) Kiefer, J. H.; Lutz, R.W. J. Chem. Phys. 1966, 44,688. Dove, J.; Teitelbaum, H. Chem. Phys. 1974,6, 431. (9) Carruthers, C.; Francoeur, K.;Teitelbaum, H. In Shock Tubes and
Waues: Proceedings of the 16th International Symposium on Shock Tubes and Waues;Grhig, H., Ed.; VCH Publishers: Weinhem, 1988; p 327. (IO) Cheng, Z.; Teitelbaum, H. To be published. (11) Mallat, S. G. Technical Report No. 412; Courant Institute of Mathematics, New York University: New York, 1988. (12) Mallat, S. G. IEEE Trans.Pattern Anal. Machine Intelligence 1989, 11, 614. (13) Daubechies, I. Commun. Pure Appl. Math. 1988, 41, 909. (14) Daubechies, I. IEEE Trans. Info. Theory 1990, 36, 961. (1 5 ) Daubechies, I. Ten Lectures on Wauelets;SIAM: Philadelphia, 1992. (16) Delprat, N.; Escudie, B.; Guillemain, P.; Kronland-Martinet, R.; Tchamitchian, P.; T o r r h n i , B. IEEE Trans. Info. Theory 1992, 38, 644. (17) Froment-Mallat wavelet code is available from the anonymous ftp site cs.nyu.edu in subdirectory pub/wave/megawave. Other wavelet code is available from: Numerical Recipes in C or Numerical Recipes in Fortran; Cambridge University Press: Cambridge, UK, 1992. There are also example books and disks of code available from Cambridge University Prcss.
