- © 2013 by the Seismological Society of America
Online Material: Color versions of spectrogram figures; R and hht code installation instructions with examples.
The Fourier transform remains one of the most popular spectral methods in time‐series analysis, so much so that the word “spectrum” is virtually equivalent to “Fourier spectrum” (Huang et al., 2001). This method assumes that a time series extends from positive to negative infinity (stationarity) and consists of a linear superposition of sinusoids (linearity). However, geophysical signals are never stationary and are not necessarily linear. This results in a trade‐off between time and frequency resolution for nonstationary signals and the creation of spurious harmonics for nonlinear signals. We present an open‐source implementation of the Hilbert–Huang transform (HHT), an alternative spectral method designed to avoid the linearity and stationarity constraints of Fourier analysis. The HHT defines instantaneous frequency as the time derivative of phase, illuminating previously inaccessible spectral details in transient signals. Nonlinear signals become frequency modulations rather than a series of fitted sinusoids, eliminating artificial harmonics in the resulting spectrogram.
In this paper, we describe the HHT algorithm and present our recently‐developed hht package for the R programming language. This package includes routines for empirical mode decomposition (EMD), ensemble empirical mode decomposition (EEMD) and Hilbert spectral analysis. It also comes with high‐level plotting functions for easy and accurate visualization of the resulting waveforms and spectra. We demonstrate this code by applying it to three signals: a synthetic nonlinear waveform, a transient signal recorded at Deception Island volcano, Antarctica, and quasi‐harmonic tremor from Reventador volcano, Ecuador. The synthetic signal shows how the EMD method breaks complex time series into simpler modes. It also illustrates how the Hilbert transforms of nonlinear signals produce frequency oscillations rather than harmonics. The transient signal demonstrates the high‐time/frequency resolution of the HHT method. The volcanic‐tremor signal has high‐frequency harmonics in the Fourier spectrogram, which are not present in the Hilbert spectrogram. The EMD of the tremor signal also reveals unexpected transient events. Data and code for each analysis are included in the electronic supplement available with this paper.
THE HILBERT–HUANG TRANSFORM (HHT)
The HHT is a combination of the Hilbert transform and the EMD algorithm.
The Hilbert Transform
The Hilbert transform for a function x(t) is the convolution: (1)
The Hilbert transform of a periodic function produces a phase shift of π/2 for positive frequencies, so: (2)results in the analytic function: (3)
We can define the instantaneous amplitude a(t) by taking the magnitude of the real and imaginary components of y(t): (4)and the instantaneous phase by taking the arctangent of the real and imaginary components of y(t): (5)
The instantaneous frequency ω(t) is the derivative of phase with respect to time: (6)
This definition is only meaningful for monocomponent signals with a zero mean (Huang et al., 1998). Otherwise, the instantaneous frequency does not reflect the actual frequency content of the signal (Fig. 1).
Empirical Mode Decomposition (EMD)
The EMD method breaks a time series into a small number of monocomponent oscillatory modes called intrinsic mode functions (IMFs; Huang et al., 1998). Once the first IMF has been calculated, it is subtracted from the original signal to produce a residual. This residual is considered a new signal, and the EMD is applied again. The process repeats until the residual no longer contains any oscillations. Thus we can represent the original signal x(t) as a sum of IMFs plus a residual: (7)
The EMD is an adaptive stepwise filter, in which each successive IMF represents the highest‐frequency mode in the remainder. The final residual r is the data trend. Because each IMF is locally symmetric and has no riding waves, it has a meaningful Hilbert transform. The instantaneous frequency of each IMF can be plotted against time to produce a high‐resolution‐ensemble spectrogram of the original signal.
Ensemble Empirical Mode Decomposition (EEMD)
IMF sets can suffer from mode mixing and lack of uniqueness. Mode mixing consists of a single IMF containing signals of different time scales, or one signal scale residing on multiple IMFs (Huang and Wu, 2008). This can result in frequency aliasing in which IMFs transition from one scale to another. The EEMD is a noise‐assisted data analysis method designed to counteract this problem. It repeatedly adds uniform white noise to the signal, performs EMD, and averages the IMFs of each trial together (Wu and Huang, 2009). As the number of trials increases, the average of noise‐perturbed copies of a signal approach the true signal. The EEMD method greatly reduces mode mixing and thus represents a significant improvement over EMD (Huang and Wu, 2008). Because each EMD trial is independent, the EEMD can be rapidly calculated on a computer cluster.
THE HHT IN SEISMOLOGY
Earthquakes often have rapidly changing spectral characteristics, and these details can be lost when seismograms are windowed during Fourier spectral analysis. However, Hilbert spectral analysis defines frequency for every single sample point, permitting the identification and examination of transient events in complicated seismograms. The sharp frequency resolution also bodes well for the development of more precise methods of calculating dispersion curves; see Chen et al. (2002) for a preliminary analysis using synthetic seismograms. The EMD method also opens the door for new and improved signal analysis techniques. For example, the IMF set can reveal frequency gliding and modulations within a single mode, details which could be lost through simple band‐pass filtering. Higher dimensional EMD methods might be able to extract polarity or phase‐arrival data from triaxial instruments (see Rehman and Mandic, 2010), but this has not yet been tested on seismic data.
A HHT analysis of the 1999 Chi‐Chi earthquake revealed high‐amplitude, low‐frequency energy, which the Fourier and wavelet transform did not show (Huang et al., 2001). Another study of the same earthquake discovered high‐amplitude, low‐frequency ground‐acceleration components in Fourier spectra of IMFs generated by the EMD method (Loh et al., 2000). Zhang, Ma, and Hartzell (2003) used the EMD method to quantify source and rupture‐propagation information from seismograms recorded during the 1994 Northridge earthquake, producing results consistent with previous studies. Zhang et al. (2003) found that the HHT can extract low‐frequency pulse‐like signals and other nonstationary features in earthquake records. Zhang (2006) compared Fourier‐based and HHT‐based site‐amplification factors for the 2001 Nisqually earthquake and found that the two methods were equivalent for linear sites. However, when site nonlinearity was high, the HHT site‐amplification factor was larger than the Fourier site‐amplification factor in the low‐ to medium‐frequency range. The increased fidelity of HHT with respect to Fourier at low frequencies has important consequences for long‐period structure design, as it is precisely these frequencies that cause the most damage.
THE HHT PACKAGE IN R
R is an open‐source programming language designed for statistical analysis and scientific computation (R Core Team, 2012). Binaries for basic R installations are available for Windows, MacOS, and many UNIX systems. In addition, R can be compiled from source. These installation options, as well as extensive documentation and additional tools, are available at http://www.r-project.org/ (last accessed April 2013). A wide variety of other resources are available for learning R, ranging in difficulty from beginning level (R for Dummies; de Vries and Meys, 2012) to ones that assume some programming experience (An Introduction to R, Venables et al., 2013).
Bundles of user‐submitted, quality‐controlled open‐source R code called packages are available for installation as well. These packages greatly extend the core functionality in R. The analyses and figures presented in this paper were prepared using the hht package (Bowman and Lees, 2013) The hht package can be installed directly from R by typing the following command into the R interpreter:
Because hht depends on the EMD package (Kim and Oh, 2013) and the fields package (Furrer et al., 2013), these are automatically included in the installation process. Alternatively, the hht package can be downloaded from the R package repository at http://cran.r-project.org/ (last accessed April 2013), in which case EMD and fields must be downloaded and installed separately. MacOS users must download a FORTRAN compiler from cran.r-project.org/bin/macosx/tools (last accessed April 2013) prior to installing hht in either case; however, UNIX and Windows users should not require additional resources. Once downloaded, the package is quite simple to use. Here we demonstrate how to generate Figure 2 in this paper.
Data and R code for each figure in this paper are included in (see supplement). There do not appear to be any MATLAB equivalents to the EMD or hht R packages, although a couple of users have uploaded simple routines to the MATLAB File Exchange website.
EXAMPLES OF NONLINEAR AND NONSTATIONARY TIME SERIES
Synthetic Nonlinear Wave
The first example consists of a low‐frequency sinusoidal carrier wave, a nonlinear Stokian wave modified from equation 8.5 in Huang et al. (1998), and random noise. Our signal is of the form (8)in which is Gaussian random noise with a mean of 0 and a standard deviation of 0.1. The Stokian component of this equation creates a frequency‐modulated signal centered at 3 Hz superimposed on a constant‐frequency signal at 1 Hz. The sinusoidal carrier wave has a frequency of 0.05 Hz. The EMD method returns eight IMFs (Fig. 2). IMFs 1 through 3 are high‐frequency noise from . IMF 4 contains the 3 Hz frequency‐modulated Stokian signal. The frequency‐modulated signal also appears in IMF 3 and 5 from time to time. IMF 5 contains the 1 Hz constant‐frequency Stokian component, but this signal sometimes switches to IMF 6. IMF 8 represents the 0.05 Hz carrier wave. Ideally, the high‐ and low‐frequency Stokian‐wave components should be returned in single IMFs, but both of them switch back and forth between multiple IMFs. This intermittency can violate local symmetry, introducing severe aliasing in the Hilbert transforms of IMFs 3 through 6.
We applied the EEMD method to reduce the effects of mode mixing. The resulting IMF set after 1000 trials at a noise amplitude of 0.1 (Fig. 3) has much less mode mixing than the raw IMF set (Fig. 2). The 3 Hz frequency‐modulated Stokian component now lies solely in IMF 4, the 1 Hz Stokian component lies in IMF 5, and the 0.05 Hz carrier wave lies in IMF 8. Although the sixth IMF still contains some intermittency, the EEMD has clearly improved the quality of the original IMF set. The Hilbert spectrum consists of all 1000 trials plotted in the time/frequency domain. We see a linear energy band at 0.05 Hz, another linear energy band at 1 Hz, and an oscillating energy band centered at 3 Hz, corresponding to the three components of the signal as described above (Fig. 4).
Our next example is a transient signal recorded by an ocean‐bottom seismometer in the flooded caldera of Deception Island volcano. The signal consists of a low‐amplitude first arrival followed by a much higher‐amplitude second arrival approximately one second later. There are three main IMFs: a high‐frequency high‐amplitude component, a low‐frequency high‐amplitude component, and a low‐frequency low‐amplitude component (Fig. 5). The high‐frequency high‐amplitude component (IMF 1) has noticeable frequency modulation. These three IMFs demonstrate the ability of the EMD method to reveal hidden time scales in the data, including frequency modulation that might have been lost had a band‐pass filter been used instead.
The Hilbert spectrogram of of the event shows a a high‐frequency phase between 10 and 20 Hz and a low‐frequency phase at approximately 4 Hz (Fig. 6). There is also a lower‐energy lower‐frequency phase present in the first second of the signal. All three components have a decreasing frequency trend in the first half second of the signal, but the 4 Hz component rises in frequency over the following 1.5 s. The complicated structure in the high‐frequency gliding phase as well as the highly scattered frequencies at the signal onset are probably a result of the sudden changes in frequency and amplitude in the signal and may not reflect reality (Huang et al., 1998).
Volcanic tremor is a long‐duration seismic (and sometimes seismoacoustic) signal produced by fluid flow in active volcanoes. The Fourier spectrogram of volcanic tremor often includes multiple harmonics and frequency gliding. We selected an example of tremor recorded at Reventador volcano (Lees et al., 2008) for analysis with the HHT. This tremor is ideal because it includes periods of harmonic and nonharmonic behavior as well as frequency gliding. Clear harmonics are visible from approximately 30 to 160 s and again from 350 to 600 s in the Fourier spectrogram. The harmonics show frequency gliding from approximately 30 to 120 s, with another episode around 450 s (Fig. 7).
The averaged IMF set produced by EEMD shows that the tremor is very complex (Fig. 8). Most energy lies in IMFs 3, 5, and 6, with the other IMFs being either very low amplitude or comprised of combinations of the three main IMFs. For example, the averaged IMF 4 is clearly a mixture of IMF 3 and 5, because it contains high‐frequency riding waves on top of lower‐frequency carrier waves, violating the definition of an IMF. However, this IMF does contain pure signal, particularly near 60 s. The three well‐resolved IMFs contain subtle structures, such as high‐frequency spindle‐shaped sequences and bursts in IMF 3, an amplitude modulated signal in IMF 5, and a lower‐frequency episodic signal in IMF 6. Although the presence of episodic high frequencies is evident in the original signal, the pulsating structure of this high‐frequency component is only revealed through EMD.
Like the Fourier spectrogram, the Hilbert spectrogram portrays this tremor signal as a series of frequency bands that change with time (Fig. 9). However, the Hilbert spectrogram shows little energy above 3 Hz and almost none above 4 Hz, whereas the Fourier spectrogram shows harmonics up to 6 Hz. The Fourier spectrogram does not show much energy below 1 Hz, whereas the Hilbert spectrogram does show significant low‐frequency energy between 300 and 600 s. The Hilbert spectrogram also shows fewer frequency modes than the Fourier spectrogram. There are two frequency bands from 40 to 200 s, as opposed to three main harmonics and over ten minor harmonics in the Fourier spectrogram over the same time interval. Also, the Hilbert frequency bands are located at lower frequencies; for example, at 100 s the Hilbert modes are at 1 and 1.75 Hz, whereas the Fourier modes are at about 1.75, 3.25, and 5 Hz. The Hilbert spectrogram bifurcates at about 350 s, with three low‐frequency modes and one higher‐frequency mode. The higher‐frequency mode is noisier but appears to glide over a greater frequency range than the other modes. The Fourier spectrogram also shows four modes, but they are higher frequency and also contain more energy in the two lower‐frequency modes, whereas the Hilbert spectrogram shows the most energy in the two higher‐frequency modes.
The Hilbert spectrogram of the tremor signal between 454 and 470 s shows three frequency bands: a low‐energy one at about 0.75 Hz, a high‐energy one at about 1.3 Hz, and a strongly modulated one at 3 Hz (Fig. 10). The lower‐frequency band dies out at 458 s, and the modulated signal grows in amplitude at that time. The middle‐frequency band loses energy just before 460 s, and the modulated signal grows even more energetic. Both lower bands regain energy at about 466 s whereas the modulated signal loses amplitude. The Hilbert spectrogram clearly portrays the complicated interaction of three different signals over a very short time frame.
The HHT combines the precision of the Hilbert transform with an adaptive signal decomposition method. Although any signal can be represented as a linear combination of arbitrary‐basis functions that span the data space, choosing the wrong basis function can greatly increase the number of terms required to fit the time series. For example, the Fourier method uses sinusoidal‐basis functions, resulting in the production of numerous (possibly infinite) harmonics when a nonsinusoidal signal is processed. Because the EMD method chooses the basis functions from the properties of each individual signal, it greatly reduces the number of terms required to represent the signal in both the time and frequency domains. Not only does this minimize the number of parameters that need to be modeled in order to reproduce the signal, it also prevents spectral energy from bleeding into nonphysical harmonics, which can obscure the true frequency content of the time series.
The HHT also has disadvantages compared to other spectral analysis methods. The EMD is computationally expensive, especially when the time series is long, has a large frequency distribution, and/or has a high sample rate. Random noise can perturb IMF sets, an effect that the EEMD method can generally correct at the expense of even more computer time. The HHT method will generally return the Fourier components of a simple linear signal. However, because the EMD method relies on extrema to separate IMFs from the signal, it can fail to split a low‐amplitude Fourier component from a high‐amplitude Fourier component if the higher‐frequency mode does not produce extrema. This may indicate that the signal is nonlinear when in fact it is a superposition of two sinusoids. Dätig and Schlurmann (2004) provide a rigorous analysis of this phenomenon. They also find that improper spline fitting can have a deleterious effect on the EMD. Kijewski‐Correa and Kareem (2006) verify the inability of the EMD to handle closely spaced waves, and they point out instances in which amplitude modulation in the time series produces nonphysical frequency modulation in the Hilbert spectrogram. They also find that the wavelet spectrogram can match the precision of the Hilbert spectrogram if the appropriate wavelet is chosen.
This research was supported under National Science Foundation Grants DGE‐1144081 and 5‐37139 as well as the UNC Martin Fund. We thank Rick Aster for his detailed review and insightful comments.