European Journal of Scientific Research ISSN 1450-216X Vol.32 No.3 (2009), pp.314-328 © EuroJournals Publishing, Inc. 2009 http://www.eurojournals.com/ejsr.htm
Comparing the Performance of Fourier Decomposition and Wavelet Decomposition for Seismic Signal Analysis
Chik, Z Department of Civil & Structural Engineering Faculty of Engineering & Built Environment Universiti Kebangsaan Malaysia 43600 UKM Bangi Selangor DarulEhsan, Malaysia E-mail: staohidul@yahoo.com or irzamri@eng.ukm.my Tel: +60(3)8921 6200/6213; Fax: +60(3) 8921 6147 Islam, T Department of Civil & Structural Engineering Faculty of Engineering & Built Environment Universiti Kebangsaan Malaysia 43600 UKM Bangi Selangor DarulEhsan, Malaysia Rosyidi, S.A Department of Civil & Structural Engineering Faculty of Engineering & Built Environment Universiti Kebangsaan Malaysia 43600 UKM Bangi Selangor DarulEhsan, Malaysia Sanusi, H Department of Civil & Structural Engineering Faculty of Engineering & Built Environment Universiti Kebangsaan Malaysia 43600 UKM Bangi Selangor DarulEhsan, Malaysia Taha, M.R Department of Civil & Structural Engineering Faculty of Engineering & Built Environment Universiti Kebangsaan Malaysia 43600 UKM Bangi Selangor DarulEhsan, MALAYSIA Mustafa, M.M Department of Civil & Structural Engineering Faculty of Engineering & Built Environment Universiti Kebangsaan Malaysia 43600 UKM Bangi Selangor DarulEhsan, Malaysia Abstract Different decomposition techniques are used to estimate the amplitude and phase of the seismic signal corresponding to frequency in geotechnical characteristics. In this paper,
Comparing the Performance of Fourier Decomposition and Wavelet Decomposition for Seismic Signal Analysis the Fourier based decomposition’s performance is evaluated and compared with the performance of Wavelet decomposition. The performance of Fourier decomposition and Wavelet decomposition on seismic signal analysis are compared through Matlab Programming to reveal their fitness. Fourier decomposition techniques have found to be promising to improve significantly the accuracy and reliability of the approximate signal by minimizing the integral square error. Wavelet decomposition leads to improvement in the signal to noise ratio (SNR) but its computational complexity is high and restrained for high level wavelet de-noising. The significance of this work is to obtain the idea about convenience and limitations of Fourier and wavelet decomposition on seismic signal for geotechnical research concentration. Keywords: Fourier decomposition, Wavelet decomposition, Seismic signal, De-noising.
315
1. Introduction
The analysis of complex time-domain signal is difficult because of the immense multi-dimensional data sets involved. Decomposition of a complex wave into a time variable-frequency format permits analysis and display of each frequency component in a unique and continuous format. Wavelet decompositions (Strang, 1993) and windowed Fourier decompositions (Unser, 1994) have emerged in the last few years as tools of great importance for the analysis of signals in which local frequencies can be extracted (Torresani, 1993). For seismic signal analysis, the fast convergence of the Fourier series allows a compact representation of the time-dependent solution. In typical cases amplitude and phase angle can accurately reconstruct the solution in the frequency domain. Even in higher-dimension, FFTs are more remarkably stable, at least if the size of the dimensions is large and the number of dimensions is small (Schatzman, 1996). Significant nonlinear interactions are manifested in the Fourier domain by strong steady streaming or by the generation of high modes. Decomposition of the complex seismic signal into a swept frequency record allows the fast generation of an accurate dispersion curve. The dispersion curve can be evaluated by their pattern of amplitudes and phase angles. The frequencydependent properties of Rayleigh-type surface wave can be utilized for imaging and characterizing the shallow subsurface. Seismic data, being non-stationary in nature, have varying frequency content in time. The conventional method of producing a time-frequency map from seismic signal analysis using the Short Time Fourier Transform (STFT) limits the time-frequency resolution by a pre-defined window length. In contrast, the continuous wavelet transform CWT method does not require preselecting a window length and does not have a fixed time-frequency resolution over the time-frequency space (Sinha et al, 2005). Wavelet analysis consists of a versatile collection of tools for the analysis and manipulation of the signals such as sound and images, as well as more general digital data sets (Yves et al, 1992). The aim of this paper is to reveal the comparison of the performance of Fourier and wavelet decomposition on seismic signal analysis regarding de-noising and reconstruction of signal.
2. Literature Review
Decomposition of recorded wavefields into a swept-frequency format permits identification of most noise by frequency phase and source offset. Decomposition can therefore be used in association with multichannel records to make adjustments to minimize noise during acquisition (Park et al, 1999.) and for separation of polarized surface waves present in the ambient seismic noise (Roueff et al, 2009). The discrete Fourier transform (DFT) has gained a wide acceptance in both academic and industry, because it is faster and it provides excellent performance on a variety of machines (Frigo, 1999).An extension of the generalized spectral decomposition method are presented for the resolution of nonlinear stochastic problems (Nouy and Maitre, 2009). In the case of Windowed Fourier Transform, the time
316
Chik, Z, Islam, T, Rosyidi, S.A, Sanusi, H, Taha, M.R and Mustafa, M.M
and frequency resolution are both fixed which make this approach particularly suitable for the signal analysis with slowly varying periodic or stationary characteristics (Unser, 1994). Due to be easier and faster analysis technique, Fourier Transform is accentuated by many researchers. Guo and Burrus (1997) proposed the algorithm to calculate FFT using discrete wavelet transform(DWT) where FFT has been used to calculate DWT previously. Wavelet theory can be viewed as a modern improvement and extension of the Fourier theory and give a flexible alternative to Fourier method in non-stationary signal analysis. Wavelet theory is one of the most successful tools to analyze, visualize, and manipulate complex non-stationary data, for which the traditional Fourier methods cannot be applied directly (Cho and Chon, 2006). Wavelets and their relatives generated a lot of interests in diverse fields ranging from astronomy to geology to biology as well as statistics and computer science. In each of these fields, the wavelets are applied for data compression, noise removal, feature extraction, classification, and regression (Saito, 2004). For the stability and Reliability, many research concentrated on wavelet decomposition in seismic signal analysis, such as using the Gaussian wavelet transform as a filtering method for ground roll removal (Corso et al, 2003), using the harmonic wavelet to obtain the dispersive phase and group velocities (Park and Kim, 2001; Kim and Park, 2002). Though wavelet transform does not tie so nicely with the concept of band limited like Fourier analysis, this lack of correspondence can be overcome using a notion of scale band signal (Odegard et al, 1992). 2.1. Fourier Analysis Fourier analysis is a mathematical technique which breaks down a signal into constituent sinusoids of different frequencies and transforms our view of the signal from time-based to frequency-based. A linear system permits the superposition which can be used to find the system response to an input signal (received signal) by decomposing the input signal into a sum of signal components, determining the system response to each component, and adding the component responses. A periodic signal can be described by a Fourier decomposition as a Fourier series, i.e. as a sum of sinusoidal and cosinusoidal oscillations. By reversing this procedure a periodic signal can be generated by superimposing sinusoidal and cosinusoidal waves. The general function is ∞ a x(t ) = 0 + ∑ (ak cos(kω0t ) + bk sin(kω0t )) (1) 2 k =1 Even if we are not concerned about system response, we may decompose a signal into a sum of signal components to show individual, important characteristics of the signal. For a more complicated signal we have to decompose it into a sum of components that have readily identifiable, simple characteristics. However, we will be able to approximate the signal, by a linear combination (weighted sum) of a set of signals, that is by the series
ˆ x(t ) = ∑ X nφn (t )
n=1 N
(2)
ˆ where x(t ) is the approximated signal of original signal x(t), signals φ n (t ) are basis signals and Xn is the weighting constant (Raghuwanshi et al, 2006).
2.2. Wavelet Analysis
Unlike a Fourier decomposition which always uses complex exponential (sine and cosine) basis functions, wavelet decomposition uses a time-localized oscillatory function as the analyzing or mother wavelet. The mother wavelet is a function that is continuous in both time and frequency and serves as the source function from which scaled and translated basis functions are constructed. The mother wavelet can be complex or real, and it generally includes an adjustable parameter which controls the properties of the localized oscillation. Wavelet analysis is more complicated than Fourier analysis since one must fully specify the mother wavelet from which the basis functions will be constructed. Nason
Comparing the Performance of Fourier Decomposition and Wavelet Decomposition for Seismic Signal Analysis
317
and Silverman (1994) revealed how to perform one dimensional and two dimensional discrete wavelet transform in the software package Wavethresh and to extract signal from noise through thresholding of wavelet coefficient. In wavelet analysis, a signal(S) is segregated into an approximation (A) and a detail (D) as Eq. (3). (3) S=A1+D1=A2+D1+D2=A3+D1+D2+D3. The approximations are the high-scale, low-frequency components of the signal. The details are the low-scale, high-frequency components. The approximation is then itself split into a second-level approximation and detail, and the process is repeated. For n-level decomposition, there are n+1 possible ways to decompose or encode the signal. In wavelet packet analysis, such the details as well as the approximations can be taken apart. Fourier analysis consists of breaking up a signal into sine waves of various frequencies. Similarly, wavelet analysis is the breaking up of a signal into shifted and scaled versions of the original (or mother) wavelet. The windowed-Fourier analysis coefficients are the doubly indexed coefficients G s ( w, t ) = ∫ s (u ) g (t − u )e − iwu du (4)
R
The analogy of this formula with that of the wavelet coefficients (Hanssen et al, 1994) is obvious 1 t −u C (a, t ) = ∫ s (u )( ) ( ψ )du (5) R a a From a theoretical viewpoint, wavelets are used to characterize large sets of mathematical functions and are used in the study of operators linked to partial differential equations. From a practical viewpoint, wavelets are used in several fields of numerical analysis, making certain complex calculations easier to handle or more precise.
3. Case of Study
3.1. Original Signal Study For this study, the seismic signal collected by Spectral Anaysis of Surface Wave (SASW) testing at Universiti Kebangsaan Malaysia (UKM) in Bangi, Selangor, Malaysia is used shown in Fig. 1. A set of impact sources of various frequencies were used to generate R waves on the pavement surface revealed in Fig. 2a.
Figure 1: SASW testing at Universiti Kebangsaan Malaysia (UKM) in Bangi, Selangor, Malaysia.
318
Chik, Z, Islam, T, Rosyidi, S.A, Sanusi, H, Taha, M.R and Mustafa, M.M
The propagation of the waves were detected using two receiving accelerometers of piezoelectric DJB A/123/E model and the analog signals were then transmitted to a Harmonie 01 dB (IEC 651–804 Type-I) acquisition box and transferred digitally a notebook computer according to Fig. 2b. Several configurations of the receiver and the source spacings were required in order to sample different depths. The measurement configuration of the SASW test used in this study is the mid point receiver spacings.
Figure 2: (a). Various wave sources (b). SASW measurement set up on the pavement surface.
(a)
(b)
The seiemic signal in Fig. 3a collected from our SASW testing is exploited to compare the performance of Fourier Decomposition with wavelet Decomposition. To simplify this complex seismic signal, it is multiplied by the shifted Gaussian Window in Fig. 3b and a simple Gaussian slice in Fig. 4 is extracted. The complex wave can be simplified according to Gaussian model as Gaussian wave shape by the Eq. (6) as
y =
∑
n +1 i=1
a ie
⎡ ⎛ x − b ⎞2 i ⎢−⎜ ⎟ ⎢ ⎝ ci ⎠ ⎣
⎤ ⎥ ⎥ ⎦
(6)
where a is the amplitude, b is the centroid (location), c is related to the peak width, n is the number of peaks to fit.
Comparing the Performance of Fourier Decomposition and Wavelet Decomposition for Seismic Signal Analysis
Figure 3: (a) Collected Seismic signal by SASW. (b) A shifted Gaussian Window.
0.04 0.03 0.02
0.7 1 0.9 0.8
319
0.01 0 -0.01 -0.02 -0.03 -0.04
0.6 0.5
Multiply
0.4 0.3 0.2 0.1 0
0
20
40
60
80
100
120
0
500
1000
1500
2000
2500
Time(ms)
Time(ms)
(a)
Figure 4: A Gaussian slice or wave packet.
0.01
(b)
0.005
0
-0.005
-0.01
-0.015
0
20
40
60
80
100
120
140
Time(ms)
In this Study, the Fourier Decomposition is performed by the Complex Gaussian Wave in Fig. 5 which is regarded as a slice of the seismic signal.
320
Chik, Z, Islam, T, Rosyidi, S.A, Sanusi, H, Taha, M.R and Mustafa, M.M
Figure 5: Complex Gaussian wave as a slice of seismic signal.
1 0.8 0.6 0.4
Magnitude
0.2 0 -0.2 -0.4 -0.6 -0.8 -5
-4
-3
-2
-1
0
1
2
3
4
5
Time(s)
3.2. Fourier Analysis System
For the complex Fourier series (Letcher, 1992) representation of the signal x(t), the basis signal chosen are (7) φ n (t ) = e j 2π ( nf1t ) n = 0,±1,±2.......... ... ˆ The measure of the closeness of x(t ) to x(t) in the approximation interval is determined by the following approximation criteria used to choose the Xn’s.
(8) This equation produces the smallest area under the square of the error signal where the ˆ minimization is over possible x(t ) made by choices of Xn’s. We have to choose the Xn’s for a given set of basis signal to minimize the integral square error given by
t1 2 ˆ min ∫ [x(t ) − x(t )] dt t2
{
}
⎡N ⎤ ˆ ε N = ∫ [x(t ) − x(t )] dt = ∫ ⎢∑ X nφ n (t ) − x(t )⎥ dt (9) t1 t1 ⎣ n =1 ⎦ The property of a basis-signal set that produces finality of coefficients is mutual orthogonality of the basis signals over the approximation interval which is desirable to avoid coefficient recomputation. The real basis signals, φi (t ) are mutually orthogonal over the time interval t10. Moreover, a generalized Fourier series is a weighted sum of orthogonal basis signals that approximates a signal over the time interval t10 and the term of the last part of eq.(14) in brackets is squared, ε N is minimized if we choose the Xn’s so that the last term is zero. That is, the value of Xn that minimizes ε N is computed as 1 t2 Xn = (15) ∫ x(t )φn (t )dt t2 2 N
2
λn
t1
For the mutual orthogonal basis signal, the parameter λn can be depicted with the help of Eq. (7) and (10) as
λn =
t1 + T1
∫
t1
e
j 2 π n f1 t
e
− j 2 π n f1 t
dt =
t1 + T1
∫
t1
d t = T1 for all n.
(16)
So, according to Eq.(14), the integral square error that remains when N terms are used to ˆ approximate x(t ) is
N t2 ⎡ 1 t2 ⎤ ε N = ∫ x (t )dt − ∑ λn ⎢ ∫ x(t )φ n (t )dt ⎥ = ∫ x 2 (t )dt − ∑ λn X n2 (17) t1 t1 n =1 n =1 ⎣ λ n t1 ⎦ For the Fourier decomposition of the complex signal, the Fourier series coefficients and integral square error can be defined according to the Eq. (15) and (17). Basis on Eq. (2), the approximate signal of the series containing 2N+1 terms is given as t2 2 N 2
ˆ x N (t ) =
n=− N
∑
N
X n e j 2 π n f1 t
(18)
The total technique of the Fourier decomposition of the complex signal is shown below by Flow-chart in Fig.6.
322
Chik, Z, Islam, T, Rosyidi, S.A, Sanusi, H, Taha, M.R and Mustafa, M.M
Figure 6: Flow-Chart of Fourier decomposition System.
Chose basis signal , φn (t ) of the original signal x(t)
Find the Parameter
λn
Find the Fourier Series Coefficient Xn
×
n =− N
∑Xφ
N
n n
ˆ Approximate signal x(t )
The integral square error,
ˆ ε n = ∫ [ x(t ) − x(t )] dt
t1
t2
2
3.3. Wavelet Analysis System
For wavelet analysis, the seismic signal in Fig. 3a collected by SASW testing is used. The wavelet decomposition of this signal is performed using MATLAB according to the following Flow chart shown in Fig. 7.
Comparing the Performance of Fourier Decomposition and Wavelet Decomposition for Seismic Signal Analysis
Figure 7: Flow-Chart of Wavelet Decomposition System.
Embed the original signal x(t) Decomposition: 1. Choose a Wavelet 2. Choose a level, N 3. Compute the Wavelet Decomposition at level N.
323
Decomposition 1. Select a threshold for each level from 1 to N. 2. Apply soft thresholding to the detail coefficients.
Find the Threshold Details Coefficient
Compute the Wavelet Reconstruction using: 1. the original approximation coefficients of level N. 2. the modified detail coefficients of levels from 1 to N.
Perceive Noise-free signal at level N.
3.4. Result and Discussions
Fourier decomposition is very mathematical and not at all obvious. Fourier decomposition is important for three reasons. First, a wide variety of signals are inherently created from superimposed sinusoids. Audio signals are a good example of this. Fourier decomposition provides a direct analysis of the information contained in these types of signals. Second, linear systems respond to sinusoids in a unique way: a sinusoidal input always results in a sinusoidal output. In this approach, systems are characterized by how they change the amplitude and phase of sinusoids passing through them. Since an input signal can be decomposed into sinusoids, knowing how a system will react to sinusoids allows the output of the system to be found. Third, the Fourier decomposition is the basis for a broad and powerful area of mathematics called Fourier analysis, and the even more advanced Laplace and ztransforms. This work accentuates the basis on the Fourier series to obtain the minimum integral square error between the original signal and approximate signal for the received complex signal. Basis on the 0. Fig.8 signal analysis and the Eq. (10), if N ∞ for all signals for which a series exists, ε N shows the complex Fourier series expansion of complex signal where x(t) is Gaussian wave over the interval -5