Comment on ‘ Group Velocity Measurement of Surface Waves by Wavelet Transform ’

1 Department of Geophysics, Graduate School of Science, Kyoto University, Sakyo, Kyoto, Japan * Corresponding author address: Prof. Taishi Okamoto, Department of Geophysics, Graduate School of Science, Kyoto University, Sakyo, Kyoto, Japan; E-mail: dollar@kugi.kyoto-u.ac.jp doi: 10.3319/TAO.2007.18.4.681(T) Yamada and Yomogida (1997) applied the discrete wavelet transform (DWT) to group velocity measurements for the first time. Although their study is one of the pioneering works in application of DWT to seismological analysis, their method gives an inaccurate value as a group velocity in some cases and requires modification. In this report, we point out the problem and propose a modified DWT method for overcoming the problem. In our method, DWT is carried out not for an analysed signal itself but for its complex envelope (Farnbach 1975). A computation algorithm for DWT coefficients for our method is given and shown to be almost the same as that by Yamada and Ohkitani (1991). The influence of the difference between the conventional method and our method on identification of group arrival times of a wave is also shown by a numerical experiment. If analysts want to identify group arrival times using DWT, our method must be adopted instead of the conventional method. (


INTRODUCTION
The discrete wavelet transform (DWT) and the continuous wavelet transform (CWT) are methods for time-frequency analysis of time history data, and a squared DWT coefficient of the data is supposed to represent energy in corresponding time and frequency range (e.g., Yamada and Ohkitani 1991), because some types of wavelet functions can form an orthonormal basis for the function space L 2 (R), the set of square-integrable functions.Yamada and Yomogida (1997) applied the DWT to group velocity measurement of dispersed surface waves based on this idea.They measure group velocities, defining a time for which an absolute value of a DWT coefficient for an observed or synthetic record has a relative maximum as a group arrival time of a wave having the same frequency as each wavelet.
In the present report, however, we show that their method leads to an inaccurate value as a group arrival time in some cases.No reports pointing to the same problem can be found thus far.We first explain what wavelet coefficients represent and show that their method adopted for the identification of the group arrival time should be modified theoretically.Next, we propose a modified DWT method overcoming the problem, and introduce a computation algorithm for DWT coefficients for our method.Finally, we make a numerical experiment, illustrating the influence of the difference between their method and our method on the identification of group arrival times of a wave and that in some cases their method leads to an inaccurate value as a group arrival time, whereas our method always gives a correct value.

Comparison with the Work by Pyrak-Nolte and Nolte (1995)
In this section, we clarify the problem in the application of the conventional DWT by Yamada and Yomogida (1997).
First, we compare their application with that by Pyrak-Nolte and Nolte (1995).The former uses the DWT and the latter the CWT for the same purpose, identification of group arrival times of dispersed surface waves.Except for the difference, the methods used in the two studies seem to be the same, because in both studies a time for which an absolute value of a wavelet coefficient for a dispersed wave has a relative maximum is interpreted as a group arrival time of a wave having the same frequency as each wavelet.However, they are by no means the same, although it is not important whether CWT is used or DWT is used.The decisive difference between the two methods is the fact that the former uses a real-valued wavelet, and the latter a complex-valued wavelet.The influence of this difference on wavelet coefficients is described below.

What Do Wavelet Coefficients Represent?
When we calculate wavelet coefficients for a signal, there are four choices in general: (A) CWT coefficients with a real-valued wavelet, (B) CWT coefficients with a complex-valued wavelet, (C) DWT coefficients with a real-valued wavelet, and (D) DWT coefficients with a complex-valued wavelet.
Among them, (C) is adopted by Yamada and Yomogida (1997) and (B) by Pyrak-Nolte and Nolte (1995).We explain the relation between the above four kinds of coefficients and what they represent.

(A) CWT Coefficients with A Real-Valued Wavelet
Consider the case (A).We assume that an analyzed signal f (t) is a real-valued function of time throughout this paper.If ψ( ) t is a real-valued analyzing wavelet whose Fourier trans- form Ψ( ) ω satisfies the admissibility condition: where ω denotes the angular frequency, then the CWT coefficient of f is defined as: where * denotes the complex conjugate and a and b scale and shift factor, respectively.Note that the equation ( 1) has the form of convolution, that ψ( ) t and ψ a b t , ( ) defined as: have a band-pass property because of the admissibility condition, and that T a b f , ( , ) ψ is a real-valued function of b for fixed a. Thus, the CWT coefficients with a real-valued wavelet represent the bandpassed signal of f (t) by ψ a b t , ( ) .

(B) CWT Coefficients with A Complex-Valued Wavelet
Next, consider the case (B).The CWT coefficient with complex-valued wavelet ψ c t ( ), which is an analytic function satisfying : where ψ( ) t is a real-valued analyzing wavelet and ψ H t ( ) is the Hilbert transform of ψ( ) t given by: and is a complex-valued function of b for fixed a.As shown in Appendix A, it follows that the equality: where T a b c f , ( , ) is the complex envelope of T a b f , ( , ) ψ , holds for all a and b.As described in (A), T a b f , ( , ) ψ is the bandpassed signal.Thus the CWT coefficient with a complex-valued wavelet is the complex envelope of the bandpassed signal.

(C) DWT Coefficients with A Real-Valued Wavelet
Now consider the case (C).The DWT coefficient of f with real-valued wavelet is defined as: where j and k denote scale and shift factor, respectively.It is clear that the DWT coefficient is connected with the CWT coefficient of f as: , which implies that the DWT coefficient with a real-valued wavelet is the sampled signal after being bandpassed.

(D) DWT Coefficients with A Complex-Valued Wavelet
Finally consider the case (D).The DWT coefficient of f with complex-valued analytic wavelet ψ c t ( ) satisfying (2) is given by: From ( 3) and ( 4), ,ψ is connected with the CWT coefficient of f as: Thus the DWT coefficient with a complex-valued wavelet is the sampled complex envelope of the signal after being bandpassed.

A Problem in the Application by Yamada and Yomogida (1997)
To identify group arrival times of a surface wave, we have to make use of an amplitude of a complex envelope of a bandpassed signal (Dziewonski et al. 1969), corresponding to (B) or (D) in the section 2.2.Hence we conclude that Yamada and Yomogida (1997)'s method, in which DWT coefficients with a real-valued (Meyer-Yamada) wavelet are used [case (C)], must be modified, unlike the work by Pyrak-Nolte and Nolte (1995) who make use of CWT coefficients with a complex-valued analytic (Morlet) wavelet [case (B)].

A MODIFIED DWT METHOD
As described in section 2, we have to use a complex-valued analytic wavelet to identify group arrival times of a wave.However, even a real-valued wavelet such as a Meyer wavelet can be used for this purpose if an analysed signal f is converted into its complex envelope f c , because the equality: holds as shown in Appendix A and substituting a = 1 / 2 j and b = k / 2 j into (5) yields: Therefore, we modify Yamada and Yomogida (1997)'s DWT method as follows: We compute DWT coefficients for a complex envelope of an original signal with a real-valued orthonormal wavelet, prepare a plot of their amplitude d j k f c , , ψ , and identify a time for which the amplitude has a relative maximum as a group arrival time.Note that DWT coefficients are not computed for an original signal with a complex-valued wavelet, and that the wavelet still has the orthonormality.And it should also be noted that Yamada and Ohkitani (1991)'s computation algorithm for DWT coefficients for Meyer wavelet is available only if an analysed signal is a real-valued function.However, as shown in Appendix B, if the analysed signal is an analytic function, the DWT coefficient for Meyer wavelet can be calculated as follows: 1 obtained by sampling a real and T-periodic function f(t) with a sampling interval of T/N, where N = 2 n and n N ∈ , are given and if we expand the complex envelope f c (t) of f as: , where ψ( ) t is a Meyer wavelet and: ∈ becomes a 2 j -periodic sequence of k for fixed j, and its fundamental period part { } , , , , ..., where F S is the N-point discrete Fourier transform of the data sequence: Note that eq. ( 6) is almost the same as eq.( 20) in Yamada and Ohkitani (1991), and therefore d j k f c , , ψ can be computed within the conventional computation scheme.

NUMERICAL EXPERIMENTS
In this section, we make a comparison between DWT coefficients obtained from the conventional method and our method proposed in section 3 by some numerical experiments.
Figure 1 shows horizontal two component broadband velocity seismograms from 2004 great Sumatra-Andaman Earthquake, observed at F-net NMR (Nemuro), north Japan.The zero time corresponds to 005240 UT 26 December 2004.
We apply both our DWT and conventional DWT method to these seismograms.The wavelet used here is a Meyer-Yamada wavelet (Yamada and Ohkitani 1991), the same as used in Yamada and Yomogida (1997).Figures 2 and 3 show absolute values of DWT coefficients for j = 8 using conventional DWT and our DWT method, respectively.The DWT coefficients are calculated for five different zero times, corresponding to 00:52:20, 00:52:30, 00:52:40, 00:52:50, and 00:53:00, respectively.Note that the level j = 8 in the DWT analysis corresponds to the period range of 96 to 384 sec.
In Figs. 2 and 3, fundamental (R1 and R2) and higher-mode (X1 to X3) Rayleigh waves and fundamental mode Love waves (G1 and G2) can be found and their theoretical group arrival times for PREM are indicated.Note that the absolute value of the DWT coefficient in Fig. 3 always has a relative maximum for the theoretical group arrival times.On the other hand, in Fig. 2, the DWT coefficient has a significantly different value for a slightly different zero time, and its absolute value does not always have a relative maximum for the theoretical group arrival times.
We try to measure the group velocity of the Rayleigh waves using the group arrival times for additional two frequency bands ( j = 7 and j = 9) by the conventional method and our method, and show the results in Fig. 4. Figures 4a and b show the estimated group velocities of fundamental mode Rayleigh waves using Yamada and Yomogida (1997)'s method and our method, respectively.The group velocities, denoted by red stars, are measured for five different zero times, which are the same as in Fig. 3.While our method always gives correct values as group velocities, the conventional method gives inaccurate values in some cases.This figure clearly shows that the conventional method may lead to an inaccurate result in the determination of a group arrival time of a dispersed surface wave and our DWT method has to be adopted.
We performed the same analysis for other two events.Figure 5 shows a transverse component broadband velocity seismogram recorded at F-net TTO (Takato) station, central Japan, for Macquarie Island Earthquake of 23 December 2004.The zero time is 145930 UT 23 December 2004, which corresponds to the centroid time.
Figures 6 and 7 show absolute values of DWT coefficients for j = 8 using conventional DWT and our DWT method, respectively.The DWT coefficients are calculated for five      different zero times, corresponding to 14:59:10, 14:59:20, 14:59:30, 14:59:40, and 14:59:50, respectively.In Figs. 6 and 7, fundamental mode Love waves (G1 and G2) can be found and their theoretical group arrival times for PREM are indicated.Note that the absolute value of the DWT coefficient in Fig. 7 always has a relative maximum for the theoretical group arrival times.On the other hand, in Fig. 6, the DWT coefficient has a significantly different value for a slightly different zero time, and its absolute value does not always have a relative maximum for the theoretical group arrival times.
Our modified DWT method also works for regional data with higher frequency contents.Figure 8 shows a transverse component broadband velocity seismogram recorded at F-net FUJ (Fujigawa) station, located on the mountainside of Mt.Fuji, central Japan, for the Noto Peninsula Earthquake of 25 March 2007.The zero time is 004200 UT 25 March 2007, which corresponds to the hypocentral time.The epicentral distance is about 270 km and the focal depth is 11 km.In Fig. 8, a prominent phase with a dominant period of 6 -7 sec can be found around 90 sec.Since group velocity of a Love wave with this period range in this region is known to be about 2.8 -3.2 km sec -1 , the phase can be considered to be a Love wave, whose energy is trapped within the crust.Figures 9 and 10 show absolute values of DWT coefficients for j = 13 using conventional DWT and our DWT method, respectively.The DWT coefficients are calculated for five different zero times, corresponding to 00:42:00, 00:42:01, 00:42:02, 00:42:03, and 00:42:04, respectively.The level j = 13 in the DWT analysis corresponds to the period range of 3 to 12 sec.Note that the difference of the zero times is only 1 sec.In Figs. 9 and 10, we give a theoretical group arrival time for the Love wave, with a group velocity of 3.2 km sec -1 .Again we find that the absolute value of the DWT coefficient in Fig. 10 always has a relative maximum for the theoretical group arrival time, whereas that in Fig. 9 has a significantly different value for a slightly different zero time and does not always have a relative maximum for the theoretical group arrival time.

SUMMARY
In this report, first we have clarified what wavelet coefficients represent theoretically and pointed out a problem in Yamada and Yomogida (1997)'s application by comparison with Pyrak-Nolte and Nolte (1995)'s work.
Next, to overcome the problem, we have modified their DWT method and given a computation algorithm for the DWT coefficients for our method.In our method, DWT coefficients are computed for a complex envelope of an original signal with a conventional real-valued orthonormal wavelet.Note that DWT coefficients are not computed for an original signal with a complex-valued wavelet, and that the wavelet still has orthonormality.The computation algorithm for our method is almost the same as the conventional one.
Finally, we have illustrated the influence of the difference between the conventional method and our method on the determination of group arrival times of a wave by some numerical experiments.Our DWT method always gives a correct value as a group arrival time of a surface wave, while the conventional DWT analysis results in inaccurate values in some cases.
Although Yamada and Yomogida (1997)'s pioneering work is of great importance in application of DWT to seismological analysis, modification is required as described in this paper.If analysts want to identify group arrival times using DWT, our method must be adopted instead of the conventional method.

APPENDIX A: PROOF OF (3) AND (5)
There are two keys to the proof.One is the Parseval's relation: where A( ) ω and B( ) ω are the Fourier transform of a(t) and b(t).The other is the fact that generally a real-valued function g(t) is connected with its complex envelope g c (t) as: , (Katukura et al. 1989) where G( ) ω and G c ( ) ω respectively denote the Fourier transform of g(t) and g c (t), and U( ) ω is the unit step function defined as: From these two keys, T a b f c , ( , ) Hence: = .
In the same way, we have:  where F s is given by ( 7) and similarly

Fig. 1 .
Fig. 1.The radial (above) and transverse (bottom) component broadband velocity seismograms observed at F-net NMR (Nemuro), north Japan from the 2004 great Sumatra-Andaman earthquake.The time is measured from 005300 UT 26 December 2004.

Fig. 2 .
Fig. 2. The absolute values of DWT coefficients (bars) for j = 8 for the (a) radial and (b) transverse component seismograms shown in Fig. 1, calculated by Yamada and Yomogida (1997)'s method.The coefficients are calculated for five different zero times.Vertical lines indicate the theoretical group arrival times for fundamental (R1 and R2) and higher-mode (X1 to X3) Rayleigh waves and fundamental mode Love waves (G1 and G2) in PREM.The DWT coefficient has a significantly different value for a slightly different zero time, and its peaks do not always correspond to the theoretical group arrival times.

Fig. 3 .
Fig. 3.The absolute values of DWT coefficients (bars) for j = 8 for the (a) radial and (b) transverse component seismograms shown in Fig. 1, calculated by our method.The coefficients are calculated for five different zero times.Vertical lines are the same as in Fig. 2. All peaks of the DWT coefficients correspond to the theoretical group arrival times.

Fig. 4 .
Fig. 4. Estimated group velocities of fundamental mode Rayleigh waves using (a) Yamada and Yomogida (1997)'s method and (b) our method.The group velocities (red stars) are measured for five different zero times, which are the same as in Fig. 3. Red curve denotes theoretical group velocity for fundamental mode Rayleigh wave in PREM.While our method (b) always gives correct values as group velocities, the conventional method (a) gives inaccurate values in some cases.

Fig. 5 .
Fig. 5.The transverse component broadband velocity seismogram observed at F-net TTO (Takato), central Japan from the Macquarie Island earthquake.The time is measured from 145930 UT 23 December 2004.

Fig. 6 .
Fig. 6.The absolute values of DWT coefficients (bars) for j = 8 for the transverse component seismograms shown in Fig. 4, calculated by Yamada and Yomogida (1997)'s method.The coefficients are calculated for five different zero times.Vertical lines indicate the theoretical group arrival times for fundamental mode Love waves (G1 and G2) in PREM.The DWT coefficient has a significantly different value for a slightly different zero time, and its peaks do not always correspond to the theoretical group arrival times.

Fig. 7 .
Fig. 7.The absolute values of DWT coefficients (bars) for j = 8 for transverse component seismograms shown in Fig. 4, calculated by our method.The coefficients are calculated for five different zero times.Vertical lines are the same as in Fig. 5.All peaks of the DWT coefficients correspond to the theoretical group arrival times.

Fig. 8 .
Fig. 8.The transverse component broadband velocity seismogram observed at F-net FUJ (Fujigawa), central Japan from the Noto Peninsula earthquake.The time is measured from 004200 UT 25 March 2007.

Fig. 9 .
Fig. 9.The absolute values of DWT coefficients (bars) for j = 13 for the transverse component seismograms shown in Fig. 7, calculated by Yamada and Yomogida (1997)'s method.The coefficients are calculated for five different zero times.The red vertical line indicates the theoretical group arrival time for a Love wave with a group velocity of 3.2 km sec -1 .The DWT coefficient has a significantly different value for a slightly different zero time, and its peak does not always correspond to the theoretical group arrival time.

Fig. 10 .
Fig. 10.The absolute values of DWT coefficients (bars) for j = 13 for transverse component seismograms shown in Fig. 7, calculated by our method.The coefficients are calculated for five different zero times.The red vertical line is the same as in Fig. 8.The peak of the DWT coefficients always correspond to the theoretical group arrival time.

T
Since f (t) is assumed to be a real and T-periodic function, its Hilbert transform f H (t) is also a real and T-periodic function, and for their discrete Fourier transforms, we have: 1970).Then using the Eq.(19) inYamada and Ohkitani (1991)  for f and f H , we can calculate d j k f c , , ψ as follows (for simplicity, T = 1 is assumed):