Tary 2014 Spectral

Reviews of Geophysics
REVIEW ARTICLE Spectral estimationWhat is new? What is next?

10.1002/2014RG000461
Jean Baptiste Tary1,2 , Roberto Henry Herrera1 , Jiajun Han1 , and Mirko van der Baan1
Key Points: 1 Department of Physics, University of Alberta, Edmonton, Alberta, Canada, 2 Now at Department of Meteorology and
Review of (non)conventional
Geophysics, University of Vienna, Vienna, Austria
techniques for time-
frequency representations
Resolution and performance of
nine techniques are tested on ve Abstract Spectral estimation, and corresponding time-frequency representation for nonstationary
benchmark signals signals, is a cornerstone in geophysical signal processing and interpretation. The last 1015 years have seen
Recently developed techniques show
the development of many new high-resolution decompositions that are often fundamentally dierent from
clear improvements in most cases
Fourier and wavelet transforms. These conventional techniques, like the short-time Fourier transform and
the continuous wavelet transform, show some limitations in terms of resolution (localization) due to the
Supporting Information:
trade-o between time and frequency localizations and smearing due to the nite size of the time series
Benchmark signals
of their template. Well-known techniques, like autoregressive methods and basis pursuit, and recently
developed techniques, such as empirical mode decomposition and the synchrosqueezing transform, can
Correspondence to:
achieve higher time-frequency localization due to reduced spectral smearing and leakage. We rst review
J. B. Tary,
jean-baptiste.tary@univie.ac.at the theory of various established and novel techniques, pointing out their assumptions, adaptability,
and expected time-frequency localization. We illustrate their performances on a provided collection of
benchmark signals, including a laughing voice, a volcano tremor, a microseismic event, and a global
Citation:
Tary, J. B., R. H. Herrera, J. Han, earthquake, with the intention to provide a fair comparison of the pros and cons of each method. Finally,
and M. van der Baan (2014), their outcomes are discussed and possible avenues for improvements are proposed.
Spectral estimationWhat is new?
What is next?, Rev. Geophys., 52,
doi:10.1002/2014RG000461.
1. Introduction
Received 9 MAY 2014 Spectral or period estimation is a fundamental analysis tool that has a wide range of applications. Accu-
Accepted 18 SEP 2014 rate and precise analysis of nonstationary spectral variations is a longstanding problem aiming at revealing
Accepted article online 23 SEP 2014
signal characteristics such as any underlying periodicities. The discrete Fourier transform (DFT) and its fast
implementation, the fast Fourier transform (FFT) [Cooley and Tukey, 1965; Cooley et al., 1969], and variants
have been well studied [e.g., Hinich and Clay, 1968; Brigham, 1988]. They provide an ecient way to estimate
the frequency content of a discrete and nite time series. The main application of the DFT to time-varying
spectra is the spectrogram, also called short-time Fourier transform (STFT). Likewise, wavelet transforms and
variants are now well-established and many review articles exist [e.g., Kumar and Foufoula-Georgiou, 1997;
Mallat, 2008].
All these methods are bound by the Heisenberg/Gabor uncertainty principle [Gabor, 1946] with a trade-o
between time and frequency localizations [Reine et al., 2009]. In this review, resolution refers to the ability
to resolve two signals close in time or frequency. It includes both localization, which corresponds to the
intrinsic ability of a technique to separate signals in the tf domain, and sampling (see section 2.2 for more
explanations). Signal windowing leads to smearing, which corresponds to the widening of the main lobe
around its central frequency, as well as sidelobe leakage which introduces spurious frequencies [Hall, 2006]
(Figure 1). Lately, various new transforms have been developed to ameliorate these issues such as empir-
ical mode decomposition (EMD) coupled with the Hilbert-Huang spectrum [Huang et al., 1998] and the
synchrosqueezing transform (SST) [Daubechies et al., 2011].
All the above-mentioned methods can be considered as nonparametric, meaning that they do not assume
any particular stationary structure prior to time-frequency (tf ) decomposition of the signal. On the other
hand, parametric methods drew a lot of attention in the 1980s and 1990s as they do not assume zero-valued
or periodic data outside the data window and hence are less restricted by smearing and leakage. Autore-
gressive and autoregressive moving-average models for example provide correct and accurate spectral
estimation if the assumed model is appropriate for the investigated time series [Makhoul, 1975; Ulrych and
Bishop, 1975; Kay and Marple, 1981].
Here we provide a comprehensive survey of the current state of parametric and nonparametric tf trans-
forms including their latest generation. Along the way, we dene the time and frequency localization
TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 1

Reviews of Geophysics 10.1002/2014RG000461
150
a) c)
Norm. amplitude
n
tio
Frequency (Hz)
Smearing
lu
100
so
re
t-f
e
bl
Leakage
ia
ar
:v
50
ST
olution
ed t-f res
0 5 10 15 STFT: fix
0
Improper 0 2 4 6 8 10
b)
freq. sampling Time (s)
Norm. amplitude
FFT, 1024 points

FFT, 50 points
AR, p=20
AR, p=10
0 5 10 15
Frequency (Hz)
Figure 1. Illustration of (a) spectral smearing and sidelobe leakage and (b) frequency sampling. An untapered signal
of 2 s consisting of two sinusoids at 5 and 10 Hz is analyzed using the fast Fourier transform (FFT) and an autoregres-
sive (AR) method (Burg algorithm) [Burg, 1972] with dierent parameters. Contrary to the autoregressive spectrum, the
Fourier spectrum with 1024 frequency points shows spectral smearing and sidelobe leakage. The Fourier spectrum with
only 50 frequency points does not show correctly the frequencies of the two sinusoids. The inuence of the autore-
gressive order p on tf localization and accurate frequency determination for autoregressive methods is shown by the
two autoregressive spectra with dierent orders. (c) Two linear chirps of 10 s are analyzed using the short-time Fourier
transform (STFT) and S transform (ST) to illustrate the concepts of xed and variable tf localization, respectively. The
instantaneous frequencies of both chirps are indicated by the dashed-dotted black lines. The short-time Fourier trans-
form, computed using a Hann taper of 0.3 s with 90% overlap, is applied to the chirp with a frequency sweep from 0 to
30 Hz, while the S transform is applied to the chirp with a frequency sweep from 0 to 160 Hz.
of each technique, as well as the dierence between spectral resolution, localization, and sample den-
sity. We present conventional nonparametric transforms used in geophysical signal processing, i.e., the
short-time Fourier transform, the continuous wavelet transform (CWT), and the S transform (ST). We omit
the Wigner-Ville distribution because of its main disadvantage: cross-terms without physical meanings
that can sometimes mask the true spectral content [Cohen, 1989]. We then describe techniques that were
developed recently, namely, the synchrosqueezing transform [Daubechies et al., 2011] and empirical mode
decomposition and variations [Huang et al., 1998; Flandrin et al., 2004; Han and van der Baan, 2013]. Finally,
new developments on the estimation of autoregressive models are presented. Within the geoscience com-
munity, Kalman lter techniques are mainly employed for time series modeling and forecasting [e.g., Bertino
et al., 2003; Segall et al., 2006; Dikpati et al., 2014] and not for spectral estimation.
We describe rst their rationale and underlying principles. We then analyze their respective advantages
and inconveniences using synthetic, speech, and seismic data to illustrate their distinctive features as well
as their complementarities. Even though the presented spectral methods are mainly illustrated on seismic
data, most were applied to other elds of earth sciences such as climate time series [Walker, 1931; Torrence
and Compo, 1998; Ghil et al., 2002; Salisbury and Wimbush, 2002], hydrology [Labat, 2005], ecology [Cazelles
et al., 2008], geology [Bolton et al., 1995; Kravchinsky et al., 2013; Thakur et al., 2013; Reager et al., 2014],
oceanography [Farge, 1992], or solar physics [Miyahara et al., 2006; Barnhart and Eichinger, 2011]. Abbrevia-
tions and template for each spectral estimation techniques are summarized in Table 1.
2. Methods
2.1. Short-Time Fourier Transform
The Fourier transform is a measure of the similarity of a signal with a family basis formed by sines and
cosines. This can be expressed as the inner product of a signal s(t) with a template (t), i.e., s(t), (t) =

s(t) (t)dt, where stands for the complex conjugate [Gao et al., 2010]. When this template is made of

Table 1. Abbreviations, Parametric or Not, and Templates of the tf Techniques Presented in This Review
Techniques Abbreviation Parametric Template

Short-time Fourier transform STFT No Sines and cosines
Continous wavelet transform CWT No Stretched or squeezed, and translated versions of the mother wavelet
S transform ST No Sines and cosines
Synchrosqueezing transform SST No Same as CWT
Matching pursuit MP No User dened dictionary of elementary waveforms
Basis pursuit BP No User dened dictionary of elementary waveforms
Short-time autoregressive method STAR Yes N/A
Kalman smoother KS Yes N/A
Empirical mode decomposition EMD No N/A
Ensemble empirical mode decomposition EEMD No N/A
Complete ensemble empirical mode decomposition CEEMD No N/A
sines and cosines, by Eulers equation for the complex-value cosine (t) = ei2ft , the inner product reduces
to the Fourier transform

FT{s(t)} = s(t), ei2ft = s(t)ei2ft dt, (1)

where t is time, f is frequency, and i = 1.
The Fourier transform is applied to a signal in its entirety and thus assumes stationarity. The short-time
Fourier transform is the most popular extension to handle nonstationary time signals. Its easy
implementation as well as its good performance explain why this is the most widely used technique for
tf representation. Indeed, this is simply the Fourier transform of successive, usually overlapping, win-
dows of the signal, each frequency distribution being associated with the central time of each window. The
continuous time short-time Fourier transform can be expressed as
( )
SF , f = s(t)w (t ) ei2ft dt, (2)

corresponding to the Fourier transform of the signal s(t) tapered by the function w(t), which is usually a
Hann or Gaussian function, and where is the time delay to the center of the taper. It is easy to see that
equation (2) is a special case of equation (1) where we use a tapered window w(t) to select which part of s(t)
will be transformed, i.e., s(t) w(t), ei2ft .
2.2. Time-Frequency Resolution

At this point, it is important to distinguish between time and frequency resolution, localization and sam-
pling (Figure 1). The sample rate in the frequency domain is inversely proportional to the size T of the taper
w(t) in the time domain (f = 1T ). There is always the option of zero-padding in the time domain to
decrease f , improving the readability of the spectrum by increasing the frequency sampling. But well sam-
pled does not mean high resolution [Kay and Marple, 1981]. By resolution, we usually mean the ability to
distinguish two spectral peaks in a similar way as we do in the time domain to distinguish between two
events. Frequency resolution then depends not only on frequency sampling but also on the localization
power of each technique. Spectral localization corresponds to the intrinsic ability of each technique to focus
individual tf components in the computed tf representations, reducing spectral smearing and sidelobe
leakage (Figure 1).
The untapered FFT of a signal is the convolution of a boxcar function with the original time series. The band-
width of the main lobe of the frequency spectrum of the boxcar function is inversely proportional to the
length of the boxcar (2T ). The sinc function corresponding to its Fourier transform introduces sidelobes
that can mask low-amplitude peaks. The result is a peak smeared around its peak frequency with decaying
sidelobes. These eects can be reduced by using tapers such as Hann, Gaussian, or Kaiser-Bessel functions
[Harris, 1978]. The attenuation of sidelobes is however associated with a widening of the main lobe [Hall,
2006]. Spectral smearing can be reduced by enlarging the taper length, yet this reduces the time localization
by imposing increased stationarity. Conversely, reducing the taper length improves the time localization at
the expense of enhancing spectral smearing, thus lowering spectral localization.

This is known as the Gabor uncertainty principle [Gabor, 1946] by analogy with the Heisenberg uncertainty
principle in quantum mechanics [Kumar and Foufoula-Georgiou, 1997]. The Gabor uncertainty principle tells
us that if we dene the standard deviation in time and frequency localization by t and f , respectively,
then there is a theoretical limit in localization given by [Hall, 2006]
1
t f . (3)
4
For signals s(t) with a duration longer than the taper w(t) support, the xed length for this taper in the
short-time Fourier transform sets the recoverable time and frequency localizations. When the signals dura-
tion is shorter than the support of the taper, then the localization is limited by the signal support itself and
not the taper.
Two alternative strategies were developed in the 1980s and 1990s to overcome the undesirable eects
introduced by xed taper lengths for the analysis window w(t) and resulting spectral smearing. Various tf
transforms such as the continuous wavelet transform use variable window lengths inversely proportional
to the central target frequency (Figure 1). Hence, long windows used for the low frequencies have good
frequency but limited time localization, whereas the short windows used for high frequencies have good
time but limited frequency localization [Kumar and Foufoula-Georgiou, 1997]. The underlying concept is that
for short-duration signals, timing errors are relatively more important than errors in frequency content and
vice-versa for long-duration signals.
The second strategy is called frequency reassignment and attempts to reduce spectral smearing [Auger
and Flandrin, 1995; Han et al., 2014]. Temporal and spectral smearing indicates that each tf amplitude is
a weighted average of neighboring points in the tf plane. This implies that a nonzero amplitude can be
retrieved even if the true signal has no component at this tf pair. The reassignment method computes
the sphere of inuence of each convolution kernel associated with an analysis window and reallocates the
energy in the tf plane to its center of gravity in the time and frequency domains, improving the readabil-
ity of the tf picture. The reassignment method has also been implemented for other methods such as
Wigner-Ville distributions or the continuous wavelet transform [Auger and Flandrin, 1995; Auger et al., 2013].
2.3. Continuous Wavelet Transform
The wavelet transform is a well-developed mathematical tool that uses a variable length for the analysis
window in equation (2), in contrast to the short-time Fourier transforms xed length window. Many excel-
lent books and tutorials describe the mathematical foundations of the wavelet transform [Daubechies,
1992; Mallat, 2008; Rioul and Vetterli, 1991; Torrence and Compo, 1998]. In this review, for the sake of simplic-
ity and space, we review briey the rationale and underlying principles to help in the understanding and
application of wavelets in geophysical problems.
In the continuous wavelet transform, the analysis taper w(t) is replaced by a mother wavelet that can be
shifted in time and stretched by a scale a, yielding equation (4).
( )
1 t
Ws (a, ) = s(t) dt. (4)

a a
In other words, w(t) = 1 ( t

a
), where the multiplication factor
1
ensures energy normalization, such
a a
that the wavelet always has unit energy at all scales:

|a, (t)|2 dt = |(t)|2 dt = 1, (5)

where (t) is the original mother wavelet and a, (t) represents the generated family of translated and
scaled wavelets.
A comparison of equations (2) and (4) shows that the continuous wavelet transform shares a similar logic
to the short-time Fourier transform; it represents the inner product of the signal with a time shifted and
stretched wavelet s(t), 1 ( t
a
). The continuous wavelet transform is simply the cross-correlation of
a
the signal s(t) with a number of wavelets that are scaled and translated versions of the original mother
wavelet. Ws (a, ) are the coecients representing a concentrated tf picture. The actual frequency f is built
from the central frequency of the mother wavelet f0 at each scale a with this relation f = f0 a [Mallat

and Zhang, 1993; Reine et al., 2009]. Practical implementations involve the discretization of the scale axis
[Thakur et al., 2013].
A mother wavelet (the starting wavelet) is chosen with some specic characteristics in order to improve
the signal representation. In all cases, the choice of the wavelet is neither unique nor arbitrary [Kumar and
Foufoula-Georgiou, 1997]. This template has to satisfy the admissibility condition [Farge, 1992; Daubechies,
1992], given by
+
C = 2 ||1 |()|
2
d < , (6)
0
where ()
is the Fourier transform of the mother wavelet (t) and is angular frequency. The condition
in equation (6) can only be satised if the frequency domain wavelet at zero frequency is zero, i.e., (0)
= 0.
In other words, the time domain wavelet has zero mean: (t)dt = 0.
The wavelet spectrum () should be continuously dierentiable, which is validated if the time domain

wavelet has a sucient time decay [Mallat, 2008]: (1 + |t|)|(t)|dt < +. This is called the regularity
condition. Thus, we have a small waveform that oscillates around zero with equal positive and negative
area (i.e., zero mean) and has a short duration, creating a wavelet. The importance of having a wavelet
that satises the admissibility condition is that the original signal, s(t), can be reconstructed by the inverse
wavelet transform.
2.4. S Transform
The S transform is halfway between the short-time Fourier transform and the continuous wavelet transform.
It uses a Gaussian template which is time-shifted by and stretched or squeezed by a factor inversely pro-
portional to the linear frequency f [Stockwell et al., 1996]. This is comparable to the mother wavelet used by
the continuous wavelet transform, leading Stockwell et al. [1996] to describe the S transform as a continuous
wavelet transform with a phase shift. The frequency-dependent Gaussian taper is dened by
2
t2
w(t) = e 2 , (7)
k
where = |f |
.
Incorporating and normalizing the Gaussian taper function lead to the expression for the
Gaussian template
( ) |f | 2 2
t f
w t, f = e 2k2 . (8)
k 2
The parameter k, by controlling the width of the Gaussian taper, can also be tuned to obtain better fre-
quency localization at the expense of reduced time localization. Substituting the Gaussian template for
(t)2 f 2
|f |
taper w(t) in equation (2), the S transform becomes the inner product s(t) e 2k2 , ei2ft , yielding
k 2
( ) |f | (t)2 f 2
SST , f = s(t)e 2k2 ei2ft dt. (9)

k 2
Like the continuous wavelet transform, the tf localization of the S transform is variable (Figure 1). The vari-
able tf localization of the S transform and continuous wavelet transform enables the analysis of signals
with frequency components separated by few orders of magnitude (e.g., 1 and 1000 Hz).
Unlike the continuous wavelet transform, the S transform keeps the properties of the Fourier transform such
as a uniform frequency sampling, spectral amplitudes independent of frequency owing to the normaliza-
tion of the Gaussian taper to unit area in equation (9), and retains the original signal phase [Stockwell, 2007].
On the other hand, the continuous wavelet transform is characterized by an exponential scale sampling (i.e.,
more scales are allocated on lower portions than higher portions of the spectrum, producing an equal num-
ber of scales per octave instead of a linear frequency sampling), locally dened phase, and show higher
amplitudes at lower frequencies.

2.5. Synchrosqueezing Transform

The synchrosqueezing transform is a derivation of the continuous wavelet transform empowered by a reas-
signment step [Daubechies et al., 2011]. It concentrates the frequency content around the instantaneous
frequencies in the wavelet domain. The synchrosqueezing transform assumes that the signal s(t) is the sum
of individual time-varying harmonic components:

K
s(t) = Ak (t) cos(k (t)) + (t), (10)
k=1
1 d
where Ak (t) is the instantaneous amplitude and fk (t) = 2 (t) is the instantaneous frequency of the signal
dt k
component k, derived from the instantaneous phase k (t). The additive noise (t) includes the contribution
of environmental and acquisition sources, and K stands for the maximum number of components in the
signal. This is a mathematical simplication, but this will enable signal separation and analysis.
The rationale behind this technique is that concentrating the tf map into the most representative instanta-
neous frequencies will decrease smearing while still allowing reconstruction. With the continuous wavelet
transform, the variable length of the mother wavelet leads to a more exible trade-o between time and
frequency localization but does not prevent spectral smearing in the tf plane. Daubechies et al. [2011]
observed that there is a physical limit to reduce the smearing eect in the tf representation using the
continuous wavelet transform. This smearing mainly occurs along scale axis; thus, if we want to get fre-
quency variations with time, we need to compute the instantaneous frequencies also in the wavelet domain.
The angular instantaneous frequency is the rate of change of the time-dependent phase ((t) = d(t) dt
)
[Taner et al., 1979]; likewise in the wavelet domain, the instantaneous frequency s (a, ) is computed as the
derivative of the wavelet transform at any point (a, ) with respect to , for all Ws (a, ) 0:
i Ws (a, )
s (a, ) = . (11)
2Ws (a, )
The instantaneous frequencies are the ridges in the tf representation [Auger et al., 2013]. To decrease the
smearing, we have to squeeze the frequencies around these ridges. This is to map the information from the
time-scale plane to the tf plane. Every point (a, ) is converted to (s (a, ), ), and this operation is called
synchrosqueezing [Daubechies et al., 2011]. Since the scale a and time are discrete values, we can have
a scaling step ak = ak1 ak for any ak where Ws (a, ) is computed. Likewise, when mapping from the
time-scale plane to the tf plane (a, ) ((a, ), ), the synchrosqueezing transform Ts (, ) is determined
only at the centers l of the frequency range [l 2, l + 2], with = l l1 :
1
Ts (l , ) = Ws (ak , )a32 ak . (12)
a
k
where the scale ak is bounded by |(ak , ) l | 2. This means that a small frequency band centered
at l and less than half of the bandwidth is used to reconstruct the concentrated frequency map at time .
The new tf representation of the signal Ts (l , ) is synchrosqueezed along the frequency (or scale) axis only
[Li and Liang, 2012].
The discretized version of Ts (l , ) in equation (12) is represented by T s (wl , tm ), where tm is the discrete time
tm = t0 + mt with t the sampling interval, m = 0, ..., n 1 and n is the total number of samples in the
discrete signal sm . The frequency localization in the synchrosqueezing transform is also limited by the Gabor
uncertainty principle. The Nyquist theorem gives the upper frequency limit fmax = 1(2t) and the mini-
mum frequency is fmin = 1(nt), with n the signal length. With these two bounds, the synchrosqueezing
transform discretizes the scale axis to compute the frequency division step l . More special considerations
are described in Thakur et al. [2013].
Summing up the previous steps, the synchrosqueezing transform assumes that the signal is a superposition
of nonstationary monochromatic wavelets and can be eciently decomposed by the continuous wavelet
transform, followed by the computation of the instantaneous frequencies plus a reassignment step to con-
centrate the energy around the ridges. An extended explanation with applications to seismic signals can be
found in Herrera et al. [2014].

The reconstruction of the individual components sk from the discrete synchrosqueezed transform T s is then
the inverse continuous wavelet transform over a small frequency band lLk (tm ) around the kth component:
( )

sk (tm ) = 2C1 T s (wl , tm ) , (13)
lLk (tm )
where C is a constant dependent on the selected wavelet. As we take the real part of the discrete syn-
chrosqueezing transform in that band, we recover the real component sk . In this review, we follow Thakur
et al. [2013] where the reconstruction is done by a standard least-squares ridge extraction method.
2.6. Matching Pursuit
In the matching pursuit (MP) algorithm, the signal is decomposed into waveforms selected from a dictionary
of tf atoms [Mallat and Zhang, 1993]. Time-frequency atoms are discrete-time elementary waveforms pop-
ulating the dictionary. This technique assumes that the signal is a linear combination of these atoms from
a dictionary that is not necessarily built with modulated sine and cosine waves atoms. Matching pursuit
uses a recognition step to identify and subsequently remove the most prominent atom via a deationary
approach and then repeats the recognition step [Chen et al., 2001]. The subtraction step prevents smearing
and leakage, creating highly localized tf decompositions.
Generally, the tf atoms are the dilations, translations, and modulations of a single window function:
( )
1 t i(t)
(t) = e , (14)
a a
where = (a, , ), i.e., the three transformations scale, modulation, and translation, respectively. Thus, the
recognition step of matching pursuit comprises the inner product of the signal s(t) and then the residual
signal R(n) s which is the original signal minus any subtracted wavelets, with a time shifted, stretched, and
modulated wavelet, i.e., s(t), (t) or R(n) s, (t).
Mallat and Zhang [1993] proposed to build a dictionary D, with an entire set of waveforms , in order to
best match the signal structures. Matching pursuit uses a greedy algorithm, which means that at each iter-
ation it chooses a waveform from the dictionary that is best adapted to the signal segment under analysis
and repeats this operation until a predened stopping criterion is reached. This produces a sparse represen-
tation of the signal, s(t), as the weighted sum of chosen waveforms ( ), picked up from the dictionary D
[Wang, 2007]:

N1
s(t) = an n (t) + R(n) s, (15)
n=0
where N is the total number of iterations, an is the amplitude of the wavelet n , and R(n) s is the residual at
iteration n, with R0 s = s.
Dictionaries are built with dierent waveforms D = { }D . The user determines the extent of the dictio-
nary and the type of waveforms used. A large dictionary allows for a more complete signal representation
but at the expense of signicantly longer computation times. If the waveforms are sines and cosines, then
it is a Fourier dictionary and the representation is analogous to a sparse Fourier transform. If a wavelet dic-
tionary is built from a mother wavelet, then the representation becomes analogous to a sparse continuous
wavelet transform. On the other hand, a tf dictionary, also known as Gabor dictionary [Mallat and Zhang,
1993], produces a decomposition of the signal s(t) into tf atoms. Computing the Wigner distribution of
each atom leads to the tf representation [Mallat and Zhang, 1993]. Matching pursuit may lead to similar
time-frequency representations compared with the previous techniques, but its deationary approach pre-
vents spectral leakage, thereby creating a much sparser representation. Unlike the previous methods, the
matching pursuit algorithm is not exploring the complete tf plane but is mapping tf atoms located at the
positions of maximum residual reduction [Addison, 2002].
2.7. Basis Pursuit
The main principle of basis pursuit (BP) is similar to that of matching pursuit as it seeks to decompose a sig-
nal into individual atoms from a predened dictionary, yet with two important dierences. First, it includes
a minimization term to reduce the number and magnitudes of retrieved atoms, yielding a sparse repre-
sentation [e.g., Chen et al., 2001]. Second, instead of applying a deationary approach where atoms are

recognized and removed sequentially, basis pursuit identies all atoms simultaneously by casting both
steps into a single inversion problem [Bonar and Sacchi, 2010; Zhang and Castagna, 2011; Vera Rodriguez
et al., 2012].
The signal s(t) is represented as the convolution of a family of wavelets (t, n) and its associated coecient
series a(t, n) as

N
s(t) = [(t, n) a(t, n)]. (16)
n=1
where N is the number of atoms, t is time, and n is the dilation of atom (t, n) determining its frequency.
Using matrix notation, equation (16) can be rewritten as
a1

a
s = (1 2 ... N ) 2 + = Da + . (17)

a
N
where n denotes the convolution matrix of (t, n) with dilation index n, D is the wavelet dictionary, and
is the noise. The tf distribution of basis pursuit corresponds then to the set of weights a associated with the
set of atoms (t, n) coming from the dictionary D [Chen et al., 2001].
At this point, two options are possible, namely, a perfect or a sparse approximate decomposition [Chen
et al., 2001; Vera Rodriguez et al., 2012]. The latter is also known as basis pursuit denoising. In the rst case,
one solves the problem: minimize ||a||1 subject to Da = s, where ||.||1 indicates the 1 norm (i.e., minimum
of absolute numbers). In the second case, one minimizes
1
J= ||s Da||22 + ||a||1 . (18)
2
The rst term of cost function J represents the data mist term based on the 2 norm, that is, the
least-squares dierence between the observed and predicted data, whereas the second term of J is the
regularization term, based on the 1 norm. is the trade-o parameter controlling the relative strength
between the data mist and the number of nonzero coecients a. A variety of solvers exist to obtain a
solution to equation (18) [Figueiredo et al., 2007; van den Berg and Friedlander, 2008; Beck and Teboulle, 2009].
Unlike matching pursuit, basis pursuit is not a greedy algorithm. It starts from an initial model and itera-
tively improves the model by swapping wavelets from the predened dictionary. Basis pursuit converges
eventually to a local optimum which is not guaranteed for matching pursuit [Chen et al., 2001].
The performance of matching pursuit and basis pursuit is strongly dependent on the predened wavelet
dictionary. Combining dierent dictionaries to make bigger, more complete dictionaries can enhance basis
pursuit decompositions [Chen et al., 2001; Rubinstein et al., 2010] but comes at the expense of making matrix
D bigger and thus prolonging computation time.
The techniques presented in the next sections, corresponding to autoregressive methods and empirical
mode decomposition and variants, are fundamentally dierent from the previous ones as they drop the
comparison with a template performed by an inner product of signal and basis functions.
2.8. Autoregressive Methods

Autoregressive methods for periodicity analysis were introduced by G. U. Yule and G. Walker in the rst half
of the twentieth century and applied to environmental time series [Yule, 1927; Walker, 1931]. Linear autore-
gression, in combination with a moving average scheme, was then used to forecast the values of economic
time series such as price indices [Wold, 1938]. The fundamental idea of autoregressive models is to reduce a
time series to a small set of parameters. Autoregressive methods are called linear prediction lters because
they use a set of past values to estimate the future ones. This has important implications for the resulting fre-
quency localization as the signal is not strictly windowed like FFT-based methods. This is the reason why the
frequency localization of autoregressive methods sometimes seems to violate the uncertainty principle
[Marple, 1982]. Using the autoregressive coecients, specic parts of the signal can then be reconstructed.

Autoregressive models consider a time series sn as a weighted sum of p past values driven by a white noise
process vn following

p
sn = ak snk + vn , (19)
k=1
where n are the time steps, ak are the coecients of the autoregressive model of order p, and vn is a random
white noise process with zero mean and variance s2 .
The power spectrum of an autoregressive lter is then given by
() s 2
SA f = , (20)
| p |2
|1 k=1 ak ei2kf fs |
| |
where fs is the sampling frequency.
Except for the method to obtain the autoregressive coecients, the main parameter of autoregressive mod-
els is the order of the lter p. The signal window length has only one main restriction: it should be longer
than twice the autoregressive order [Ulrych and Ooe, 1983]. The number of coecients has a critical inu-
ence on the estimated spectrum as it sets the number of recoverable frequencies independent of signal
s(t). For the dierent applications, we use a grid-search algorithm to nd the minimum in a least-squares
sense of the spectral comparison between Fourier and autoregressive spectra to determine the optimum
autoregressive order [Tary et al., 2014].
As autoregressive models are all-pole lters, they are characterized by peaked spectra. Therefore, only
signals with such spectra should be good candidates for this method [Kay and Marple, 1981]. The poles
dening autoregressive lters correspond to the dierent modes of a resonator. In the following sections,
we will describe two dierent autoregressive methods able to capture time-varying frequency content.
2.8.1. Short-Time Autoregressive Method
As the name suggests, the short-time autoregressive method (STAR) is an application of autoregressive
models on short windows. Like the short-time Fourier transform, the short-time autoregressive method
deals with nonstationary spectra by assuming that signals are piecewise stationary in short, untapered win-
dows. The autoregressive spectrum, calculated from the autoregressive coecients determined using the
Burg method [Burg, 1972] and equation (20), is computed for short, overlapping windows and then merged
sequentially to obtain the nal picture. The time localization is then the same as the short-time Fourier trans-
form but is independent of the frequency localization, although the maximum autoregressive order p is still
limited by the window length. The frequency localization depends mainly on the signal-to-noise ratio and
the autoregressive order [Marple, 1982; Quirk and Liu, 1983], even though spectral estimates using the Burg
method employ the total data window [e.g., Muthuswamy and Thakor, 1998].
2.8.2. Continuous Autoregressive Model: Kalman Smoother
Rather than assuming stationarity in short windows, one can modify equation (19) to use time-varying
autoregressive coecients dened at each time step n,

p
sn = akn snk + vn , (21)
k=1
where akn are the nonstationary autoregressive coecients.

The nonstationary coecients are usually determined using stochastic or deterministic regression methods
[Kaipio and Karjalainen, 1997]. For the deterministic method, a set of basis functions with variable ampli-
tudes is used in order to model the time-varying frequency content [e.g., Hall et al., 1983]. In the case of
stochastic methods, a recursive algorithm is generally employed to determine the autoregressive coe-
cients at instant n from the autoregressive coecients at instant n 1 [e.g., Baziw and Weir-Jones, 2002].
In the following, we present a stochastic method, the Kalman smoother (KS) [Khan and Dutt, 2007; Tary
et al., 2014].
The autoregressive model assumes that the observation at instant n, yn , is a linear combination of a vector
of previous measurements mn weighted by the vector of autoregressive coecients xn , also called the state
vector, plus some random noise n . The measurement equation is then

sn = m n x n + n , (22)

where the prime denotes matrix transpose.

The Kalman lter uses a recursive scheme to predict the next estimate from the previous one. The state
variable at instant n, xn , is then obtained from the state variable at n 1, xn1 , by the transition equation
xn = Txn1 + n , (23)
where T is the state transition matrix and n is the state noise. In addition, the following assumptions
are made
1. The measurement and state noises are independent, white noise processes which are normally dis-
tributed with zero mean and known covariance matrices.
2. The state vector xn is a Gauss-Markov process independent on both measurement and state noises.
3. The initial state is dened by its mean 0|0 and covariance matrix P0|0 , given by
[ ]
0|0 = E x0 ,
[( )( ) ]
P0|0 = E x0 x 0|0 x0 x 0|0 , (24)
where the hat refers to estimates, E[] is the mathematical expectation, and the notation n|n1 corre-
sponds to the estimate at time step n given the estimates at n 1. Initial estimates for 0|0 , P0|0 and x 0|0
are generally obtained using a training data set [Tary et al., 2014].
Optimal estimates of the state vector are obtained by minimizing the covariance of the estimates at the
instant n + 1 given the set of observations s1n+1 = [s1 , ..., sn+1 ]. This is achieved by using the prediction and
update equations for the Kalman lter [Grewal and Andrews, 2001]. The estimates of the state vector can be
further smoothed by applying the Rauch-Tung-Striebel backward recursions [Rauch et al., 1965]. Using the
forward recursions of the Kalman lter together with the smoothing backward recursions corresponds to
the Kalman smoother. Additional details on this algorithm can be found in Tary et al. [2014].
2.9. Empirical Mode Decomposition and Variants
Empirical mode decomposition, developed by Huang et al. [1998], is a powerful signal analysis technique
to model nonstationary and nonlinear signal systems. Through the extraction of intrinsic mode functions
(IMFs), empirical mode decomposition captures the nonstationary feature of seismic data. A tf represen-
tation is obtained by combining empirical mode decomposition with the Hilbert transform to compute
instantaneous frequencies [Taner et al., 1979; Magrin-Chagnolleau and Baraniuk, 1999]. This is sometimes
called the Hilbert-Huang spectrum [Huang et al., 1998].
The instantaneous frequency [Taner et al., 1979] holds the promise of the highest possible time localization
since it produces a frequency at each time sample, but at the expense of a very limited frequency localiza-
tion, due to the Gabor uncertainty principle, equation (3). This often leads to rapid uctuations from time
sample to time sample. In addition, it can even produce negative frequencies which hold uncertain physical
meaning. Instantaneous frequency is thus a somewhat controversial variable [Barnes, 2007; Fomel, 2007; Han
and van der Baan, 2013]. Empirical mode decomposition aims at solving the predicament of instantaneous
frequency, as each IMF after decomposition is guaranteed to be a symmetric, narrow-band waveform,
which ensures that the instantaneous frequency of each IMF is smooth and positive.
The IMFs are computed recursively, starting with the most oscillatory one. The decomposition method uses
the envelopes dened by the local maxima and minima of the input seismic trace. Once the maxima of the
original signal are identied, a spline is used to interpolate all the local maxima and construct the upper
envelope. The same procedure is used for local minima to obtain the lower envelope. The mean of the upper
and lower envelopes is then subtracted from the initial signal. This interpolation process is pursued on the
remaining signal. This procedure is called sifting, and it terminates when the mean envelope is reasonably
close to zero everywhere. The resultant signal is designated as the rst IMF.
The rst IMF is then subtracted from the data, and the dierence is treated as a new signal on which the
same sifting procedure is applied to calculate the next IMF. The decomposition is stopped when the last IMF
has a small amplitude or becomes monotonic [Huang et al., 1998; Bekara and Van der Baan, 2009; Han and
van der Baan, 2013]. As the name implies, empirical mode decomposition is an empirical decomposition
in that no a priori decomposition basis is chosen such as sines and cosines for the Fourier transform or a
mother wavelet for the Wavelet transform.

Empirical mode decomposition has sev-

3 a
eral interesting properties that makes it an
Amplitude
attractive tool for signal analysis. It results in
0
complete signal decomposition; i.e., the origi-
nal signal is reconstructed by summing all IMFs.
3
No loss of information is incurred. The empiri-
150
b
Frequency (Hz) Frequency (Hz)

s5 cal mode decomposition is a quasi-orthogonal
s4 decomposition in that the cross-correlation
100
s3 coecients between the dierent IMFs are
50 s1 always close to zero. This minimizes energy
s2
leakage between the IMFs [Bekara and Van der
0
c Baan, 2009].
100 Even though empirical mode decomposition
oers several promising properties, some fea-
50
tures encumber its direct application: mode
STFT mixing and splitting, aliasing, and end-point
0
0 2 4 6 8 10 artefacts [Mandic et al., 2013]. Empirical mode
Time (s) decomposition is also relatively inecient
Figure 2. Synthetic signal constituted by straight, for at signals due to the extrema interpo-
frequency-modulated, and gliding harmonic components and lation step and for signals with frequency
one Morlet wavelet. (a) Time series. (b) Instantaneous frequen- components that are not well-separated.
cies (see text for equations corresponding to the signal labels).
(c) Short-time Fourier transform (STFT) computed using a Hann Two variants were recently introduced to
taper of 0.4 s with 90% overlap. overcome some of the negative features asso-
ciated with empirical mode decomposition,
namely, ensemble empirical mode decom-
position (EEMD) and complete ensemble empirical mode decomposition (CEEMD). Ensemble empirical
mode decomposition, briey speaking, is empirical mode decomposition combined with noise stabilization.
Using the injection of controlled zero mean, Gaussian white noise, ensemble empirical mode decomposi-
tion eectively reduces mode mixing [Wu and Huang, 2009; Tong et al., 2012; Mandic et al., 2013]. Adding
white Gaussian noise helps perturb the signal and enables the empirical mode decomposition algorithm
to visit all possible solutions in the nite neighborhood of the nal answer, and it also takes advantage of
the zero mean of the noise to cancel aliasing [Wu and Huang, 2009]. However, ensemble empirical mode
decomposition leaves two problems: rst, dierent noise realizations may end up with dierent numbers of
IMFs, which could mix up the IMFs after each decomposition in the frequency domain; second, ensemble
empirical mode decomposition does not maintain the complete decomposition feature of empirical mode
decomposition [Torres et al., 2011; Han and van der Baan, 2013]. In other words, contrary to empirical mode
decomposition, the sum of the IMFs obtained via ensemble empirical mode decomposition does not neces-
sarily reconstruct the original signal. The reconstruction error of ensemble empirical mode decomposition is
often acceptable when the injected noise is of small amplitude.
An important improvement of ensemble empirical mode decomposition is complete ensemble empiri-
cal mode decomposition, proposed by Torres et al. [2011], which obtains again the nal IMFs sequentially,
contrary to ensemble empirical mode decomposition. The rst IMF component is obtained in the same
way as ensemble empirical mode decomposition. The next IMF components can be calculated by the
following equation:
1
N
IMF(k+1) = E [r + Ek [vn ]]. (25)
N n=1 1 k
where Ek is the empirical mode decomposition operator which produces the k-th IMF component, rk is the
k-th decomposition residue, is a xed percentage of injected Gaussian white noise vn , and N is the number
of realizations. Complete ensemble empirical mode decomposition is a robust extension of empirical mode
decomposition methods. It leads to complete signal reconstructions.
The properties of empirical mode decomposition and its variants are not fully understood. In particular, its
exact time and frequency localizations are not predictable, although Han and van der Baan [2013] provide

150
100
50
CWT ST
0
0 2 4 6 8 10
SST MP
150
Frequency (Hz)
100
50
STAR BP
0
0 2 4 6 8 10
150
Frequency (Hz)
100
50
KS EEMD
0
0 2 4 6 8 10
Time (s)
Figure 3. Time-frequency representations of the synthetic signal of Figure 2 using the ST, STAR, KS, CWT, SST, EMD,
BP, and MP. The dierent techniques are computed using the following parameters: ST, parameter k=10; STAR, win-
dow of 0.3 s with 90% overlap and an autoregressive order of 21; KS, autoregressive order of 21; CWT and SST, Morlet
wavelet with 64 voices per octave; EEMD, 100 realizations and random noise injection with 10% of the signal maximum
amplitude; MP, Gabor dictionary; BP, Morlet wavelet dictionary; and regularization parameter =0.5.
some hypotheses. Given the uncertainty principle, the instantaneous frequency produces the highest possi-
ble time localization at each time sample, which necessarily comes along with poor frequency localization.
Based on the constant-Q band-pass feature of empirical mode decomposition [Flandrin et al., 2004], ensem-
ble empirical mode decomposition and complete ensemble empirical mode decomposition extract the
IMFs in each octave more accurately with the help of noise stabilization. The inherent frequency localiza-
tion of each individual IMF, given a white-noise input signal, is one octave with a time localization inversely
proportional to the center frequency of this octave. With increasing IMFs, each octave is centered around a
decreasing base frequency with a narrower linear frequency range. The obtained IMFs thus have an increas-
ing frequency localization at the expense of a decreasing time localization. In other words, the rst IMF has
the highest time localization and the lowest frequency localization, and vice-versa for the last IMF.
3. Applications
For each benchmark signal, the important parameters used with each technique are mentioned in the
caption of the gures.
3.1. Toy Examples
3.1.1. Synthetic Example
In order to test the performances of the aforementioned techniques on a noise-free signal, we build a
synthetic signal, shown in Figure 2, consisting of two harmonic components at 15 and 35 Hz (s1 and s2),

one frequency-modulated harmonic around 65 Hz

(s3), one harmonic gliding from 15 to 158 Hz (s4),
and one Morlet wavelet with a central frequency of
113 Hz (s5). The equations and associated instanta-
neous frequencies (IF, in brackets) corresponding to
these components are
s1(t) = 0.6 cos(70t) [IF 35]
s2(t) = 0.8 cos(30t) [IF 15]
s3(t) = 0.7 cos(130t+5 sin(2t)) [IF 65 + 5 cos(2t)]
( t
) t
Figure 4. Speech signal corresponding to a sample of 8100 8 100 8
a voice laughing and its time-frequency representation s4(t) = sin [IF ]
ln(100) 2
using the STFT. The STFT is computed using a Hann taper
of 0.025 s with 90% overlap.
( )
t2 5 fs
s5(t) = 3 exp cos(5t) [IF |t|]
2 2 log2 (N5 )
1 d(t)
where instantaneous frequencies are calculated as 2 dt
with (t) the argument of the signal components,
ln and log2 are the natural and binary logarithms, respectively, N5 is the number of samples of signal s5, and
fs is the sampling frequency.
This synthetic signal includes straight, frequency-modulated, sharp, and progressive changes in spectral
lines, as well as a Morlet wavelet with localized tf support. The decomposition results for all techniques are
presented in Figure 3.
The tf representations of this noise-free synthetic signal clearly shows the expected performances of the
conventional techniques: short-time Fourier transform (STFT), continuous wavelet transform (CWT), and S
transform (ST). Spectral lines for the short-time Fourier transform have a constant width, reecting the xed
tf localization of this method, while the variable tf localization of the continuous wavelet transform and
S transform involves the widening of the spectral lines toward higher frequencies together with the better
temporal localization of the high-frequency Morlet atom. Both the short-time Fourier transform and S trans-
form seem to retain the amplitude of the dierent components, unlike the continuous wavelet transform
which show a decrease in maximum amplitude toward higher frequencies due to the frequency widening
of the wavelet support and the wavelet energy conservation required by this technique.
The autoregressive techniques (short-time autoregressive method STAR and Kalman Smoother KS), basis
pursuit (BP), synchrosqueezing transform (SST), and ensemble empirical mode decomposition (EEMD) show
higher tf localizations close to the theoretical instantaneous frequencies, with the exception of short-time
autoregressive method and ensemble empirical mode decomposition which exhibits somewhat noisy
frequency components for those varying over time. Ensemble empirical mode decomposition maps the
continuous spectrum of the Morlet atom to a single point, contrary to the more conventional methods. Like
autoregressive methods, empirical mode decomposition is challenged when handling continuous spec-
tra [e.g., Kay and Marple, 1981]. Even though the localization of this atom is fully recovered by the Kalman
smoother, its time position is slightly shifted due to some smoothing. The tf smoothing also leads to the
introduction of spurious components such as the one at 80 Hz and 3.9 s. The separation of the Morlet
atom into two spectral peaks by the short-time autoregressive method illustrates another pitfall of autore-
gressive methods, which is spectral line splitting. Splitting occurs essentially for high autoregressive orders,
compared to the number of frequency components in the signal, and noise-free signals [Chen and Stegen,
1974; Kay and Marple, 1979]. In this case, the presence of small amplitude noise would be benecial, as it
prevents allocation of multiple autoregressive coecients to higher amplitude spectral peaks. Additional

4000
CWT ST
3000
2000
1000
0
0 0.04 0.08 0.12 0.16
SST MP
4000
STAR BP
Frequency (Hz)
3000
2000
1000
0
0 0.04 0.08 0.12 0.16
4000
KS EEMD
Frequency (Hz)
3000
2000
1000
0
0 0.04 0.08 0.12 0.16
Time (s)
Figure 5. Time-frequency representations of the speech signal shown in Figure 4 using the ST, STAR, KS, CWT, SST, EMD,
BP, and MP. The computing parameters are as follows: ST, parameter k=32; STAR, window of 0.025 s with 90% overlap and
an autoregressive order of 14; KS, autoregressive order of 14; CWT and SST, Morlet wavelet with 64 voices per octave;
EEMD, 100 realizations and random noise injection with 15% of the signal maximum amplitude; MP: Gabor dictionary;
BP, Morlet wavelet dictionary and regularization parameter =0.5.
autoregressive coecients are sometimes necessary to handle the appearance of frequency components
over time.
The sparse representation imposed by the matching pursuit (MP) algorithm and the choice of the Gabor
dictionary lead to discontinuous spectral lines and the localized mapping of the Morlet atom. We select the
Morlet wavelet dictionary for analyzing the benchmark signals with basis pursuit which is very similar to the
Gabor atoms dictionary used with matching pursuit [Addison, 2002]. The better performance of basis pursuit
compared with matching pursuit is mainly due to the optimization scheme for nding the optimal solution
in equation (18) ensuring sparsity of the solution [Chen et al., 2001; Vera Rodriguez et al., 2012]. Noticeably,
the tf representations of autoregressive methods, basis pursuit, synchrosqueezing transform, and
ensemble empirical mode decomposition are not showing the dierences in amplitude between the various
signal components.
3.1.2. Seminal Example: Speech Signal
Spectral analysis is commonly used to characterize speech signal and for pattern recognition [Rabiner and
Juang, 1993]. The dierent methods are applied to a short sample of a speech signal, specically a laughing
voice (Figure 4). The signal, sampled at 8000 Hz, is mainly composed of three harmonics at 600, 1050, and
1570 Hz (Figure 5). Its signal-to-noise ratio (SNR), based on the maximum amplitude of the noise and the
signal, is greater than 3.

Norm. amplitude
Explosion
0.5
0.5
1
STFT
Frequency (Hz)
40 Volcano-tectonic event
Gliding tremor
30
20
10
0
0 100 200 300 400 500 600
Time (s)
Figure 6. Volcanic pre-explosion recordings from the Redoubt volcano, Alaska [Hotovec et al., 2013] (top), and
time-frequency representation obtained using the STFT (bottom). The STFT is computed using a Hann taper
of 5.12 s with 90% overlap. The inset is an enlargement of two volcano-tectonic events showing their main
frequency components.
All techniques are able to pick up the three main spectral components. An additional low-amplitude spectral
line is discernible around 2050 Hz only for some of the techniques, namely, the short-time Fourier trans-
form, S transform, basis pursuit, short-time autoregressive method, and Kalman smoother. The decrease
in frequency localization towards higher frequencies prevents the continuous wavelet transform and syn-
chrosqueezing transform from detecting this low-amplitude spectral component. Even though the principle
of the S transform is very similar to the continuous wavelet transform, the frequency localization has been
adjusted using the parameter k (see equation (9)) to depict the low-amplitude line. Time-frequency rep-
resentations of nonconventional methods (i.e., autoregressive methods, synchrosqueezing transform,
empirical mode decomposition, matching pursuit, and basis pursuit) show less noise than the conven-
tional methods, which are prone to smearing and leakage. The highest frequency components are however
oscillating between 1000 and 1600 Hz for ensemble empirical mode decomposition and hence cannot be
separated. For basis pursuit, the value of the regularization parameter is the same as for the synthetic sig-
nal (0.5) which gives approximately the same weight to data mist and nonzero coecients a in the cost
function (equation (18)). This choice seems well adapted for the representation of spectral lines.
One can wonder if the rapid oscillations present in the synchrosqueezing transform and Kalman smoother
decompositions correspond to the correct picture of the signal tf content, as no other techniques show
these oscillations. Comparing in detail the synchrosqueezing transform and Kalman smoother represen-
tations, we can see that most of these oscillations are similar in both pictures, especially for the spectral
line around 1570 Hz. The fact that both techniques based on entirely dierent methodologies produce
very close results likely means that these oscillations are indeed correct. Several reasons can explain why
the other techniques are missing these rapid frequency modulations, smearing for the short-time Fourier
transform, continuous wavelet transform and S transform, averaging over a relatively long window for the
short-time autoregressive method, irregular uctuations of instantaneous frequencies for ensemble empir-
ical mode decomposition, sparsity, and size of the operator (i.e., atom family) for matching pursuit, and
amount of smoothing in the inversion for basis pursuit. A high localization in both time and frequency is
here the key to their detection.

50
CWT ST
40
30
20
10
0
0 200 400 600
SST MP
600
50
STAR BP
Frequency (Hz)
40
30
20
10
0
0 200 400 600
50
KS CEEMD
Frequency (Hz)
40
30
20
10
0
0 200 400 600
Time (s)
Figure 7. Time-frequency representations of the signal recorded at Redoubt volcano, Alaska (See Figure 6), using the ST,
STAR, KS, CWT, SST, EMD, BP, and MP. It includes numerous volcano-tectonic earthquakes as well as a gliding tremor. The
computing parameters are as follows: ST, parameter k=32; STAR, window of 6 s with 90% overlap and an autoregressive
order of 97; KS, autoregressive order of 97; CWT and SST, Morlet wavelet with 32 voices per octave; CEEMD, 100 realiza-
tions and random noise injection with 10% of the signal maximum amplitude; MP, Gabor dictionary; BP, Morlet wavelet
dictionary and regularization parameter = 3. Insets are enlargements of the same two volcano-tectonic events as those
enlarged in Figure 6 for the STFT and indicated by the dashed rectangle in the CWT plot.
3.2. Geophysical Data Examples

3.2.1. Volcanic Tremor
Volcanic tremors are observed on many volcanoes worldwide [e.g., Chouet, 1996]. They distinguish them-
selves from their relatives, the long-period events, by the absence of strong onsets with broad spectra. They
are characterized by long durations and harmonics mainly between 0.5 and 10 Hz. Potential causes range
from the resonance of uid-lled cracks [Aki et al., 1977; Frehner and Schmalholz, 2010], nonlaminar ow in
sinuous conduits [Julian, 1994; Benson et al., 2008], or repetitive events with quasi-perfect periodicity [Lees et
al., 2004].
The tremor sample presented here was recorded in 2009 at Redoubt volcano, Alaska (event 12 in Hotovec
et al. [2013]), by a short-period station (natural frequency 12 Hz, sampling frequency 100 Hz). This tremor
was preceded by swarms of volcano-tectonic events with decreasing interevent periodicity and followed by
an explosive episode in the volcano. The signal is constituted by one fundamental harmonic plus one over-
tone gliding toward higher frequencies over time, as well as numerous volcano-tectonic events and some
broadband noise up to 15 Hz (Figure 6). Hotovec et al. [2013] suggest that the tremor cause is the repeti-
tion of periodic volcano-tectonic events with decreasing interevent periodicity owing to their large number
before the tremor. Hotovec et al. [2013] note also that increasing the size of the signal window improves

the readability of harmonics and over-

tones in short-time Fourier transform
tf decompositions, coming at the
expense of time localization. This limita-
tion of the short-time Fourier transform
could be alleviated by some of the tech-
niques presented in this study, such as
autoregressive methods, empirical mode
decomposition, basis pursuit, and the
synchrosqueezing transform.
With this signal, the dierent techniques
are challenged in terms of noise sensitiv-
ity and adaptability to dierent types of tf
signatures. Volcano-tectonic events have
SNR ranging between 1 to 5, while the
tremor is approximately at the noise level.
Figure 8. (top) Raw recordings and (bottom) time-frequency repre- The decomposition results are presented
sentation obtained with the STFT of a microseismic event recorded in Figures 6 and 7. Overall, conventional
by a borehole geophone (vertical component) during a hydraulic methods are performing similarly on this
fracturing treatment. Two seismic phases are indicated, a P-to-S
converted wave, likely converted very close to the receiver, and the signal. The frequency content is mainly
S wave. The STFT is computed using a Hann taper of 0.03 s with below 15 Hz, with the exceptions of the
90% overlap. gliding tremor (from 365 s and 6.2 Hz to
507 s and 27 Hz) and its overtone (from
364 s and 11.7 Hz to 470 s and 28.6 Hz).
Volcano-tectonic events correspond to the broad vertical lines with a frequency content mainly below
15 Hz. Some of them seem to consist of several frequency components (7), especially those after the
appearance of the gliding tremor. This feature is better seen on the tf representations of the S transform,
continuous wavelet transform, synchrosqueezing transform, basis pursuit and short-time autoregressive
method (Figure 7).
On the other hand, no signicant improvements are achieved by the nonconventional techniques for the
delineation of the gliding tremors compared with the conventional ones. The tf picture for the short-time
Fourier transform presented in Figure 6 has a higher frequency localization than the one presented in
Hotovec et al. [2013], even though the same window size was used. The main parameters to improve the
short-time Fourier transform resolution and readability are the window size, the type of taper used, and the
number of points of the Fourier transform. The good performances of the continuous wavelet transform on
this signal are transferred to the synchrosqueezing transform, which shows the highest tf localization of
all techniques.
The high-amplitude noise degrades the performance of complete ensemble empirical mode decomposition
(CEEMD) technique such that the tremor and the volcano-tectonic events are not, or barely, identiable
before and after around 400 s, respectively. The noise also degrades the representation of this signal by
matching pursuit due to the absence of a separation mechanism between signal and noise components in
this algorithm. The iterative procedure and the compromise entrusted to the regularization parameter ,
which is equal to 3 for this signal, enable basis pursuit to reduce contributions from the noise.
The spectral localization as well as the performance of both autoregressive methods depends on the
SNR [Marple, 1982; Quirk and Liu, 1983]; this is particularly clear on the tf representation of the Kalman
smoother. For the present signal, a high autoregressive order is needed to account for the broadband char-
acter of the volcano-tectonic events and the background noise. The autoregressive methods are tting
peaked spectra to these broadband frequency components, leading to the horizontally hatched aspect of
this Kalman smoother representation.
Both the short-time Fourier transform and continuous wavelet transform represent well broad contin-
uous spectra, whereas autoregressive methods, the synchrosqueezing transform, and empirical mode
decomposition try to condense such spectra to a few single frequencies [Kay and Marple, 1981; Auger et al.,

CWT 450 ST
350
250
150
0 0.1 0.2 0.3
SST MP
450 STAR BP
Frequency (Hz)
350
250
150
0 0.1 0.2 0.3
450 KS CEEMD
Frequency (Hz)
350
250
150
0 0.1 0.2 0.3
Time (s)
Figure 9. Time-frequency representations of the microseismic event in Figure 8), using the ST, STAR, KS, CWT, SST, EMD,
BP, and MP. The computing parameters are as follows: ST, parameter k = 4; STAR, window of 0.04 s with 90% overlap and
an autoregressive order of 11; KS, autoregressive order of 11; CWT and SST, Morlet wavelet with 64 voices per octave;
CEEMD, 100 realizations and random noise injection with 15% of the signal maximum amplitude; MP, Gabor dictionary;
BP, Morlet wavelet dictionary and regularization parameter = 0.05.
2013; Mandic et al., 2013]. Hence, they display the gliding tremor fairly well but not necessarily the
volcano-tectonic events.
3.2.2. Microseismic Event
Microseismic events are small brittle failure events occurring when uids are injected inside a rock volume,
be it for enhanced oil recovery, geothermal operations, or carbon capture and storage for example. These
events are usually recorded by geophones deployed in boreholes to decrease surface noise contamination
and be closer to the point of uid injection into the reservoir. Their magnitude are usually negative, ranging
from 3 to 1 [Cipolla et al., 2011]. The tf representation of these events is complicated mainly by their
very short durations generally between 0.1 and 1 s. A clear separation of the dierent seismic phases is then
dicult to obtain due to the limits in time and frequency localizations of conventional tf methods [Herrera
et al., 2015].
In this section, we analyze a high-quality microseismic event recorded during a 24-stage hydraulic fracturing
experiment that took place in Rolla, British Columbia, Canada, in 2011 [Eaton et al., 2013]. The fracturing is
monitored by a toolstring of six downhole short-period geophones, with a sampling frequency of 2000 Hz.
The signal shown in Figure 8 corresponds to the recordings of a magnitude 1.7 microseismic event by the
vertical component of the deepest geophone. Its SNR is superior to 1.5 allowing an easy identication of
two wave arrivals. The two wavepackets correspond to a P-to-S converted wave, likely converted close to

the receiver location, and the S wave.

On the short-time Fourier transform
picture (Figure 8), the converted wave
has a frequency content around 320 Hz
while the S-wave consists of two main
components at 320 and 210 Hz.
The tf representations of this micro-
seismic event are presented in Figures
8 and 9. Due to the short signal win-
dow imposed by the short duration of
these kind of events, the picture of con-
ventional methods suers heavily from
spectral smearing and the dierent
Figure 10. (top) Raw recordings and (bottom) time-frequency represen- components are spread out over a fre-
tation obtained with the STFT of the Mw 9 Tohoku earthquake by the quency interval of approximately 50 Hz.
broad-band station KDAK (vertical component) of the IRIS IDA network.
Apart from matching pursuit and basis
The main seismic phases are also indicated. The STFT is computed using
a Hann taper of 100 s with 90% overlap. pursuit, the other nonconventional
methods are converting the dierent
zones into several localized lines. The
synchrosqueezing transform, basis pursuit, Kalman smoother, and complete ensemble empirical mode
decomposition representations oer the best performance in terms of time and frequency localizations, the
short-time autoregressive method representation being noisier.
On the other hand, Kalman smoother and complete ensemble empirical mode decomposition represen-
tations show rapid amplitude modulations, together with frequency modulations for complete ensemble
empirical mode decomposition, which are not corroborated by the other representations. Matching pur-
suit clearly localizes the main three frequency components but show higher smearing than the other
nonconventional methods. The basis pursuit algorithm detects approximately the same components but
focus them around their central frequencies. The regularization parameter used for this signal is equal
to 0.05 which is favoring the reduction of data mist over the number of nonzero coecients a. Matching
pursuit, basis pursuit, and complete ensemble empirical mode decomposition are able to represent the
monofrequency tails in the event coda between 0.15 and 0.3 s.
3.2.3. Megathrust Earthquake: Tohoku 2011, Mw 9
The Mw 9 Tohoku earthquake occurred on 11 March 2011, oshore the northwest coast of Japan where the
Pacic plate subducts under the Okhotsk plate [Tajima and Kennett, 2012]. The seismogram presented in
Figure 10 was recorded by the borehole station KDAK from the IRIS IDA network located in Kodiak Island
on the Aleutian trench, South Alaska. This station has a three-component broadband seismometer with
sampling frequency at 20 Hz giving an usable frequency band between around 0.003 and 10 Hz. The main
signal from the earthquake lasts for about half an hour, with a high SNR superior to 10.
The tf representations of this seismogram using the presented techniques are given in Figures 10 and
11. P waves arrive at 265 s and have a broadband frequency content up to the Nyquist frequency (not
shown in gures). S waves, arriving at 695 s, have a more localized frequency content around 0.015 Hz. The
late high-amplitude arrival corresponds to surface waves. Their dispersion curve is clearly visible on the tf
representations, starting from around 0.02 Hz at 1115 s to reach 0.05 Hz at later times.
The broadband character of the P waves gives to the short-time Fourier transform and S transform repre-
sentations a blurry aspect. Together with the continuous wavelet transform representation, they show some
smearing around the dispersion curve of the surface waves. While the S wave arrival is discernible on the
short-time Fourier transform and continuous wavelet transform representations, this is not the case for the S
transform. This arrival is actually also visible only on complete ensemble empirical mode decomposition and
perhaps the short-time autoregressive method, Kalman smoother, and basis pursuit tf plots (Figure 11).
Except for the short-time autoregressive method, nonconventional techniques elude the P waves arrival
with its broad spectrum in their tf representations.
The surface waves are well characterized by the nonconventional techniques, especially the synchrosqueez-
ing transform and Kalman smoother. The matching pursuit technique shows the dispersion curve relatively

CWT ST
0.4
0.3
0.2
0.1
0
0 1000 2000 3000
SST MP
3000
STAR BP
Frequency (Hz)
0.4
0.3
0.2
0.1
0
0 1000 2000 3000
KS CEEMD
Frequency (Hz)
0.4
0.3
0.2
0.1
0
0 1000 2000 3000
Time (s)
Figure 11. Time-frequency representations of the Mw 9 Tohoku earthquake recorded by station KDAK (vertical compo-
nent) of the IRIS network (See Figure 10), using the ST, STAR, KS, CWT, SST, EMD, BP, and MP. The computing parameters
are as follows: ST, parameter k = 4; STAR, window of 150 s with 90% overlap and an autoregressive order of 6; KS,
autoregressive order of 6; CWT and SST, Bump wavelet with 64 voices per octave; CEEMD, 100 realizations and random
noise injection with 10% of the signal maximum amplitude; MP, Ricker dictionary; BP, Morlet wavelet dictionary and
regularization parameter = 1.
well, even though it seems that this method has more diculties to delineate time-varying spectral lines, as
previously seen also on the synthetic signal (Figure 3). The high amplitude of the surface waves allows for
an easy recovery of this component by autoregressive methods, even with a low number of autoregressive
coecients (six autoregressive coecients in Figure 11).
4. Discussion
Table 2 gives the key parameters, a summary of the main tf features, and the main references for each
tf technique.
4.1. Parametrization and Adaptability

Depending on the tf template, or an a priori model, the resulting tf representations will be very contrast-
ing as illustrated by the presented applications. Obtained tf estimates are then not absolute measures
but actually reect the contents of their template. All techniques are however able to decompose and fully
reconstruct the signal under investigation, given some basis, with the exception of autoregressive methods
which are appropriate only for peaked spectra. The fact that these transforms are fully reversible indicates
that they all lead to correct tf representations. However, one method might provide more insight into
a particular signal than another one. Results from dierent techniques should be compared in order to

Table 2. Summary of Key tf Parameters and Features for the tf Techniques Presented in This Reviewa
Techniques Key Parameters Key tf Features Main References

STFT Window length and overlap, taper type, Fixed tf resolution, smearing of t f components Cohen [1989]
number of frequency bins
CWT Type of mother wavelet and its central frequency, Variable tf resolution, exponential frequency Daubechies [1992] and
number of voices and octaves axis, smearing of tf components Kumar and Foufoula-Georgiou [1997]
ST Localization parameter k, number of frequency Variable tf resolution, linear frequency Stockwell et al. [1996] and
bins axis, smearing of tf components Reine et al. [2009]
SST Same as CWT plus computation of instantaneous Reassignment of CWT tf distribution on Daubechies et al. [2011],
frequencies and reassignment the instantaneous frequencies ridges leading Thakur et al. [2013], and
to narrow-band frequency components Herrera et al. [2014]
MP Dictionary denition, algorithm stopping Sparse tf representations, deationary Mallat and Zhang [1993] and
criterion approach potentially leading to non Wang [2007]
optimum t f decompositions
BP Dictionary denition, regularization parameter Sparse tf representations, inversion Chen et al. [2001] and
scheme leading to locally optimum tf decompositions Vera Rodriguez et al. [2012]
STAR Autoregressive order, window length and Peaked spectra, no windowing and hence Makhoul [1975],
overlap, autoregressive method smearing reduction Kay and Marple [1981], and
Leonard and Kennett [1999]
KS Autoregressive order, parameters for the Smoother peaked spectra, time-varying Khan and Dutt [2007] and
training data set (noise amplitude, number autoregressive coecients Tary et al. [2014]
of realizations)
EMD Stopping criteria for IMFs extraction, Computation of instantaneous frequencies Huang et al. [1998],
extrema interpolation algorithm from IMFs, narrow-band and uctuating Flandrin et al. [2004], and
frequency components, problems of mode Bekara and Van der Baan [2009]
mixing and splitting, aliasing, end-points
artefacts
EEMD Same as EMD plus amplitude of Gaussian Same as EMD but reduces mode mixing Wu and Huang [2009],
white noise injected and number of realizations and aliasing, possibility of IMFs mixing Mandic et al. [2013], and
in the frequency domain and incomplete Tong et al. [2012]
decomposition of the signal
CEEMD Same as EEMD Same as EEMD without IMFs mixing in Torres et al. [2011] and
the frequency domain, complete decomposition Han and van der Baan [2013]
of the signal
a IMF, intrinsic mode function.
gain condence on their tf decomposition, particularly for techniques that are strongly model dependent
(e.g., matching and basis pursuits and autoregressive methods). While the various tf techniques are more
or less suitable for dierent kind of signals, these techniques can also serve two main rationales which are
time series representation and time series manipulation [Auger et al., 2013].
In addition, the presented techniques not only rely on dierent assumptions such as periodicity or peaked
spectra but also on the determination on various parameters that are critical for the tf picture. The
necessary parameters depend on the technique considered. For the short-time Fourier transform, four main
parameters inuence the tf representation: window size, taper type, overlap, and number of frequency
bins in the fast Fourier transform [Harris, 1978; Kay and Marple, 1981]. Wavelet- and dictionary-based tech-
niques both depend on the choice of the basis atom (Morlet, Ricker, bump wavelets, etc.), even though
techniques exist to estimate this basis atom from the data themselves [Mesa, 2005], or to merge several dic-
tionaries into a super-dictionary [Chen et al., 2001; Rubinstein et al., 2010]. Additional parameters aect
the tf decomposition such as the number of octaves and voices for the continuous wavelet transform
[Daubechies, 1992] and the regularization parameter for basis pursuit [Vera Rodriguez et al., 2012]. As the syn-
chrosqueezing transform depends on the continuous wavelet transform results, it is governed by the same
parameters plus those for instantaneous frequencies computation [Thakur et al., 2013].
The main parameters for empirical mode decomposition are the stopping criteria for the sifting process and
the extrema interpolation scheme (splines in this study) [Han and van der Baan, 2013; Mandic et al., 2013].
For the dierent examples presented in the application section, ensemble empirical mode decomposition
and complete ensemble empirical mode decomposition employed with the same parameters show similar

tf representations. On the other hand, the ensemble empirical mode decomposition implementation tends
to extract smoother IMFs with increasing IMF numbers due to the decrease in Gaussian noise uctuations,
which is not the case in the complete ensemble empirical mode decomposition implementation. From our
tests, complete ensemble empirical mode decomposition is more computationally intensive than ensemble
empirical mode decomposition which contradicts the conclusions of Torres et al. [2011]. We show ensemble
empirical mode decomposition results when no signicant improvement is achieved by using complete
ensemble empirical mode decomposition.
Autoregressive methods mainly depend mainly on the autoregressive order. If the autoregressive order is
too low, then the spectrum will be very smooth and some spectral peaks will be lost or misplaced, whereas
if the autoregressive order is too high, then spurious peaks can appear and spectral line splitting can occur
[Kay and Marple, 1981]. Various methods have been proposed to determine the autoregressive order, such
as spectral comparison [Tary et al., 2014], information theory-based criteria (e.g., Akaike information crite-
rion, nal prediction error) [Priestley, 1994] which weigh the mist of the model with the data against the
model complexity and higher order statistics [Lesage, 2008]. In the case of signal with noise, information
theory criteria usually do not converge toward an optimum value but decrease with increasing autore-
gressive orders. However, the curve of Akaike information criterion against autoregressive order can be
useful to determine when no signicant improvements are made by increasing the autoregressive order
and hence nd the optimum order. The other methods provide an autoregressive order that seem to give
reasonable results.
The various techniques show dierent levels of adaptability to all possible time series. Techniques that are
based on a xed basis (e.g., mother wavelet and dictionaries) or present strong model assumptions are
naturally less adjustable to dierent kinds of signals than techniques completely data driven such as empir-
ical mode decomposition. On the downside, empirical mode decomposition algorithm seems very sensitive
to the SNR as no criterion to separate noise and signal is dened in this kind of algorithm. At the same time,
empirical mode decomposition coupled with the Hilbert-Huang spectrum, which involves the computa-
tion of instantaneous frequencies, seems to limit its tf representation potential to narrow-band signals,
as illustrated by the broadband signals of the Tohoku mega-earthquake (Figure 11) and the Redoubt
volcano (Figure 7).
Algorithms of conventional methods for tf analysis (short-time Fourier transform, continuous wavelet
transform, S transform) show more adaptability to dierent spectral bandwidths owing to the exibility
of their templates. A complete representation of the tf content of a signal, including narrow-band and
broadband components, is achieved by the conventional techniques at the cost of a lower time-frequency
localization. Nonconventional techniques all rely on various assumptions to increase their localization or
separating power, at the expense of rejecting parts of the signal. The biggest challenge is then to obtain
an adaptable technique with both high time and frequency localizations. None of the presented techniques
succeeds unequivocally in this task, even though the basis pursuit and synchrosqueezing transform give
promising results.
4.2. Performances
Conventional methods, subject to the uncertainty principle, present some smearing for all examples, espe-
cially in the case of the microseismic event (Figure 9). The S transform performances are also limited by the
total number of samples of the signal N as the resulting tf array will have a size N N. Allowing a lower
number of frequencies, like for the fast Fourier transform, and decimating the signal in time can help reduce
the size of the tf array. Matching pursuit performances are found to be generally lower than those of the
other methods, including the conventional ones. The sparse representation of matching pursuit might be
inadequate for the applications shown in this study. On the other hand, matching pursuit may be more suit-
able for signal compression. Conversely, basis pursuit shows that superlocalization achieved by sparsity
constrained inversion is appropriate for most of the examples shown in this study. In addition, the com-
promise between data mist and sparsity through the resolution parameter increases the exibility of
this technique.
The reassignment method applied on instantaneous frequencies for the synchrosqueezing transform
improves the tf representation of the continuous wavelet transform in all the examples shown. How-
ever, some limitations of the continuous wavelet transform, such as the exponential frequency sampling,
that sometimes compromises its tf representations are then transferred to the tf pictures of the

synchrosqueezing transform. Empirical mode decomposition and variants ensemble empirical mode
decomposition and complete ensemble empirical mode decomposition oer high-localized tf represen-
tation for signals characterized by good SNR (e.g., the microseismic event, surface waves of the Tohoku
earthquake). On the other hand, instantaneous frequencies computed from empirical mode decomposition
IMFs have an inconsistent aspect in the presence of high-amplitude broadband noise (e.g., Redoubt volcano
example) and for close spectral lines (e.g., speech signal) [see also Bowman and Lees, 2013]. In these con-
texts, instantaneous frequencies calculated from the continuous wavelet transform picture (synchrosqueez-
ing transform) show a more continuous aspect, expected for these geophysical and audio signals, and
lower noise sensitivity. Finally, instantaneous frequencies, if coupled with decomposition into subsignals
like the synchrosqueezing transform and empirical mode decomposition, may be the paradigm of the
ultimate frequency localization, but examples like the P and S waves arrivals of Tohoku mega-earthquake
seem to indicate that some natural signals with broadband spectra are not well represented by
instantaneous frequencies.
Autoregressive methods are shown to be particularly ecient for spectral line delineation, for the speech
signal or the dispersion curve of surface waves for example. For the Kalman smoother method, a careful
attention should be given to potential spurious frequency components (e.g., synthetic signal), and poor tf
pictures are obtained for signals with low SNR (e.g., Redoubt volcano example). These methods combine
high time-frequency localization with an absence of smearing, leading to clearer tf representation over
conventional techniques. Improvements obtained by using autoregressive methods are more signicant
for signals with high SNR [Marple, 1982]. This could explain the relative absence of improvements for the
Redoubt volcano signal (Figure 7).
All considered methods allow for signal reconstruction. In this review, we focus on signal representation, not
signal manipulation. The fact that all methods permit reconstruction indicates that they all can be used for
signal manipulation. We anticipate that their suitability for manipulation could be quite dierent from the
visual impression provided by their signal representation.
4.3. Possible Improvements

An automatic, data-driven selection of the main parameters for the continuous wavelet transform (mother
wavelet, central frequency, and bandwidth) [Mesa, 2005] could be combined with the synchrosqueez-
ing transform algorithm in order to remove the usual trial-and-error procedure. This could also improve
the resulting tf representation of the continuous wavelet transform. On the other hand, the continuous
wavelet transform could be replaced by another time-frequency decomposition such as the curvelet trans-
form [Yang and Ying, 2014] or the S transform. The latter shows equal or better performances than the
continuous wavelet transform and also requires less tunable parameters. Data-driven selection of the right
dictionary could also improve the performance of methods such as matching pursuit and basis pursuit
[Hou and Shi, 2011].
The selection of the optimal autoregressive order is a longstanding issue for autoregressive methods that
is left open. Autoregressive methods relying on stochastic estimates usually involve high computational
costs that limit their real-time implementation. Modern computational tools and programming languages
should help mitigate this problem; this is also true for other methods like empirical mode decomposition
and variants.
Algorithms imposing sparsity on their tf decomposition, such as constrained inversion schemes [Vera
Rodriguez et al., 2012], might increase their tf localization in an articial way for broad-band signals. In
other words, these signals have a broad spectrum which is compacted by the decomposition method, pro-
viding in essence an incorrect representation. Hence, there is a need for algorithms that can handle both
narrowband and broadband tf representations depending on the actual underlying signal characteristics.
In a similar way as the synchrosqueezing transform is combining the continuous wavelet transform, instan-
taneous frequencies, and reassignment in a single method, mixing algorithms [e.g., Puryear et al., 2012] or
dictionaries [e.g., Rubinstein et al., 2010] might be alternatives to address these important challenges.
Finally, all considered algorithms in this review are single-trace and thus one-dimensional methods. How-
ever, multidimensional signals are very common (e.g., photos and shot gathers in reection seismology).
Likewise, multicomponent data are also abundant. In a multicomponent recording, dierent types of signals
are recorded at the same spatial position and time (e.g., three-component seismic signals, omni-directional

microphones, and red-green-blue channels in photos). Many of the presented methods have their multi-
dimensional counterparts including multidimensional wavelet transforms [Kumar and Foufoula-Georgiou,
1997], 2-D empirical mode decomposition [Nunes and Delchelle, 2009], and multidimensional autoregres-
sive models [McClellan, 1982; Hansen and Chellappa, 1988]. Both matching pursuit and basis pursuit can
easily be extended to handle multidimensional signals too [Bergeaud and Mallat, 1995; Chen et al., 2001;
Mendels et al., 2007; Bruckstein et al., 2009].
We anticipate that many of the other novel high-localization methods will also be transformed to higher
dimensions. The topic of high-localization multicomponent spectral analysis has however received relatively
less attention for geophysical applications, despite its practical importance. Generally, each component is
decomposed separately (e.g., the vertical component in a seismic section or the red channel in a photo)
and subsequently individually or jointly analyzed. Ideally however, multicomponent, multidimensional data
should be decomposed jointly, instead of individually, with a single algorithm. We nevertheless anticipate
new developments in this area such that signal characteristics present on one or more recorded compo-
nents can be successfully jointly decomposed and subsequently analyzed instead of trying to infer their
presence from the individual component decompositions. The recent introduction of the quaternion-based
Fourier transform oers much promise to achieve multicomponent, multidimensional decompositions [Ell
and Sangwine, 2007].
5. Conclusion
Time-frequency analysis is a cornerstone of geophysical signal analysis. Close examination of the theoret-
ical basis of the presented conventional and nonconventional techniques reveals their assumptions and
restrictions. These also often determine their time-frequency localization. Conventional techniques, i.e., the
short-time Fourier transform, continuous wavelet transform, and S transform, are based on a quantitative
Acknowledgments
The authors would like to thank the
comparison of the signal with a template. This is also true for some of the nonconventional techniques like
sponsors of the Microseismic Industry matching pursuit and basis pursuit, while the synchrosqueezing transform actually combines it with instan-
consortium for nancial support, and taneous frequencies determination and the reassignment method. Other nonconventional techniques
Arc Resources, Spectraseis, Nanomet-
rics, and ESG Solutions are thanked
are directly based on data modeling such as autoregressive methods and empirical mode decomposition
for their support of the eld project and variants.
of the Rolla experiment. We acknowl-
edge the IRIS consortium which The applications to a synthetic signal, a speech signal, and various geophysical time series show the
was used for access to waveform capabilities of the various techniques in terms of their adaptability to dierent frequency content,
and metadata of the KDAK and REF
time-frequency localization, and noise sensitivity. Conventional techniques perform similarly on the various
stations of the IRIS IDA and AV net-
works, respectively. All benchmark signals, clearly showing their main shortcomings which are spectral smearing and the trade-o between
signals are provided as supplemen- time and frequency localizations. The nonconventional techniques address them more or less eciently in
tary materials to this review. We
dierent ways. Out of all the presented techniques, the synchrosqueezing transform and basis pursuit are
thank the Editor-in-Chief and Matt
Hall for their thorough reviews which the only one providing high-localization and low noise sensitivity, even though the synchrosqueezing trans-
improved the quality of this review. form results are broadly dictated by those of the continuous wavelet transform. Matching pursuit and basis
The authors of the original codes
pursuit results are also bound by the suitability of the dictionary to the analyzed signal. Kalman smoother
for some techniques are also greatly
acknowledged. For the synchrosqueez- and empirical mode decomposition show high time-frequency localizations but also a high sensitivity to the
ing transform, the code uses the main signal-to-noise ratio. The actual performance of each technique is then ultimately conferred by their original
functions of the synchrosqueezing
design (e.g., matching pursuit for signal compression and autoregressive method for time series prediction).
transform toolbox (https://web.math.
princeton.edu/~ebrevdo/synsq/). For Recent developments like data-driven procedures (empirical mode decomposition), algorithmic combina-
the Kalman smoother, the main rou-
tines are adapted from those of M.E. tions (synchrosqueezing transform), or inversion schemes (basis pursuit) show potential ways to improve
Khan (http://www.cs.ubc.ca/~emtiyaz/ time-frequency representations. They might hold the keys to help circumvent longstanding hurdles such as
publications.html). Empirical mode the high-localization representation of both narrow-band and broadband spectra at the same time.
decomposition codes are modied
from those of P. Flandrin (http://perso.
ens-lyon.fr/patrick.andrin/emd.html), References
and M.E. Torres (http://bioingenieria.
Addison, P. S. (2002), The Illustrated Wavelet Transform Handbook: Introductory Theory and Applications in Science, Engineering, Medicine
edu.ar/grupos/ldnlys/metorres/re_
and Finance, IOP, Bristol, U. K.
inter.htm). The source codes for the
Aki, K., M. Fehler, and S. Das (1977), Source mechanism of volcanic tremors: Fluid-driven crack model and their application to the 1963
basis pursuit technique were kindly
Kilauea eruption, J. Volcanol. Geotherm. Res., 2(3), 259287.
provided by M. Sacchi, D. Bonar, and I.
Auger, F., and P. Flandrin (1995), Improving the readability of time-frequency and time-scale representations by the reassignment
Vera Rodriguez.
method, IEEE Trans. Signal Process., 43(5), 10681089.
The Editor on this paper was Fabio Auger, F., P. Flandrin, Y.-T. Lin, S. McLaughlin, S. Meignen, T. Oberlin, and H.-T. Wu (2013), Time-frequency reassignment and syn-
Florindo. He thanks Matt Hall and one chrosqueezing: An overview, IEEE Signal Process. Mag., 30(6), 3241, doi:10.1109/MSP.2013.2265316.
anonymous reviewer. Barnes, A. (2007), A tutorial on complex seismic trace analysis, Geophysics, 72(6), W33W43, doi:10.1190/1.2785048.

Barnhart, B. L., and W. E. Eichinger (2011), Analysis of sunspot variability using the HilbertHuang transform, Sol. Phys., 269, 439449,
doi:10.1007/s11207-010-9701-6.
Baziw, E., and I. Weir-Jones (2002), Application of kalman ltering techniques for microseismic event detection, Pure Appl. Geophys., 159,
449471.
Beck, A., and M. Teboulle (2009), A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imaging Sci., 2(1),
183202, doi:10.1137/080716542.
Bekara, M., and M. Van der Baan (2009), Random and coherent noise attenuation by empirical mode decomposition, Geophysics, 74(5),
V89V98, doi:10.1190/1.3157244.
Benson, P. M., S. Vinciguerra, P. G. Meredith, and R. P. Young (2008), Laboratory simulation of volcano seismicity, Science, 322, 249252.
Bergeaud, F., and S. Mallat (1995), Matching pursuit of images, Proc. IEEE Int. Conf. Image Process., 1, 5356, doi:10.1109/ICIP.1995.529037.
Bertino, L., G. Evensen, and H. Wackernagel (2003), Sequential data assimilation techniques in oceanography, Int. Stat. Rev., 71, 223241,
doi:10.1111/j.1751-5823.2003.tb00194.x.
Bolton, E. W., K. A. Maasch, and J. M. Lilly (1995), A wavelet analysis of Plio-Pleistocene climate indicators: A new view of periodicity
evolution, Geophys. Res. Lett., 22(20), 27532756, doi:10.1029/95GL02799.
Bonar, D., and M. Sacchi (2010), Complex spectral decomposition via inversion strategies, SEG Annu. Meeting Abstr., 14081412,
doi:10.1190/1.3513105.
Bowman, D. C., and J. M. Lees (2013), The HilbertHuang transform: A high resolution spectral method for nonlinear and nonstationary
time series, Seismol. Res. Lett., 84(6), 10741080, doi:10.1785/0220130025.
Brigham, E. O. (1988), The Fast Fourier Transform and Its Applications, Signal Processing Series, Prentice-Hall, Englewood Clis, N. J.
Bruckstein, A., D. Donoho, and M. Elad (2009), From sparse solutions of systems of equations to sparse modeling of signals and images,
SIAM Rev., 51(1), 3481.
Burg, J. P. (1972), The relationship between maximum entropy spectra and maximum likelihood spectra, Geophysics, 37, 375376.
Cazelles, B., M. Chavez, D. Berteaux, F. Mnard, J. O. Vik, S. Jenouvrier, and N. C. Stenseth (2008), Wavelet analysis of ecological time series,
Oecologia, 156(2), 287304, doi:10.1007/s00442-008-0993-2.
Chen, S. S., D. L. Donoho, and M. A. Saunders (2001), Atomic decomposition by basis pursuit, SIAM Rev., 43(1), 129159.
Chen, W. Y., and G. R. Stegen (1974), Experiments with maximum entropy power spectra of sinusoids, J. Geophys. Res., 79(20), 30193022,
doi:10.1029/JB079i020p03019.
Chouet, B. A. (1996), Long-period volcano seismicity: Its source and use in eruption forecasting, Nature, 380, 309316.
Cipolla, C., S. Maxwell, M. Mack, and R. Downie (2011), A practical guide to interpreting microseismic measurements, paper presented at
SPE-144067-MS International Conference, Soc. Pet. Eng., The Woodlands, Tex., 14-16 June.
Cohen, L. (1989), Time-frequency distributionsA review, Proc. IEEE, 77(7), 941981.
Cooley, J. W., and J. W. Tukey (1965), An algorithm for the machine calculation of complex Fourier series, Math. Comput., 19(90), 297301.
Cooley, J. W., P. A. W. Lewis, and P. D. Welch (1969), The fast fourier transform and its applications, IEEE Trans. Educ., 12(1), 2734.
Daubechies, I. (1992), Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 61, Soc. for Indust. and
Appl. Math. (SIAM), Philadelphia, Pa.
Daubechies, I., J. Lu, and H.-T. Wu (2011), Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool, Appl.
Comput. Harmonic Anal., 30(2), 243261, doi:10.1016/j.acha.2010.08.002.
Dikpati, M., J. L. Anderson, and D. Mitra (2014), Ensemble Kalman lter data assimilation in a Babcock-Leighton solar dynamo model:
An observation system simulation experiment for reconstructing meridional ow speed, Geophys. Res. Lett., 41, 53615369,
doi:10.1002/2014GL061077.
Eaton, D., M. van der Baan, J. B. Tary, B. Birkelo, N. Spriggs, S. Cutten, and K. Pike (2013), Broadband microseismic observations from a
Montney hydraulic fracture treatment, northeastern B.C., Canada, CSEG Recorder, 38(3), 4453.
Ell, T. A., and S. J. Sangwine (2007), Hypercomplex Fourier transforms of color images, IEEE Trans. Image Process., 16, 2235.
Farge, M. (1992), Wavelet transforms and their applications to turbulence, Annu. Rev. Fluid Mech., 24(1), 395457,
doi:10.1146/annurev.uid.24.1.395.
Figueiredo, M., R. Nowak, and S. Wright (2007), Gradient projection for sparse reconstruction: Application to compressed sensing and
other inverse problems, IEEE J. Sel. Top. Signal Process., 1(12), 568598.
Flandrin, P., G. Rilling, and P. Gonalvs (2004), Empirical mode decomposition as a lter bank, IEEE Signal Process. Lett., 11, 112114.
Fomel, S. (2007), Local seismic attributes, Geophysics, 72(3), A29A33, doi:10.1190/1.2437573.
Frehner, M., and S. M. Schmalholz (2010), Finite-element simulations of Stoneley guided-waves reection and scattering at the tips of
uid-lled fractures, Geophysics, 75(2), T23T36, doi:10.1190/1.3340361.
Gabor, D. (1946), Theory of communication, J. Inst. Electr. Eng., 93(26), 429457, doi:10.1049/ji-3-2.1946.0074.
Gao, R., R. Yan, and SpringerLink (2010), Wavelets: Theory and Applications for Manufacturing, Springer, Bcher, New York.
Ghil, M., et al. (2002), Advanced spectral methods for climatic time series, Rev. Geophys., 40(1), 1003, doi:10.1029/2000RG000092.
Grewal, M. S., and A. P. Andrews (2001), Kalman Filtering: Theory and Practice Using MATLAB, 2nd ed., Wiley-Intersci., New York.
Hall, M. (2006), Resolution and uncertainty in spectral decomposition, First Break, 24, 4347, doi:10.3997/1365-2397.2006027.
Hall, M. G., A. V. Oppenheim, and A. S. Willsky (1983), Time-varying parametric modeling of speech, Signal Process., 5, 267285.
Han, J., and M. van der Baan (2013), Empirical mode decomposition for seismic time-frequency analysis, Geophysics, 78(2), O9O19,
doi:10.1190/geo2012-0199.1.
Han, L., M. D. Sacchi, and L. Han (2014), Spectral decomposition and de-noising via time-frequency and space-wavenumber reassign-
ment, Geophys. Prospect., 62(2), 244257, doi:10.1111/1365-2478.12088.
Hansen, R. R., Jr., and R. Chellappa (1988), Two-dimensional robust spectrum estimation, IEEE Trans. Acoust. Speech Signal Process., 36(7),
10511066, doi:10.1109/29.1628.
Harris, F. (1978), On the use of windows for harmonic analysis with the discrete Fourier transform, Proc. IEEE, 66(1), 5183,
doi:10.1109/PROC.1978.10837.
Herrera, R. H., J. Han, and M. van der Baan (2014), Applications of the synchrosqueezing transform in seismic time-frequency analysis,
Geophysics, 79(3), V55V64, doi:10.1190/geo2013-0204.1.
Herrera, R. H., J. B. Tary, M. van der Baan, and D. W. Eaton (2015), Body wave separation in the time-frequency domain, IEEE Geosci.
Remote Sens. Lett., 12(2), 364368, doi:10.1109/LGRS.2014.2342033.
Hinich, M. J., and C. S. Clay (1968), The application of the discrete fourier transform in the estimation of power spectra, coherence, and
bispectra of geophysical data, Rev. Geophys., 6(3), 347363.
Hotovec, A. J., S. G. Prejean, J. E. Vidale, and J. Gomberg (2013), Strongly gliding harmonic tremor during the 2009 eruption of Redoubt
Volcano, J. Volcanol. Geotherm. Res., 259, 8999.

Hou, T. Y., and Z. Shi (2011), Adaptive data analysis via sparse time-frequency representation, Adv. Adaptive Data Anal., 3(1), 128.
Huang, N. E., Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu (1998), The empirical mode decom-
position and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. London A, 454, 903995,
doi:10.1098/rspa.1998.0193.
Julian, B. R. (1994), Volcanic tremor: Nonlinear excitation by uid ow, J. Geophys. Res., 99(B6), 11,85911,877.
Kaipio, J. P., and P. A. Karjalainen (1997), Estimation of event-related synchronization changes by a new TVAR method, IEEE Trans. Biomed.
Eng., 44, 649656.
Kay, S., and S. Marple (1979), Sources of and remedies for spectral line splitting in autoregressive spectrum analysis, IEEE Int. Conf. ICASSP
Acoustics Speech Signal Process., 4, 151154.
Kay, S., and S. Marple (1981), Spectrum analysisA modern perspective, Proc. IEEE, 69(11), 13801419, doi:10.1109/PROC.1981.12184.
Khan, M. E., and D. N. Dutt (2007), An expectation-maximization algorithm based kalman smoother approach for event-related
desynchronization (ERD) estimation from EEG, IEEE Trans. Biomed. Eng., 54, 11911198.
Kravchinsky, V. A., C. G. Langereis, S. D. Walker, K. G. Dlusskiy, and D. White (2013), Discovery of Holocene millennial climate cycles
in the Asian continental interior: Has the sun been governing the continental climate?, Global Planet. Change, 110, 386396,
doi:10.1016/j.gloplacha.2013.02.011.
Kumar, P., and E. Foufoula-Georgiou (1997), Wavelet analysis for geophysical applications, Rev. Geophys., 35(4), 385412,
doi:10.1029/97RG00427.
Labat, D. (2005), Recent advances in wavelet analyses: Part 1. A review of concepts, J. Hydrol., 314(14), 275288,
doi:10.1016/j.jhydrol.2005.04.003.
Lees, J. M., E. I. Gordeev, and M. Ripepe (2004), Explosions and periodic tremor at Karymsky volcano, Kamchatka, Russia, Geophys. J. Int.,
158(3), 11511167.
Leonard, M., and B. L. N. Kennett (1999), Multi-component autoregressive techniques for the analysis of seismograms, Phys. Earth Planet.
Int., 113, 247263.
Lesage, P. (2008), Automatic estimation of optimal autoregressive lters for the analysis of volcanic seismic activity, Nat. Hazards Earth
Syst. Sci., 8, 369376.
Li, C., and M. Liang (2012), A generalized synchrosqueezing transform for enhancing signal time-frequency representation, Signal Process.,
92(9), 22642274, doi:10.1016/j.sigpro.2012.02.019.
Magrin-Chagnolleau, I., and R. Baraniuk (1999), Empirical mode decomposition based time-frequency attributes, Proc. 69th SEG Meeting,
19491952.
Makhoul, J. (1975), Spectral linear prediction: Properties and applications, IEEE Trans. Acoust. Speech Signal Process., 23(3), 283296.
Mallat, S. (2008), A Wavelet Tour of Signal Processing, The Sparse Way, 3rd ed., Acad. Press, Burlington, Mass.
Mallat, S., and Z. Zhang (1993), Matching pursuits with time-frequency dictionaries, IEEE Trans. Signal Process., 41(12), 33973415,
doi:10.1109/78.258082.
Mandic, D. P., R. U. Rehman, Z. Wu, and N. E. Huang (2013), Empirical mode decomposition-based time-frequency analysis of multivariate
signals: The power of adaptive data analysis, IEEE Signal Process. Mag., 30(6), 7486, doi:10.1109/MSP.2013.2267931.
Marple, S. L. J. (1982), Frequency resolution of Fourier and maximum entropy spectral estimates, Geophysics, 47(9), 13031307.
McClellan, J. H. (1982), Multidimensional spectral estimation, Proc. IEEE, 70(9), 10291039.
Mendels, F., P. Vandergheynst, and J.-P. Thiran (2007), Matching pursuit-based shape representation and recognition using scale-space,
Int. J. Imaging Syst. Technol., 16, 162180, doi:10.1002/ima.20078.
Mesa, H. (2005), Adapted wavelets for pattern detection, in Progress in Pattern Recognition, Image Analysis and Applications, Lecture Notes
in Computer Science, vol. 3773, pp. 933944, Springer, Berlin, Heidelberg.
Miyahara, H., K. Masuda, Y. Muraki, H. Kitagawa, and T. Nakamura (2006), Variation of solar cyclicity during the Spoerer Minimum,
J. Geophys. Res., 111, A03103, doi:10.1029/2005JA011016.
Muthuswamy, J., and N. V. Thakor (1998), Spectral analysis methods for neurological signals, J. Neurosci. Methods, 83, 114.
Nunes, J.-C., and E. Delchelle (2009), Empirical mode decomposition: Applications on signal and image processing, Adv. Adapt. Data
Anal., 1(1), 125175, doi:10.1142/S1793536909000059.
Priestley, M. B. (1994), Spectral Analysis and Time Series, Acad. Press, London, U. K.
Puryear, C. I., O. N. Portniaguine, C. M. Cobos, and J. P. Castagna (2012), Constrained least-squares spectral analysis: Application to seismic
data, Geophysics, 77(5), V143V167, doi:10.1190/geo2011-0210.1.
Quirk, M., and B. Liu (1983), On the Resolution of Autoregressive Spectral Estimation, IEEE ICASSP, Boston, Mass.
Rabiner, L., and B.-H. Juang (1993), Fundamentals of Speech Recognition, Prentice Hall, Englewood Clis, N. J.
Rauch, H. E., C. T. Striebel, and F. Tung (1965), Maximum likelihood estimates of linear dynamic systems, AIAA J., 3, 14451450.
Reager, J. T., B. F. Thomas, and J. S. Famiglietti (2014), River basin ood potential inferred using GRACE gravity observations at several
months lead time, Nat. Geosci., 7, 588592, doi:10.1038/ngeo2203.
Reine, C., M. van der Baan, and R. Clark (2009), The robustness of seismic attenuation measurements using xed- and variable-window
time-frequency transforms, Geophysics, 74, WA123WA135, doi:10.1190/1.3043726.
Rioul, O., and M. Vetterli (1991), Wavelets and signal processing, IEEE Signal Process. Mag., 8(4), 1438, doi:10.1109/79.91217.
Rubinstein, R., A. M. Bruckstein, and M. Elad (2010), Dictionaries for sparse representation modeling, Proc. IEEE, 98(6), 10451057,
doi:10.1109/JPROC.2010.2040551.
Salisbury, J. I., and M. Wimbush (2002), Using modern time series analysis techniques to predict ENSO events from the SOI time series,
Nonlinear Process. Geophys., 9, 341345, doi:10.5194/npg-9-341-2002.
Segall, P., E. K. Desmarais, D. Shelly, A. Miklius, and P. Cervelli (2006), Earthquakes triggered by silent slip events on Kilauea volcano,
Hawaii, Nature, 442, 7174, doi:10.1038/nature04938.
Stockwell, R. G. (2007), Why use the S-transform?, in Pseudo-Dierential Operators: Partial Dierential Equations and Time-Frequency Analy-
sis, Fields Institute Communications, vol. 52, edited by L. Rodino, B.-W. Schulze, and M. W. Wong, pp. 279309, Am. Math. Soc., Toronto,
Ontario, Canada.
Stockwell, R. G., L. Mansinha, and R. P. Lowe (1996), Localization of the complex spectrum: The S-transform, IEEE Trans. Signal Process.,
44(4), 9981001, doi:10.1109/78.492555.
Tajima, F., and B. L. N. Kennett (2012), Interlocking of heterogeneous plate coupling and aftershock area expansion pattern for the 2011
Tohoku-Oki Mw9 earthquake, Geophys. Res. Lett., 39, L05307, doi:10.1029/2011GL050703.
Taner, M. T., F. Koehler, and R. E. Sheri (1979), Complex seismic trace analysis, Geophysics, 44(11), 10411063.
Tary, J. B., R. H. Herrera, and M. van der Baan (2014), Time-Varying Autoregressive Model for Spectral Analysis of Microseismic Experiments
and Long-Period Events, 196(1), 600611, doi:10.1093/gji/ggt400.

Thakur, G., E. Brevdo, N. S. Fuckar, and H.-T. Wu (2013), The Synchrosqueezing algorithm for time-varying spectral analysis: Robustness
properties and new paleoclimate applications, Signal Process., 93(5), 10791094, doi:10.1016/j.sigpro.2012.11.029.
Tong, W., Z. Mingcai, Y. Qihao, and Z. Huyuan (2012), Comparing the applications of EMD and EEMD on time-frequency analysis of
seismic signal, J. Appl. Geophys., 83, 2934, doi:10.1016/j.jappgeo.2012.05.002.
Torrence, C., and G. P. Compo (1998), A practical guide to wavelet analysis, Bull. Am. Meteorol. Soc., 79(1), 6178.
Torres, M., M. Colominas, G. Schlotthauer, and P. Flandrin (2011), A complete ensemble empirical mode decomposition with adaptive
noise, in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 41444147, IEEE, Prague, Czech Republic.
Ulrych, T. J., and T. N. Bishop (1975), Maximum entropy spectral analysis and autoregressive decomposition, Rev. Geophys. Space Phys.,
13(1), 183200.
Ulrych, T. J., and M. Ooe (1983), Autoregressive and Mixed Autoregressive-Moving Average Models and Spectra, Springer, Berlin, Germany.
van den Berg, E., and M. P. Friedlander (2008), Probing the Pareto frontier for basis pursuit solutions, SIAM J. Sci. Comput., 31(2), 890912,
doi:10.1137/080714488.
Vera Rodriguez, I., D. Bonar, and M. D. Sacchi (2012), Microseismic data denoising using a 3C group sparsity constrained time-frequency
transform, Geophysics, 77, V21V29, doi:10.1190/geo2011-0260.1.
Walker, G. T. (1931), On periodicity in series of related terms, Proc. R. Soc. London Ser. A, 131, 518532.
Wang, Y. (2007), Seismic time-frequency spectral decomposition by matching pursuit, Geophysics, 72(1), V13V20, doi:10.1190/1.2387109.
Wold, H. (1938), A Study in the Analysis of Stationary Time Series, Almquist and Wiksell, Stockholm, Sweden.
Wu, Z., and N. E. Huang (2009), Ensemble empirical mode decomposition: A noise-assisted data analysis method, Adv. Adapt. Data Anal.,
1(1), 141.
Yang, H., and L. Ying (2014), Synchrosqueezed curvelet transform for Two-dimensional mode decomposition, SIAM J. Math. Anal.,
46(3), 20522083.
Yule, G. U. (1927), On a method of investigating periodicities in disturbed series, with special reference to Wolfers sunspot numbers,
Philos. Trans. R. Soc. London Ser. A, 226, 267298.
Zhang, R., and J. Castagna (2011), Seismic sparse-layer reectivity inversion using basis pursuit decomposition, Geophysics, 76(6),
R147R158, doi:10.1190/geo2011-0103.1.

Tary 2014 Spectral

Uploaded by

Copyright:

Available Formats

Tary 2014 Spectral

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Tary 2014 Spectral

Uploaded by

Copyright:

Available Formats

Reviews of Geophysics

REVIEW ARTICLE Spectral estimationWhat is new? What is next?

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 1

FFT, 1024 points

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 2

Techniques Abbreviation Parametric Template

2.2. Time-Frequency Resolution

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 3

In other words, w(t) = 1 ( t

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 4

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 5

2.5. Synchrosqueezing Transform

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 6

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 7

2.8. Autoregressive Methods

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 8

where akn are the nonstationary autoregressive coecients.

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 9

where the prime denotes matrix transpose.

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 10

Empirical mode decomposition has sev-

Frequency (Hz) Frequency (Hz)

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 11

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 12

one frequency-modulated harmonic around 65 Hz

s1(t) = 0.6 cos(70t) [IF 35]

s2(t) = 0.8 cos(30t) [IF 15]

s3(t) = 0.7 cos(130t+5 sin(2t)) [IF 65 + 5 cos(2t)]

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 13

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 14

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 15

3.2. Geophysical Data Examples

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 16

the readability of harmonics and over-

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 17

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 18

the receiver location, and the S wave.

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 19

4.1. Parametrization and Adaptability

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 20

Techniques Key Parameters Key tf Features Main References

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 21

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 22

4.3. Possible Improvements

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 23

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 24

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 25

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 26

TARY ET AL. 2014. American Geophysical Union. All Rights Reserved. 27

You might also like