The Fourier decomposition

[Link] method for nonlinear and

non-stationary time series
Research
analysis
Pushpendra Singh1,2 , Shiv Dutt Joshi1 , Rakesh Kumar
Cite this article: Singh P, Joshi SD, Patney RK,
Saha K. 2017 The Fourier decomposition Patney1 and Kaushik Saha3
method for nonlinear and non-stationary time
1 Department of Electrical Engineering, Indian Institute
Downloaded from [Link] on 11 August 2023

series analysis. Proc. R. Soc. A 473: 20160871.

[Link] of Technology Delhi, Delhi, India
2 Department of ECE, Jaypee Institute of Information Technology

Received: 25 November 2016 Noida, Noida, India

Accepted: 14 February 2017 3 Samsung R&D Institute India, Delhi, India

PS, 0000-0001-5615-519X
Subject Areas:
algorithmic information theory, computer For many decades, there has been a general perception
in the literature that Fourier methods are not suitable
modelling and simulation, electrical
for the analysis of nonlinear and non-stationary data.
engineering In this paper, we propose a novel and adaptive
Fourier decomposition method (FDM), based on
Keywords: the Fourier theory, and demonstrate its efficacy for
Fourier decomposition method, Fourier the analysis of nonlinear and non-stationary time
intrinsic band functions, analytic Fourier series. The proposed FDM decomposes any data
intrinsic band functions, zero-phase filter into a small number of ‘Fourier intrinsic band
functions’ (FIBFs). The FDM presents a generalized
bank-based multivariate Fourier
Fourier expansion with variable amplitudes and
decomposition method, empirical mode variable frequencies of a time series by the Fourier
decomposition method itself. We propose an idea of zero-phase
filter bank-based multivariate FDM (MFDM), for the
Author for correspondence: analysis of multivariate nonlinear and non-stationary
time series, using the FDM. We also present an
Pushpendra Singh
algorithm to obtain cut-off frequencies for MFDM. The
e-mail: spushp@[Link]; proposed MFDM generates a finite number of band-
pushpendrasingh@[Link] limited multivariate FIBFs (MFIBFs). The MFDM
preserves some intrinsic physical properties of the
multivariate data, such as scale alignment, trend
and instantaneous frequency. The proposed methods
provide a time–frequency–energy (TFE) distribution
that reveals the intrinsic structure of a data. Numerical
computations and simulations have been carried out
and comparison is made with the empirical mode
Electronic supplementary material is available decomposition algorithms.
online at [Link]
figshare.c.3716299.

2017 The Author(s) Published by the Royal Society. All rights reserved.
 1. Introduction 2
The time–frequency representation (TFR) of a signal is a well-established powerful tool for the

[Link] Proc. R. Soc. A 473: 20160871

...................................................
analysis of time series signals. It maps a one-dimensional signal of time into a two-dimensional
signal of time and frequency, generally, by using time-varying frequency representations. There
exist three types of TFRs [1]: linear (e.g. the short-time Fourier transform (FT) and wavelet
transform), quadratic (like the Wigner distribution and the ambiguity function), nonlinear and
non-quadratic (e.g. Cohen’s non-negative distribution) TFRs. TFR is obtained by a formulation
often referred to as time–frequency distribution (TFD) and provides insight into the complex
structure of a signal consisting of several components. There are various methods for non-
stationary data processing and analysis such as the spectrogram, wavelet analysis, the Wigner–
Ville distribution, evolutionary spectrum [2] and the empirical orthogonal function expansion
Downloaded from [Link] on 11 August 2023

(EOF) (or principal component analysis or singular value decomposition). Although these
approaches have many useful applications, however, the analysis of non-stationary signals are
not well presented by these methods.
The recently proposed empirical mode decomposition (EMD) [3] provides a general method
for examining the TFD. The EMD is an adaptive signal decomposition algorithm for the analysis
of non-stationary and nonlinear signals (i.e. signals generated from nonlinear systems). The EMD
has become an established method for signal and other data analysis in various applications such
as medical studies [4–7], meteorology [3], geophysical studies [8] and image analysis [9]. The EMD
decomposes any given data into a set of finite number of narrow band intrinsic mode functions
(IMFs) which are derived directly from the data, whereas other signal decomposition techniques
such as Fourier and wavelet transforms incorporate predefined fixed basis for signal modelling
and analysis. The ensemble EMD (EEMD) is a noise-assisted data analysis method developed
in [10] to overcome the time-scale separation problem of EMD. The MEMD developed in [11] is a
generalization of the EMD for multichannel data analysis. The compact EMD (CEMD) algorithm
is proposed in [12] to reduce mode mixing, end effect and detrend uncertainty present in the
EMD. The IMFs generated by the EMD algorithm are dependent on distribution of local extrema
of signal and the type of spline used for upper and lower envelope interpolation. The traditional
EMD uses the cubic spline for upper and lower envelopes interpolation. The EMD algorithm,
proposed in [13] to reduce mode mixing and detrend uncertainty, uses non-polynomial cubic
spline interpolation to obtain upper and lower envelopes which improves orthogonality among
IMFs [14].
The property of energy preservation is important for any kind of transformation, and it
is obtained by the orthogonal decomposition of a signal in various transforms such as the
Fourier, wavelet and Fourier–Bessel representation. The energy preserving property is especially
important for the accurate and faithful analysis of three-dimensional TFE distribution of a signal.
In order to preserve energy in the signal decomposition, energy preserving EMD (EPEMD)
algorithms are proposed in [15] which ensure orthogonality among IMFs or generate a set of
IMFs which are linearly independent, non-orthogonal yet energy preserving (LINOEP) vectors.
Despite considerable success, all of the EMD algorithms are based on empirical, heuristic and
ad hoc procedures that make them hard to analyse mathematically, and EMD may suffer from
mode mixing, detrend uncertainty, aliasing and end effect artefacts [16]. There is also a lack of
mathematical understanding of the EMD algorithms, e.g. dependence of IMFs on the number
of sifting, and the stopping criteria, convergence property and stability to noise perturbation.
Despite all these limitations, the EMD is a widely used non-stationary data analysis method.
Therefore, in this paper, EMD is used as a reference to establish the validity, reliability and
calibration of the proposed methods.
There has been a general understanding in the literature (e.g. [3,16,17]) for many decades
that Fourier methods are not suitable for the analysis of nonlinear and non-stationary data, and
various reasons (e.g. linearity, periodicity or stationarity) are provided to support it. The FT is
valid under very general Dirichlet conditions (i.e. the signal is absolutely integrable with finite
number of maxima and minima, and finite number of finite discontinuities in any finite interval)
 and thus includes analysis of nonlinear and non-stationary signals as well. Therefore, in this
3
study, we explore and provide algorithms to analyse nonlinear and non-stationary data by the
Fourier method termed the Fourier decomposition method (FDM), which generates a set of a small

[Link] Proc. R. Soc. A 473: 20160871

...................................................
number of band limited Fourier intrinsic band functions (FIBFs). It is already well established that
Fourier theory-based methods are the best tool for spectrum analysis, and this study demonstrates
that the Fourier methods are also the best tool for time–frequency analysis and processing of any
signal. The power of the FT can also be realized from the fact that the analytic representation
and, hence, the Hilbert transform (HT) of a signal are inherently present in the complex Fourier
representation.
We also propose a method which captures the essential features of the MEMD, using a
zero-phase filter bank (ZPFB) approach to construct MFIBFs and the residue components. This
multivariate FDM (MFDM) generates frequency band-matched MFIBFs and residue through
Downloaded from [Link] on 11 August 2023

zero-phase filtering (ZPF). Thus, we obtain an adaptive, data-driven, ZPFB-based time–frequency

analysis method.
For the adaptive data analysis approach, the most difficult challenge has been to establish a
general adaptive decomposition method without a priori bases. In this study, we propose the FDM
and MFDM general adaptive data analysis methods that are inspired by the EMD algorithms and
their filter bank properties [18,19]. The main contributions of this study are as follows:

(i) Introduction of the FIBFs which are complete, adaptive, local, orthogonal and
uncorrelated by the virtue of construction.
(ii) Introduction of a novel FDM, completely based on the Fourier theory, to decompose given
data into a set of analytic FIBFs.
(iii) Introduction of a novel MFDM, which is based on the zero-phase filter-bank approach
that can be realized by the Fourier as well as filter theory, for multichannel data analysis.
(iv) An algorithm is also presented to obtain cut-off frequencies required for the
decomposition of data into a set of FIBFs via MFDM.

Thus, in this study a generalized Fourier expansion of a signal is obtained by the Fourier method
itself. The representation of a signal by a generalized Fourier expansion is also the main objective
of all the EMD algorithms and other data analysis methods such as synchrosqueezed wavelet
transforms (SSWTs) [20], variational mode decomposition (VMD) [21], eigenvalue decomposition
(EVD) [22], sparse TFR [23], time-varying vibration decomposition [24] and resonance-based
signal decomposition approach [25].
The rest of this paper is organized as follows: In §2, the analytic signal and EMD algorithm
are briefly presented. The proposed FDM is explained in §3. We propose the ZPFB-based MFDM
algorithm in §4. Simulation results are presented in §5. Finally, conclusions are presented in §6.

2. Brief overview of the analytic signal and empirical mode decomposition

algorithm
The EMD is a well-established adaptive signal decomposition algorithm that decomposes a
nonlinear and non-stationary data into a set of finite band-limited IMFs and residue through a
sifting process. The decomposition of a signal x(t) is expressed as the sum of  IMFs and residue

 +1

x(t) = yi (t) + r (t) = yi (t), (2.1)
i=1 i=1

where yi (t) is the ith IMF and r (t) = y+1 (t) is final residue.
All the IMFs, {yi (t)}i=1 , must satisfy two basic conditions [3]: (i) in the complete duration of
time series, the number of extrema (i.e. maxima and minima) and the number of zero crossings
are equal or differ at most by one. (ii) At any point of time in the complete duration of time
series, the average of the upper and lower envelopes, obtained by the interpolation of local
 maxima and the local minima, is zero. The first condition ensures that IMFs are narrow band
4
signals and the second condition is necessary to ensure that the IF does not fluctuate excessively
because of asymmetry of waveforms [3]. Thus, the IMFs admit amplitude–frequency modulated

[Link] Proc. R. Soc. A 473: 20160871

...................................................
(AM–FM) representation [20] (i.e. yi (t) ≈ ai (t) cos(φi (t)), with ai (t), dφi (t)/ dt = φi (t) > 0 ∀t) and
well-behaved HT [3]. For any IMF yi (t), its HT ŷi (t) is defined as convolution of yi (t) and 1/π t,
∞
i.e. ŷi (t) = (1/π )p.v. −∞ (yi (τ )/(t − τ )) dτ , where p.v. stands for the Cauchy principal value of the
integral. It is to be noted that the HT is a non-causal and unstable (i.e. not absolutely integrable)
linear time-invariant (LTI) system model. The HT is an ideal operator that in practice cannot
be physically realizable. Practically, the quadrature method [26], which under certain conditions
is equivalent to the HT, using the FT obtain analytic signal. Although the HT is global, it
emphasizes the properties of the function at the local time t. Thus, the HT is used to examine and
reveal the local properties of the function yi (t) and hence signal x(t) in (2.1). An analytic signal
Downloaded from [Link] on 11 August 2023

zi (t) can be represented by zi (t) = yi (t) + jŷi (t) = ai (t) exp ( jφi (t)), where ai (t) = [y2i (t) + ŷ2i (t)]1/2 and
φi (t) = tan−1 [ŷi (t)/yi (t)] are instantaneous amplitude (IA) and instantaneous phase (IP) of yi (t),
respectively. The instantaneous frequency (IF) of yi (t) is defined as the time derivative of IP, i.e.
ωi (t) = φi (t) = (ŷi (t)yi (t) − ŷi (t)yi (t))/(ŷ2i (t) + y2i (t)). The physical meaning of IF ωi (t) constrains that
φi (t) must be a mono-component function (i.e. an increasing or a non-decreasing function) of time.
The Bedrosian and Nuttall theorems [27,28] further impose non-overlapping spectra constraints
on the pair [ai (t), cos(φi (t))] of a signal yi (t) = ai (t) cos(φi (t)).

3. The Fourier decomposition method

We propose a class of functions termed as the FIBFs, belonging to C∞ [a, b], here with the following
formal definition.

Definition 3.1. Let x(t) be an arbitrary signal, defined in the interval [a, b], following the
Dirichlet conditions. A set of functions {yi (t) : yi (t) ∈ C∞ [a, b], 1 ≤ i ≤ M} is called a FIBF set of x(t),
if the following conditions are satisfied:

(1) x(t) = Mi=1 yi (t) + a0 , where a0 is the mean value of x(t);
b
(2) the FIBFs are zero mean functions, that is a yi (t) dt = 0, ∀i;
b
(3) the FIBFs are orthogonal functions, that is a yi (t)yl (t) dt = 0, for i = l;
(4) the FIBF admit the Analytic FIBFs (AFIBFs) representation: yi (t) + jŷi (t) = ai (t) exp( jφi (t)),
with IF ωi (t) = (d/dt)φi (t) ≥ 0, ∀t and amplitude ai (t) ≥ 0, ∀t, where ŷi (t) is obtained by the
complex exponential Fourier representation and it is equivalent to the HT of FIBF yi (t).

Thus, the AFIBFs are monocomponent signals because, physically, IF has meaning only
for monocomponent signals consisting of a single-frequency component or a narrow range of
frequencies varying as a function of time [17]. Hence, the FIBF is sum of zero mean sinusoidal
functions of consecutive frequency bands.
The main objective of this study is to develop a novel and adaptive decomposition method,
completely based on the Fourier theory, to obtain a unique representation of multicomponent
signal as a sum of the mean-value and non-stationary monocomponent signals satisfying the
properties presented in definition 3.1. The necessary conditions [3] for a set of basis vectors to
represent a nonlinear and non-stationary time series are completeness, orthogonality, locality
and adaptiveness. The FIBFs, intrinsically, follow all the necessary conditions by virtue of the
proposed decomposition.

(a) Continuous time Fourier decomposition method

In this section, we propose a continuous time FDM (CT-FDM) which decomposes a signal into a
set of FIBFs. The available data are usually of finite duration, non-stationary and generated from
systems that may be nonlinear. Let x(t) be a real-valued and time-limited signal in the interval
 
t1 ≤ t ≤ t1 + T0 . We construct a periodic extension of this signal as xT0 (t) = ∞
k=−∞ x(t − kT0 ) such 5
that x(t) = xT0 (t)w(t), where w(t) = 1, for t1 ≤ t ≤ t1 + T0 and zero otherwise. The Fourier series
expansion of xT0 (t) is given by

[Link] Proc. R. Soc. A 473: 20160871

...................................................
∞

xT0 (t) = a0 + [ak cos(kω0 t) + bk sin(kω0 t)], (3.1)
k=1
 t +T  t +T
where ω0 = 2π/T0 rad s−1 , a0 = (1/T0 ) t11 0 xT0 (t) dt, ak = (2/T0 ) t11 0 xT0 (t) cos(kω0 t) dt and
 t1 +T0
bk = (2/T0 ) t1 xT0 (t) sin(kω0 t) dt. Using the complex exponential representation of cos(kω0 t) =
[exp( jkω0 t) + exp(−jkω0 t)]/2 and sin(kω0 t) = [exp( jkω0 t) − exp(−jkω0 t)]/(2j) functions, we can
write (3.1) as
∞
1
Downloaded from [Link] on 11 August 2023

xT0 (t) = a0 + [ck exp ( jkω0 t) + c∗k exp (−jkω0 t)], (3.2)
2
k=1

where ck = (ak − jbk ) and c∗k = (ak + jbk ). From (3.2), it is clear that

xT0 (t) = a0 + Re{zT0 (t)}, (3.3)

where the analytic function

∞

zT0 (t)  ck exp ( jkω0 t) (3.4)
k=1

is complex conjugate of z̃T0 (t)  ∞ ∗
k=1 ck exp (−jkω0 t) and Re{zT0 (t)} denotes the real part of zT0 (t).
In order to obtain a set of analytic FIBFs (AFIBFs), we write zT0 (t) as


M
zT0 (t) = ai (t) exp ( jφi (t)), (3.5)
i=1

where, in forward search (low to high frequency scan (LTH-FS)) of AFIBFs, a1 (t) exp
 1  2
( jφ1 (t)) = N ck exp ( jkω0 t), a2 (t) exp ( jφ2 (t)) = N k=(N1 +1) ck exp ( jkω0 t), . . . , aM (t) exp ( jφM (t)) =
∞ k=1
c
k=(NM−1 +1) k exp ( jkω 0 t). Thus, in general, we can write


Ni
ai (t) exp ( jφi (t)) = ck exp ( jkω0 t), for i = 1, . . . , M, (3.6)
k=Ni−1 +1

with N0 = 0 and NM = ∞. The FIBFs are the real parts of AFIBFs presented in (3.6). In order to
obtain a minimum number of AFIBFs in LTH-FS, for each i, start with (Ni−1 + 1) and append
more terms until the maximum value of Ni is reached such that (Ni−1 + 1) ≤ Ni ≤ ∞ and
dφi (t)
ai (t) ≥ 0, ωi (t) = ≥ 0, ∀t, (3.7)
dt
where ai (t) and ωi (t) = 2π fi (t) are IA and IF of ith FIBF, respectively. It is easy to observe that such
a decomposition is always possible.
Similarly, in reverse search (from high to low frequency scan (HTL-FS)) of AFIBFs, we
 (N1 −1)
obtain a1 (t) exp ( jφ1 (t)) = ∞ k=N1 ck exp ( jkω0 t), a2 (t) exp ( jφ2 (t)) = k=N2 ck exp ( jkω0 t), . . . , aM (t)
(N −1)
exp ( jφM (t)) = k=1M−1 ck exp ( jkω0 t), and the lower and upper limits of sum in (3.6) would
change to k = Ni to (Ni−1 − 1), respectively, with N0 = ∞ and NM = 1. Here, we start with
(Ni−1 − 1), decrease and select minimum value of Ni such that 1 ≤ Ni ≤ (Ni−1 − 1) and (3.7) is
satisfied for i = 1, . . . , M.
We observe that (3.5) has precisely the form which in the literature [3] is termed as a
generalized Fourier expansion. Moreover, this representation is complete, orthogonal, local,
adaptive and purely Fourier based. Thus, we have presented a generalized signal specific Fourier
expansion of a time series in (3.5) by the Fourier method itself. The variable amplitude and the
IF have improved the efficiency of the expansion by expanding the signal into finite number of
analytic FIBFs in (3.5), and enabled the expansion to accommodate non-stationary data. Thus, we
 have obtained a variable amplitude and frequency representation, whereas the classical Fourier
6
expansion provides constant amplitudes and fixed-frequencies representation.
For each of the FIBFs, the amplitude ai (t) and frequency fi (t) are functions of time; therefore,

[Link] Proc. R. Soc. A 473: 20160871

...................................................
we define the three-dimensional {t, fi (t), ai (t)} TFD of amplitude as the Fourier–Hilbert spectrum
(FHS) H( f , t). The marginal Hilbert spectrum (MHS), derived from Hilbert spectrum, is defined
T
in [3]. Similarly, we define the marginal FHS (MFHS) from the FHS as H( f ) = 0 0 H( f , t) dt. The
MFHS offers a measure of total IA (or energy) contribution from each frequency in an average
sense. The frequency in either H( f , t) or H( f ) has a different meaning from the Fourier spectral
analysis [3]. The presence of energy at each frequency in MFHS H( f ) means that, in the total
duration of the signal, there is a higher likelihood for such a wave (FIBF) to have appeared
locally. The frequency in the MFHS indicates only the likelihood that an oscillation with such
a frequency exists [3]. The exact occurrence time of that oscillation is given in the full FHS.
Downloaded from [Link] on 11 August 2023

The instantaneous energy (IE) density, which can be used to measure the variation of energy
f
with time, can be defined as E(t) = 0M H2 ( f , t) df , where fM is the maximum frequency of signal.
From (3.1), we obtain the energy of signal x(t) (or power of signal xT0 (t)) by the Parseval’s theorem

as Ex = a20 + 12 ∞ [a2 + b2n ], and from (3.4) energy of the analytic signal (or power of signal zT0 (t))
∞ 2 i=1 2 n
as Ez = n=1 [an + bn ]. Therefore, relation between Ex and Ez is given by Ex = a20 + Ez /2. Hence,
the energy of zero mean signal is half of the energy of its analytic signal.
We observe that the FDM follows a nonlinear time-invariant (NTI) system model, and hence it
is a nonlinear time-invariant ZPF operations to decompose a signal into a set of FIBFs. The FDM
is nonlinear because it does not follow the principle of superposition, that is there exist n signals
 
{xl }nl=1 : FDM[ nl=1 xl ] = nl=1 FDM[xl ]. We present the following counterexample to prove that
the FDM is a nonlinear system model.

Example. Let x1 (t) = sin(10π t) and x2 (t) = sin(100π t) in time interval 0 ≤ t ≤ 1. Clearly, both
signal x1 (t) and x2 (t) are FIBFs, FDM[x1 (t)] = x1 (t) and FDM[x2 (t)] = x2 (t), thus FDM[x1 (t)] +
FDM[x2 (t)] = x1 (t) + x2 (t). However, it is very easy to verify, by the FDM algorithms,
that FDM[x1 (t) + x2 (t)] = {x1 (t), x2 (t)} which implies that FDM[x1 (t) + x2 (t)] = FDM[x1 (t)] +
FDM[x2 (t)]. It is also easy to show that FDM is a time-invariant system model as FDM[x1 (t − τ )] =
x1 (t − τ ), FDM[x2 (t − τ )] = x2 (t − τ ) and FDM[x1 (t − τ ) + x2 (t − τ )] = {x1 (t − τ ), x2 (t − τ )}.

(b) Discrete Fourier decomposition method

In practice, the continuous time signals are, generally, discretized for further processing by a
computing device; therefore, we present the FDM for discrete signal. Let, x[n], be a discrete signal
of length N. Using the discrete FT (DFT), we can write x[n] as


N−1  
j2π kn
x[n] = X[k] exp , (3.8)
N
k=0

where X[k] = (1/N) N−1n=0 x[n] exp(−j2π kn/N) is the DFT of the signal x[n]. Let N be an even
number (we can proceed in the similar fashion when N is an odd number), then X[0] and X[N/2]
are real numbers; and we can write x[n] as


N/2−1     
N−1  
j2π kn N j2π kn
x[n] = X[0] + X[k] exp +X exp( jπ n) + X[k] exp . (3.9)
N 2 N
k=1 k=N/2+1

N/2−1
Since x[n] is real, therefore, z1 [n]  k=1 X[k] exp(j2π kn/N) is complex conjugate of z2 [n] 
N−1
k=N/2+1 X[k] exp(j2π kn/N) and we can write (3.9) as
 
N
x[n] = X[0] + 2Re{z1 [n]} + X (−1)n , (3.10)
2
 Table 1. The DT-FDM algorithmic summary (LTH-FS) to obtain AFIBFs, for i = 1, . . . , M with N0 = 0, NM = (N/2 − 1) (or
7
NM = (N − 1)/2 if N is odd).
..........................................................................................................................................................................................................

[Link] Proc. R. Soc. A 473: 20160871

STEP 1. Obtain X[k] = FFT{x[n]}.

...................................................

..........................................................................................................................................................................................................
Ni j2πkn
STEP 2. Set AFIBFi = k=(Ni−1 +1) X[k] exp( N ) = ai [n] exp( jφi [n]), obtain maximum value of Ni such that
(Ni−1 + 1) ≤ Ni ≤ ( N2 − 1) and phase φi [n] of AFIBFi is a monotonically increasing function, that is, ωi [n] =

φi [n+1]−φi [n−1]
2
≥ 0, ∀n. A simple implementation of this is presented in algorithm 1.
..........................................................................................................................................................................................................

where Re{z1 [n]} denote the real part of z1 [n]. Now, we write analytic signal z1 [n] as
 
Downloaded from [Link] on 11 August 2023


N/2−1
j2π kn 
M
X[k] exp = ai [n] exp( jφi [n]), (3.11)
N
k=1 i=1
 1
where, in forward search (LTH-FS) of AFIBFs, we obtain a1 [n] exp ( jφ1 [n]) = N k=1 X[k] exp
N2  N/2−1
(j2π kn/N), a2 [n] exp ( jφ2 [n]) = k=(N1 +1) X[k] exp(j2π kn/N), . . ., aM [n] exp ( jφM [n]) = k=(NM−1 +1)
X[k] exp(j2π kn/N), and, in general,


Ni  
j2π kn
ai [n] exp( jφi [n] = X[k] exp , (3.12)
N
k=Ni−1 +1

with N0 = 0 and NM = (N/2 − 1). In order to obtain a minimum number of FIBFs in LTH-FS, for
each i, we scan from (Ni−1 + 1) to (N/2 − 1), obtain a maximum value of Ni such that (Ni−1 + 1) ≤
Ni ≤ (N/2 − 1) and phase φi [n] is a monotonically increasing function, that is an estimate of the
discrete IF
ωi [n] = (φi [n + 1] − φi [n]) ≥ 0, ∀n (3.13)
or  
φi [n + 1] − φi [n − 1]
ωi [n] = ≥ 0, ∀n (3.14)
2
and amplitude ai [n] ≥ 0, ∀n and for i = 1, . . . , M. We observe that such a decomposition always
exists. An estimation of the discrete IF by (3.14) has the following advantages: (i) it is unbiased
and has zero group delay (GD) for linear frequency modulated signals [17], and (ii) it corresponds
to the first moment in frequency of a number of TFDs [29–31].
N/2−1
Similarly, in HTL-FS for FIBFs, we obtain a1 [n] exp ( jφ1 [n]) = k=N1 X[k] exp(j2π kn/N),
 1 −1 (NM−1 −1)
a2 [n] exp ( jφ2 [n]) = N k=N2 X[k] exp(j2π kn/N), . . . , aM [n] exp ( jφM [n]) = k=1 X[k] exp( j2π kn/
N). The lower and upper limits of the sum in (3.12) would change to k = Ni to (Ni−1 − 1),
respectively, with N0 = N/2, NM = 1. In this case, for each i, we scan from (Ni−1 − 1) to 1,
obtain the minimum value of Ni such that 1 ≤ Ni ≤ (Ni−1 − 1) and phase φi [n] is a monotonically
increasing function.
It is interesting to observe that FDM provides two different views of TFE distribution, namely
low to high frequency and high to low frequency views of a signal. Depending on the signal,
both views may be the same or they may reveal two different kinds of features of the signal. The
DT-FDM is summarized in tables 1 and 2, its implementation is presented in algorithms 1 and
2. The FDM with FT and FDM with discrete time FT (DTFT) are summarized in the electronic
supplementary material.

4. Multivariate Fourier decomposition method

From (3.6) and (3.12), we observe that the operation that generates the FIBFs is actually the result
of the Fourier method based ZPF. This is another motivation, in addition to the filter bank (FB)
properties of IMFs, to use the Fourier or other methods of ZPF to decompose any data into a set
of FIBFs. The ZPF of a real valued signal x[n] by a zero-phase filter hz [n] can be obtained by two
 Table 2. The DT-FDM algorithmic summary (HTL-FS) to obtain AFIBFs, for i = 1, . . . , M with N0 = N/2 (or N0 = (N + 1)/2
8
if N is odd) and NM = 1.
..........................................................................................................................................................................................................

[Link] Proc. R. Soc. A 473: 20160871

STEP 1. Obtain X[k] = FFT{x[n]}.

...................................................
 
..........................................................................................................................................................................................................
(Ni−1 −1) j2πkn
STEP 2. Set AFIBFi = k=Ni X[k] exp N = ai [n] exp( jφi [n]), obtain minimum value of Ni such that
1≤ Ni ≤ (Ni−1 − 1) and phase φi [n] of AFIBFi is a monotonically increasing function, that is, ωi [n] =
φi [n+1]−φi [n−1]
2
≥ 0, ∀n. A simple implementation of this is presented in algorithm 2.
..........................................................................................................................................................................................................

Algorithm 1: An FDM algorithm (LTH-FS) to decompose a signal into a set of AFIBFs.

Downloaded from [Link] on 11 August 2023

input : A signal x[n] of size N

output: {AFIBFi [n], ωi [n]}M
i=1
begin
N0 = 0; i = 0; n = [0, 1, . . . , (N − 1)]; k = [0, 1, . . . , (N − 1)]; X[k] = FFT(x[n]);
if N is even then
NM = (N/2) − 1;
else
NM = (N − 1)/2;
while (k < NM ) do
i = i + 1; tmp[n] = 0;
for k = (Ni−1 + 1) to NM do 
j2πkn
tmp[n] = tmp[n] + X[k] exp N ; tmpφ[n] = unwrap_phase(tmp[n]);
(tmpφ[n+1]−tmpφ[n−1])
tmpω[n] = 2 ;
if (tmpω[n] ≥ 0, ∀n) then
AFIBFi [n] = tmp[n]; ωi [n] = tmpω[n]; Ni = k;

end

Algorithm 2: An FDM algorithm (HTL-FS) to decompose a signal into a set of AFIBFs.

input : A signal x[n] of size N
output: {AFIBFi [n], ωi [n]}M i=1
begin
i = 0; n = [0, 1, . . . , (N − 1)]; k = [0, 1, . . . , (N − 1)]; X[k] = FFT(x[n]);
if N is even then
N0 = N/2;
else
N0 = (N + 1)/2;
while (k > 1) do
i = i + 1; tmp[n] = 0;
for k = (Ni−1 − 1) to 1 do
j2πkn
tmp[n] = tmp[n] + X[k] exp( N ); tmpφ[n] = unwrap_phase(tmp[n]);
(tmpφ[n+1]−tmpφ[n−1])
tmpω[n] = 2 ;
if (tmpω[n] ≥ 0, ∀n) then
AFIBFi [n] = tmp[n]; ωi [n] = tmpω[n]; Ni = k;

end

methods: (i) the convolution method, y[n] = x[n]  hz [n] ⇒ Y[k] = X[k]Hz [k], where hz [n] = h[n] 
h[−n] which implies Hz [k] = |H[k]|2 , for a real sequence h[n], (ii) the Fourier method where we set
the frequency response of the zero-phase filter Hz [k] = 1 at desired frequency band and Hz [k] = 0
 
otherwise, and obtain output by the inverse DFT (IDFT) as y[n] = N−1 k=0 X[k]Hz [k] exp( j2π kn/N), 9

where X[k] = (1/N) N−1 k=0 x[n] exp(−j2π kn/N) is the DFT of x[n].
We use ZPF, which doesn’t shift the essential features of a signal, and propose MFDM

[Link] Proc. R. Soc. A 473: 20160871

...................................................
algorithm to generate MFIBFs and residues. We apply zero-phase high pass filtering (ZP-HPF)
with cut-off frequency fc1 to each of the components of the P-variate (P-channel) time series
{xp (t)}Pp=1 and obtain the first set of MFIBFs yp1 (t). The first set of residues is obtained as follows:
rp1 (t) = xp (t) − yp1 (t), p = 1, 2, . . . , P. (4.1)

Next, we apply ZP-HPF with cut-off frequency fc2 to first set of residues rp1 (t) and obtain
second set of MFIBFs yp2 (t). The second set of residues is obtained as

rp2 (t) = rp1 (t) − yp2 (t) p = 1, 2, . . . , P. (4.2)

Downloaded from [Link] on 11 August 2023

We can repeat this ZP-HPF procedure  times and obtain the final set of MFIBFs yp (t) and
residues (with cut-off frequency fc )
rp (t) = rp(−1) (t) − yp (t) p = 1, 2, . . . , P. (4.3)

Through the addition of (4.1)–(4.3) we obtain an expression, similar to (2.1), for a P-variate time
series as


xp (t) = ypi (t) + rp (t), p = 1, 2, . . . , P. (4.4)
i=1
When we use the Fourier-based ZPF to obtain MFIBFs, as in (3.6), the first three conditions
(definition 3.1) of FIBFs are fully satisfied and the fourth one is approximately satisfied. Obviously,
the fourth condition of FIBFs cannot be guaranteed simultaneously for all the P-channel data. This
is similar to the problem encountered in multivariate EMD (MEMD) algorithm in which, the first
condition of IMF is not imposed in the derivation of multivariate IMFs [11].
The question is, how can cut-off frequencies (CFs) fc1 , fc2 , . . . , fc corresponding to zero-phase
high pass filters hz1 (t), hz2 (t), . . . , hz (t) be obtained? There are various ways in which these cut-
off frequencies can be selected, e.g. dyadic (fc1 = fM /2, fc2 = fM /22 , . . . , fc = fM /2 , where fM is the
maximum frequency of a signal x(t) and for the sampled signal, maximum frequency is (Fs /2) half
of the sampling frequency), non-dyadic, uniform and non-uniform CFs. We can take the FT of a
signal x(t) to obtain its spectrum details and develop a strategy to decide CFs.
For a narrowband signal, we define the ratio of centre frequency ( f˜Ci ) to bandwidth (BW) as

f˜Ci ( fHi + fLi )

m= , f˜Ci = , (4.5)
( fHi − fLi ) 2
where fHi is the highest frequency and fLi is the lowest frequency of ith band of a filter bank. From
(4.5), we obtain  
(2m − 1) 1
fLi = fHi , m > . (4.6)
(2m + 1) 2
From (4.5) and (4.6), we observe that the ratios, for the consecutive ith and (i + 1)th bands, of
centre frequencies f˜Ci , cut-off frequencies fci and BWs ( fHi − fLi ) can be taken as a constant, that is

f˜Ci fci ( fHi − fLi )

= = = l. (4.7)
f˜Ci+1 fci+1 ( fHi+1 − fLi+1 )
From (4.5) to (4.7), we obtain l = (2m + 1)/(2m − 1) or m = (1/2)(l + 1)/(l − 1) with l > 1, and as
[m → ∞, l → 1], [m → 1/2, l → ∞]. We have the the liberty to select any suitable value of l or m,
and the greater the value of m (or lesser the value of l) the narrower the band. In the case of a
dyadic FB, l = 2 and m = 1.5 are fixed values. If required, we can vary the value of m (or l) for
each band rather than taking the fixed value. Thus, we propose a compact way to decide cut-off
frequencies, as summarized in table 3.
In MFDM, we can use zero-phase low pass filtering (ZP-LPF) in place of ZP-HPF to decompose
a signal in order of residue to the first MFIBFs, that is rp (y), yp (t), . . . , yp1 (t). We use ZPF as it
 Table 3. An algorithm to obtain cut-off frequencies fci .
10
..........................................................................................................................................................................................................

1. Select suitable value of m and set fH1 = Fs /2.

[Link] Proc. R. Soc. A 473: 20160871

..........................................................................................................................................................................................................

...................................................
2. Set fci = [(2m − 1)/(2m + 1)]fHi .
..........................................................................................................................................................................................................

3. Set fHi+1 = fci .

..........................................................................................................................................................................................................

4. Repeat step 2 to 3 for i = 1, 2, . . . , .

..........................................................................................................................................................................................................

Table 4. Comparisons among the Fourier, wavelet, EMD-Hilbert and proposed FDM in data analysis.

Fourier wavelet EMD-Hilbert FDM

a priori → adaptive
Downloaded from [Link] on 11 August 2023

basis a priori a priori adaptive

..........................................................................................................................................................................................................

frequency convolution: global convolution: regional differentiation: local differentiation: local

..........................................................................................................................................................................................................

uncertainty yes yes no no

..........................................................................................................................................................................................................

presentation frequency–energy frequency–time– frequency–time– frequency–time–

energy energy energy
..........................................................................................................................................................................................................

nonlinear no no yes yes

..........................................................................................................................................................................................................

non-stationary no yes yes yes

..........................................................................................................................................................................................................

harmonics yes yes no yes → no

..........................................................................................................................................................................................................

theoretical base complete complete empirical complete

..........................................................................................................................................................................................................

preserves salient features (such as maxima and minima) in the filtered time waveform exactly at
the time where these features occur in the unfiltered waveform. It is pertinent to note that the
conventional (non-zero phase) filtering shifts feature in a signal and hence cannot be employed.
The ZPF of time series can be obtained through the non-causal finite impulse response (FIR) or
infinite impulse response (IIR) filters.
As for the MEMD and noise-assisted MEMD (NA-MEMD) [19], the proposed MFDM
algorithm produces equal numbers of scale-aligned MFIBFs for all the channels, preserving joint
channel properties that make it suitable for direct multichannel modelling. The proposed FDM
does not suffer from mode mixing, detrend uncertainty and end effect artefacts as extraction of
FIBFs does not depend on distribution of local extrema across the range of signal. A comparison
of the EMD algorithm with Fourier and Wavelet is presented in [32] (the first four columns of the
table 4) using the basis function, frequency computation, uncertainty principle, presentation of
results, nonlinearity and non-stationarity of data, harmonics present in the signal representation
and theoretical base of the methods. We use the same parameters to present comparisons among
the Fourier, Wavelet, EMD and FDM for data analysis in table 4.

5. Simulation and numerical results

Matlab software of the MEMD ([Link]
[Link]), EMD and EEMD ([Link] available online, has been
used for simulation results. The objectives of the simulations are (i) to compare the results
obtained by the decomposition of signals (widely used in the literature) from the proposed
algorithms as well as from the EMD, EEMD and MEMD algorithms; (ii) to demonstrate the
efficacy of the proposed FDM for the TFE distribution analysis of simulated as well as real-
life (such as earthquake and speech) signals, and compare results with the EMD, EEMD and
wavelet transform; (iii) to show that the proposed FDM, like EMD, is not limited by the
uncertainty principle. In order to prove the last point, we use the unit sample sequence, which
is concentrated in time, and by the proposed method we show that it is also concentrated
in frequency. A Matlab code of the proposed FDM has been made publicly available at
 (a) 1 11
x (t)
0
−1

[Link] Proc. R. Soc. A 473: 20160871

...................................................
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

(b)
0.05
y1–y3

0
−0.05
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
0.02
y4–y5

0
−0.02
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
1
y6–y7
Downloaded from [Link] on 11 August 2023

0
−1
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

(c)
1
0
y1

−1
−2
−3
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
2
0
y2

−2

0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

2
y4–y5

0
−2
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
0.05
y6–y10

0
−0.05
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

(d) 0.01
0
y1

−0.01
−0.02
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
0.04
0.02
y2

0
−0.02
−0.04
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
0.2
0.1
y4–y5

0
−0.1
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
1
y6–y10

−1
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
time

Figure 1. Decompositions of a signal x(t) (a) which is a sinusoid mixed with intermittent interference through (b) proposed
FDM (c) EMD and (d) EEMD algorithm. (Online version in colour.)

([Link]
osition_Method_FDM).

(a) Intermittency and mode mixing

The intermittency in the time series is a main cause of mode mixing, as shown for signal x(t) in
figure 1a, and mode splitting in the EMD algorithm [10]. These issues are mitigated by EEMD and
NA-MEMD. The decomposition of a signal x(t) through the FDM algorithm is shown in figure 1b
 x1(t) + n1(t) x2(t) + n2(t) x3(t) + n3(t) x4(t) + n4(t)
5 5 5 5 12
0 0 0 0
−5 −5 −5 −5
0 0.5 1.0 0 0.5 1.0 0 0.5 1.0 0 0.5 1.0

[Link] Proc. R. Soc. A 473: 20160871

1 1 1 1

...................................................
0 0 0 0
−1 −10 100 200 300 −1 0 100 200 300 −10 100 200 300
0 100 200 300
1 1 1 1
0 0 0 0
−1 −1 −1 −1
0 100 200 300
1 10 100 200 300 1 0 100 200 300 10 100 200 300
0 0 0 0
−1 −1 −1 −1
0 100 200 300
1 10 100 200 300 1 0 100 200 300 10 100 200 300
0 0 0 0
−1 −1 −1 −1
10 100 200 300 10 100 200 300 1 0 100 200 300 10 100 200 300
0 0 0 0
Downloaded from [Link] on 11 August 2023

−1 −1 −1 −1
0 100 200 300
1 10 100 200 300 1 0 100 200 300 10 100 200 300
0 0 0 0
−1 −1 −1 −1
0 100 200 300 0 100 200 300 0 100 200 300 0 100 200 300

Figure 2. MFDM applied to a quadri-variate tone-noise mixture (top row) generates perfectly aligned intrinsic modes in all the
four channels (fourth to seventh row); the second and third row contain the noise components. (Online version in colour.)

without end effect artefacts. The EMD and EEMD generate 10 IMFs (y1 to y10 ), FDM yields seven
FIBFs (y1 to y7 ), where the sum of y1 to y3 , y4 to y5 , y6 to y10 is represented by y1 –y3 , y4 –y5 and
y6 –y10 , respectively. The FDM is able to localize the mono-component sinusoid within a single
FIBF, outperforming the EMD (figure 1c) and EEMD (figure 1d) in terms of reconstruction error,
end effect artefacts and orthogonalilty of generated components. On the same machine (Intel core
i3 CPU 530 @ 2.93 GHz), the computation time for FDM, EMD and EEMD are 0.56 s, 0.21 s and
77.18 s, respectively. The ensemble size for EEMD was 500 with a 16.94 dB signal-to-noise ratio.

(b) Multivariate data decomposition

We have used quadri-variate time series signal, which is a summation of sinusoids (with
combination of frequencies f1 = 4 Hz, f2 = 8 Hz, f3 = 16 Hz, f4 = 32 Hz) and Gaussian white noise

with zero mean and standard deviation of 0.2, that is xj (t) = 4i=1 sin(2π fi t) + nj (t) (for j = 1, . . . , 4).
The 32 Hz sinusoid is present in the first, second and fourth channels; the 16 Hz sinusoid is present
in the first, second and third channels; the 8 Hz sinusoid is common to all channels; the 4 Hz
sinusoid is present in the first, third and fourth channels. On the same machine (Intel core i3
CPU 530 @ 2.93 GHz), the computation time for MFDM is 0.45 s and for MEMD is 69.5 s in this
simulation. The MFDM algorithm, similar to MEMD, generates perfectly aligned intrinsic bands,
as is clear from the figure 2, in all the four channels, whereas MFDM is computationally more
efficient.

(c) Time–frequency–energy analysis of mixture of linear chirp and frequency modulated

sinusoid
Figure 3 shows TFE estimates for a non-stationary signal, which is a mixture of a linear chirp and
FM sinusoid, obtained using the FDM and EMD algorithm. There is an enhanced TFE tracking
when using FDM with LTH-FS. However, the TFE result with HTL-FS is similar to that obtained
with EMD algorithm, as can be seen from figure 3.

(d) Intrawave frequency modulation

We first decompose the following signal that has intrawave frequency modulation and it is
considered challenging because the IF itself has very high frequency modulation [33]
1 cos(32π t + 0.2 cos(64π t))
x(t) = + . (5.1)
1.2 + cos(2π t) 1.5 + sin(2π t)
 (a) FDM: time−frequency−energy estimate of FIBFs (LTH-FS)
13
3000

[Link] Proc. R. Soc. A 473: 20160871

2500

...................................................
frequency (Hz)

2000
1500
1000
500

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
time (s)
Downloaded from [Link] on 11 August 2023

(b) FDM: time−frequency−energy estimate of FIBFs (HTL-FS)

3000
frequency (Hz)

2500

2000

1500

1000

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
time (s)
(c) EMD: time−frequency−energy estimate of IMFs
2000
1800
1600
frequency (Hz)

1400
1200
1000
800
600
400
200
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
time (s)

Figure 3. The TFE analysis of non-stationary signal which is mixture of linear chirp and FM sinusoid: FDM (a,b) and EMD (c).
(Online version in colour.)

Next, we consider a model wave

x(t) = cos(ωt +  sin ωt) (5.2)

that satisfies the following highly nonlinear differential equation [3]

d2 x(t)
+ [ω + ω cos(ωt)]2 x(t) − [ω2 sin(ωt)] 1 − x2 (t) = 0 (5.3)
dt2
with ω = 1 and  = 0.5. We demonstrate through figures 4–7 that our method applies well to
these challenging cases with good accuracy. For input signal (5.1), FDM yields two FIBFs and
TFE distribution without end effects; however, EMD generates three IMFs and TFE distribution
has end effects (figure 6). The main difference between the performance of the FDM and the
 6
14
5

[Link] Proc. R. Soc. A 473: 20160871

...................................................
4

3
amplitude

0
Downloaded from [Link] on 11 August 2023

−1

−2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
time (s)

Figure 4. The signal x(t) (solid line), sum of the mean value and lowest frequency FIBF (dashed line) of (5.1). (Online version
in colour.)

(a)
4
amplitude

2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
(b)
5
amplitude

0
−5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
(c)
2
amplitude

−2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
(d) ×10−17
5
amplitude

−5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
time (s)

Figure 5. The mean value, FIBF1, FIBF2 and highest frequency component by FDM (a–d). The mean value and highest frequency
component correspond to term X[0] and X[N/2](−1)n of (3.10), respectively. Sum of the mean value, FIBF1, FIBF2 and highest
frequency component exactly synthesize signal x(t) given by (5.1). (Online version in colour.)

EMD algorithm, for signal (5.2), is the end effects (figure 7). Moreover, both methods are
able to reveal intrawave frequency modulation and nonlinearity present in the signals. These
examples clearly demonstrate that the proposed FDM can indeed be applied for the analysis of
nonlinear signals.
 (a) FDM: time−frequency−energy estimate of FIBFs (LTH-FS)
15
18
16

[Link] Proc. R. Soc. A 473: 20160871

...................................................
frequency (Hz) 14
12
10
8
6
4
2
0 0.1 0.2 0.3 0.4
0.5 0.6 0.7 0.8 0.9 1.0
time (s)
Downloaded from [Link] on 11 August 2023

(b) EMD: time−frequency−energy estimate of IMFs

30
25
frequency (Hz)

20
15
10
5

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
time (s)

Figure 6. The TFE analysis of signal x(t) given by (5.1): FDM (a) and EMD (b). Notice end effect artefacts in TFE plot by EMD.
(Online version in colour.)

(a) FDM: time−frequency−energy estimate of FIBFs (LTH-FS)

0.19
0.18
frequency (Hz)

0.17
0.16
0.15
0.14
0.13

0 5 10 15 20 25 30
time (s)

(b) EMD: time−frequency−energy estimate of IMFs

0.20
frequency (Hz)

0.15

0.10

0.05

0 5 10 15 20 25 30
time (s)

Figure 7. The TFE analysis of signal x(t) given by (5.2): FDM (a) and EMD (b). Notice the end effect artefacts in TFE plot by EMD.
(Online version in colour.)
 (a) FDM: time−frequency−energy estimate of FIBFs (LTH-FS) 16
45

[Link] Proc. R. Soc. A 473: 20160871

40

...................................................
frequency (Hz)

35
30
25
20
15
10
5
0
0 1 2 3 4 5 6 7 8 9 10
Downloaded from [Link] on 11 August 2023

time (s)

(b) FDM: time−frequency−energy estimate of FIBFs (HTL-FS)

45
40
frequency (Hz)

35
30
25
20
15
10
5
0
0 1 2 3 4 5 6 7 8 9 10
time (s)
(c) EMD: time−frequency−energy estimate of IMFs
45
40
35
frequency (Hz)

30
25
20
15
10
5
0
0 1 2 3 4 5 6 7 8 9 10
time (s)

Figure 8. The TFE analysis of the Gaussian white noise (with zero mean and unit variance): FDM (a,b) and EMD (c). The FDM
is able to decompose data into distinct frequency bands, whereas bands are not clearly separated by EMD. (Online version
in colour.)

(e) Analysis of a white Gaussian noise

Figure 8 shows the TFE analysis of a white Gaussian noise (with zero mean, unit variance,
1024 samples and sampling frequency Fs = 100 Hz) obtained from the FDM and EMD algorithm.
Clearly, both LTH-FS and HTL-FS views of TFE are similar and the entire data are decomposed
into a set of FIBFs of distinct frequency bands. Figure 9 shows the power spectral density (PSD)
plot of the same white Gaussian noise with the FDM and the EMD algorithm. The FDM has
decomposed the entire data into narrowband and orthogonal FIBFs; however, EMD is unable to
separate the bands and there is severe energy leakage among the IMFs. Both LTH-FS and HTL-FS
views of PSD look similar but FIBFs have different frequency bands. For example, approximately
 (a) FDM: power spectral density of FIBFs (LTH-FS)
50 17
power/frequency (dB/Hz) 0

[Link] Proc. R. Soc. A 473: 20160871

−50

...................................................
−100
−150
−200
−250
−300
−350
−400
0 5 10 15 20 25 30 35 40 45 50
Downloaded from [Link] on 11 August 2023

frequency (Hz)

(b) FDM: power spectral density of FIBFs (HTL-FS)

50
power/frequency (dB/Hz)

0
−50
−100
−150
−200
−250
−300
−350
−400
0 5 10 15 20 25 30 35 40 45 50
frequency (Hz)

(c) EMD: power spectral density of IMFs

50
power/frequency (dB/Hz)

0
−50
−100
−150
−200
−250
0 5 10 15 20 25 30 35 40 45 50
frequency (Hz)

Figure 9. The PSD analysis of the Gaussian white noise (with zero mean and unit variance): FDM (a,b) and EMD (c). The FDM
is able to decompose data into distinct frequency bands, whereas bands are not clearly separated by EMD. (Online version
in colour.)

40 Hz is the cut-off frequency for one of the bands in PSD (LTH-FS), whereas it is the centre
frequency of one band in the other PSD (HTL-FS) view. There are enhanced TFE and PSD tracking
when using FDM.

(f) Time–frequency–energy analysis of unit sample sequence

The unit sample sequence is defined as δ[n − n0 ] = 1 at n = n0 and zero otherwise. By using
π
the relation z[n] = (1/π ) 0 X(ω) exp( jωn) dω, we obtain the analytic representation of x[n] =
 (a)
1.0 18
Re{z[n]}
0.5

[Link] Proc. R. Soc. A 473: 20160871

...................................................
0
0 50 100 150 200 250 300 350 400
(b)
1.0
Im{z[n]}

0
−1.0
0 50 100 150 200 250 300 350 400
(c)
1.0
|z[n]| = a[n]

0.5
0
Downloaded from [Link] on 11 August 2023

0 50 100 150 200 250 300 350 400

Figure 10. The analytic representation of δ[n − n0 ] (with, n0 = 199, sampling frequency Fs = 100 Hz, length N = 400): real
part of z[n] (a), imaginary part of z[n] (b) and absolute value of z[n] (c). (Online version in colour.)

2
amplitude

0
–2
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
time (s)
0.02
amplitude

–0.02
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
time (s)
1
amplitude

0
–1
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
time (s)
×10–3
5
amplitude

–5
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
time (s)

Figure 11. The mean-value, FIBF1, FIBF2 and highest frequency component plots by the FDM of unit sample sequence
δ[n − n0 ] (with, n0 = 199, sampling frequency Fs = 100 Hz, length N = 400). (Online version in colour.)

δ[n − n0 ] ⇔ X(ω) = exp(−jωn0 ) as

sin(π (n − n0 )) + j[1 − cos(π (n − n0 ))]
z[n] = = a[n] exp( jφ[n]), (5.4)
π (n − n0 )
where the real part of z[n] is δ[n − n0 ] = sin(π (n − n0 ))/π (n − n0 ), a[n] = |sin((π/2)(n − n0 ))/
(π/2) (n − n0 ))|, φ[n] = (π/2)(n − n0 ) and hence ω[n] = π/2 which corresponds to half of the
Nyquist frequency, that is Fs /4 Hz. Figure 10 shows the plots of the real part, imaginary part
and absolute value of z[n] with n0 = 199, sampling frequency Fs = 100 Hz, length N = 400.
Theoretically, this clearly indicates that most of the energy of signal δ[n − n0 ] is concentrated at
time t = 1.99 s (n0 = 199) and frequency f = 25 Hz. Figure 11 shows the plots of FIBFs obtained by
the FDM and figure 12 shows the TFE analysis of unit sample sequence δ[n − n0 ] from the FDM,
EEMD and continuous wavelet transform (CWT) methods. These figures clearly indicate that the
 (a) FDM: time–frequency–energy estimate of FIBFs (LTH-FS)
19
30

[Link] Proc. R. Soc. A 473: 20160871

...................................................
25

frequency (Hz)
20

15

10

5
Downloaded from [Link] on 11 August 2023

0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

time (s)
EEMD: time–frequency–energy estimate of IMFs
(b)
40
35
30
frequency (Hz)

25
20
15
10
5
0
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
time (s)
(c) time–frequency–energy estimate by CWT
50
45
40
pseudo−frequency (Hz)

35
30
25
20
15
10
5
0
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
time (s)

Figure 12. The TFE analysis of unit sample sequence δ[n − n0 ] (with, n0 = 199, sampling frequency Fs = 100 Hz, length
N = 400) by the: FDM (a), EEMD (b) and continuous wavelet transform (CWT) (c). The energy is concentrated in both time
and frequency in TFE plot by the FDM (i.e. it is not limited by uncertainty principle), whereas there is spreading of the energy
in TFE plots by EEMD and CWT. (Online version in colour.)

TFE plot obtained by the FDM is the same as the theoretical estimation, whereas there is energy
spread over a range of frequencies and lack of accuracy in the TFE plots obtained by the EEMD
and CWT methods.
The IF is defined through differentiation of phase rather than integration and hence overcomes
the restriction of the uncertainty principle [34]. One can obtain arbitrary precision in time and
 (a) El Centro earthquake 18 May 1940. North−south component
0.4 20
0.3

[Link] Proc. R. Soc. A 473: 20160871

...................................................
0.2

acceleration (G)
0.1

−0.1

−0.2

−0.3
0 10 20 30 40 50 60
Downloaded from [Link] on 11 August 2023

time (s)
(b) ×10−4 power spectral density
1.2

1.0
spectral density

0.8

0.6

0.4

0.2

0
0 5 10 15 20 25
frequency (Hz)

marginal spectrum by FDM

(c) 1.4

1.2

1.0
spectral density

0.8

0.6

0.4

0.2

0
0 5 10 15 20 25
frequency (Hz)

(d) 0.7
marginal spectrum by EMD

0.6

0.5
energy density

0.4

0.3

0.2

0.1

0
0 5 10 15 20 25
frequency (Hz)

Figure 13. Plots of (a) El Centro earthquake 18 May 1940 north–south component, (b) Fourier based power spectral density
(PSD), (c) marginal spectrum by FDM and (d) marginal spectrum by EMD. (Online version in colour.)
 (a) E(t) by FDM
0.50 21
0.45

[Link] Proc. R. Soc. A 473: 20160871

...................................................
0.40
0.35
energy density 0.30
0.25
0.20
0.15
0.10
0.05
Downloaded from [Link] on 11 August 2023

0
0 10 20 30 40 50 60
time (s)

(b) E(t) by EMD

0.9
0.8
0.7
energy density

0.6
0.5
0.4
0.3
0.2
0.1
0
0 10 20 30 40 50 60
time (s)

Figure 14. The instantaneous energy density E(t) plots of El Centro earthquake data using the: (a) FDM and (b) EMD. (Online
version in colour.)

frequency resolution subject only to the sampling rate in the case of discrete time signals. The
uncertainty principle, in the context of signal analysis, states that the finer the time resolution
one wants, the cruder would be the resulting frequency resolution. However, this example
clearly demonstrates that the FHS obtained by FDM is indeed not limited by the uncertainty
principle and a signal can be highly concentrated in both time and frequency domain, as shown
in figure 12a.
Discussion: It is to be noted that the FHS provides time and average frequency distribution. In
order to explain the average frequency effects, we consider a sum of sinusoids of equal amplitudes
 N
x(t) = N k=1 A cos(ω0 kt). Its analytic representation is given by z(t) = k=1 A exp( jω0 kt) =
(A sin(ω0 (N/2)t)/sin(ω0 t/2)) exp( jω0 ((N + 1)/2)t) = a(t) exp( jφ(t)), which implies phase φ(t) =
(ω0 ((N + 1)/2)t) and hence IF f (t) = (1/2π )ω0 ((N + 1)/2). This is what we observe in these figures
(especially figure 12a). It is also to be noted that if amplitudes of the constituent sinusoids are
not equal, then the resultant IF would not be a constant (average) frequency, that is it would be
a time-varying frequency.

(g) Application to earthquake data analysis

Earthquake time series signals are nonlinear and non-stationary in character. The El Centro
earthquake data (sampled at Fs = 50 Hz) has been taken from ([Link]
[Link]) and is shown in figure 13. The critical frequency range that matters in the structural
 (a) FDM: time–frequency–energy estimate of FIBFs (HTL-FS)
25 22

[Link] Proc. R. Soc. A 473: 20160871

...................................................
20

frequency (Hz)
15

10

5
Downloaded from [Link] on 11 August 2023

0
0 5 10 15 20 25 30 35 40 45 50
time (s)
(b) EMD: time–frequency–energy estimate of IMFs

20

15
frequency (Hz)

10

0
0 5 10 15 20 25 30 35 40 45 50
time (s)
(c) scalogram percentage of energy for each wavelet coefficient
25
15.3846
9.46746
pseudo–frequency (Hz)

5.82613
3.58531
2.20634
1.35775
0.835538
0.514178
0.316417
0.194718
0.119827
0 5 10 15 20 25 30 35 40 45 50
seconds

Figure 15. The TFE plots of El Centro earthquake data using the (a) FDM, (b) EMD and (c) CWT. (Online version in colour.)

design is less than 10 Hz, and from the Fourier-based PSD, marginal spectrum by FDM and
marginal spectrum by EMD, one can observe that almost all the energy in this data, as shown
in figure 13, lies within 10 Hz. The IE fluctuations of the El Centro earthquake data estimated
by the FDM and EMD methods are shown in figure 14. The IE estimated by the EMD has some
transient and rapid fluctuations, whereas the IE estimated by the FDM appears smoother than
the EMD case. The TFE distribution by the FDM, EMD and CWT methods are shown in figure 15,
and all three methods indicate maximum energy concentration around 1.7 Hz and 2 s. There is
 (a) (b) 23
0.5 0.10

|X(f)|
I 0 0.05

[Link] Proc. R. Soc. A 473: 20160871

−0.5 0

...................................................
0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0 1000 2000 3000 4000 5000 6000 7000 8000
time (s) frequency (Hz)
×10−7
5 0.02
y1 y2
0 0
−5 −0.02

0.05 0.05
y3 0 y4 0
−0.05 −0.05
0.1 0.2
y5 0 y6 0
−0.1 −0.2
0.5 0.1
y7
Downloaded from [Link] on 11 August 2023

0 y8 0
−0.5 −0.1
×10−3 ×10−4
2 2
y9 0 y10 0
−2 −2
×10−4 2
2
y11 0 y12 0
−2 −2

Figure 16. Long vowel ‘I’ sound with sampling frequency Fs = 16000 Hz (a) and its Fourier spectrum (b). The FDM generates
highest frequency component (y1 ), FIBFs (y2 –y11 ) and mean-value component (y12 ). The fundamental frequency F0 is captured
accurately in FIBF y8 . (Online version in colour.)

enhanced TFE tracking by the FDM as it provides better and finer details of how the different
waves arrive from the epical centre to the recording station, for example the compression waves of
small amplitude but higher frequency range of 10–20 Hz, the shear and surface waves of strongest
amplitude and lower frequency range of below 5 Hz which does most of the damage, and other
body shear waves which are present over the full duration of the data span.

(h) Application to speech signal analysis

Finally, to illustrate the advantages of FDM in speech signal analysis, we decompose a quasi-
periodic voice signal (long vowel ‘I’). In figure 16, the FIBFs y2 –y5 seem to extract most of the
noise, while y6 to y7 seem to express most of the signal information at high frequencies. Moreover,
the FIBF y8 captures accurately F0 , the fundamental frequency (inverse of glottal pulse or pitch
period length) of the voice signal. By contrast, the EMD (figure 17a) presents mode mixing
between modes 2 and 3, between modes 3 and 4 and between modes 4 and 5. Also, the EMD
is not able to capture F0 in any mode. In figure 17b, the EEMD is able to capture F0 in mode
y5 , although it can be observed that appearance of positive minima and negative maxima in the
mode y4 and y7 , which is a violation of the first IMF condition. The EEMD modes y3 and y4 are a
asymmetric waveform which is a violation of the second IMF condition. Besides these problems,
the major drawbacks in this case are huge computation complexity for this ensemble size and
the reconstruction error that may be significant. The computation complexity can be reduced by
selecting the smaller ensemble size but that increases the reconstruction error. Figure 18 shows
the plots of TFE distribution obtained by the FDM, EMD and EEMD algorithms for same data
(long vowel ‘I’). Clearly, there is enhanced TFE tracking when using the FDM. The results in this
application as well are in clear favour of the Fourier method proposed in this paper.

6. Conclusion
In this paper, we have proposed (i) a novel and adaptive FDM for nonlinear and non-stationary
time series analysis, which decomposes any data into a small number of band-limited FIBFs. The
FDM is a generalized Fourier expansion with variable amplitudes and variable frequencies of a
time series by the Fourier method itself. (ii) The zero-phase filter bank-based MFDM algorithm,
 (a) EMD
0.2 0.5
24
y1 0 y2 0

[Link] Proc. R. Soc. A 473: 20160871

...................................................
−0.2 −0.5
0.5 0.2
y3 0 y4 0
−0.5 −0.2
0.05 0.02
y5 0 y6 0
−0.05 −0.02
0.02 0.01
y7 0 y8 0
−0.02 −0.01
×10–3 ×10–3
5 1
Downloaded from [Link] on 11 August 2023

y9 0 y10 0

−5 −1

×10–5
1 0.02
y11 0 y12 0.01
−1 0

(b) EEMD
0.1 0.05
y1 0 y2 0
−0.1 −0.05
0.5 0.1
y3 0 0
y4
−0.5 −0.1
0.1 0.02
y5 0 0
y6
−0.1 −0.02
×10–3 ×10–3
5 5
y7 0 y8 0
−5 −5
×10–3 ×10–3
5 2
y9 0 y10 0
−5 −2
×10–4
5 0.012
y11 0 y12 0.010

−5 0.008

Figure 17. IMFs (y1 –y12 ) of speech signal generated by (a) EMD and (b) EEMD (with zero mean and 0.2 standard deviation of the
added noise of normal distribution and the ensemble size of 300) algorithms for long vowel ‘I’ with Fs = 16 000 Hz. The EMD is
not able to catch F0 in any mode, and the EEMD is able to capture F0 accurately in IMF y5 .

for the analysis of multivariate nonlinear and non-stationary time series, which generates a
finite number of band limited multivariate FIBFs (MFIBFs). (iii) An algorithm to obtain cut-off
frequencies required in MFDM for zero-phase high or low pass filtering of multivariate signals.
The fundamental and conceptual contributions of this study are the Fourier based adaptive
signal decomposition method (FDM) and the introduction of the FIBFs. The FIBFs form the
basis of the decomposition that are complete, orthogonal, local and adaptive. The instantaneous
frequencies of the FIBFs yield TFE distribution of any signal. The TFE distribution of a signal
is used in various fields of science and engineering for analysis of physical phenomena and
engineering systems. The proposed methods produce the final presentation of the results as a TFE
distribution that reveals the imbedded structures of a signal. Unlike the various EMD algorithms,
the FDM and MFDM are mathematically well defined, supported by the well-established theories
of a zero-phase filter and the FTs. The FDM and MFDM do not suffer from the limitations
with regard to mode mixing, detrend uncertainty and end effect artefacts of EMD algorithms.
Simulation results clearly demonstrate the efficacy of the proposed methods.
 (a) FDM: time–frequency–energy estimate of FIBFs (HTL-FS)
25

6000

[Link] Proc. R. Soc. A 473: 20160871

...................................................
frequency (Hz) 5000

4000

3000

2000

1000
Downloaded from [Link] on 11 August 2023

0
0 0.05 0.10 0.15 0.20 0.25 0.30 0.35
time (s)
(b) EMD: time–frequency–energy estimate of IMFs
7000

6000

5000
frequency (Hz)

4000

3000

2000

1000

0
0 0.05 0.10 0.15 0.20 0.25 0.30 0.35
time (s)
(c) EEMD: time–frequency–energy estimate of IMFs

7000

6000
frequency (Hz)

5000

4000

3000

2000

1000

0
0 0.05 0.10 0.15 0.20 0.25 0.30 0.35
time (s)

Figure 18. The TFE plots of speech signal (long vowel ‘I’) by FDM (a), EMD (b) and (c) EEMD. There is enhanced TFE tracking
when using FDM and it captures F0 accurately in voice signal. (Online version in colour.)

Data accessibility. The Matlab software of the MEMD ([Link]

[Link]), EMD and EEMD ([Link] the proposed FDM ([Link]
[Link]/publication/274570245_Matlab_Code_Of_The_Fourier_Decomposition_Method_FDM) are
available online.
 Authors’ contributions. P.S., S.D.J., R.K.P. and K.S. designed the study and organized the entire research program.
26
P.S. wrote the manuscript, MATLAB code and performed the simulations. S.D.J., R.K.P. and K.S. corrected
and revised the manuscript, supervised the modelling and simulation. All authors discussed the results,

[Link] Proc. R. Soc. A 473: 20160871

...................................................
commented on the manuscript and approved it.
Competing interests. We declare we have no conflicts of interest.
Funding. The authors thank the Department of Electrical Engineering, IIT Delhi for funding extra page charges.
Acknowledgements. The authors thank JIIT Noida for permitting research to be carried out at IIT Delhi.

References
1. Hlawatsch F, Boudreaux-Bartels GF. 1992 Linear and quadratic time-frequency signal
representations. IEEE Signal Process. Mag. 9, 21–67. (doi:10.1109/79.127284)
Downloaded from [Link] on 11 August 2023

2. Priestley MB. 1965 Evolutionary spectra and non-stationary processes. J. R. Stat. Soc. B 27,
204–237.
3. Huang NE, Shen Z, Long S, Wu M, Shih H, Zheng Q, Yen N, Tung C, Liu H. 1988 The
empirical mode decomposition and Hilbert spectrum for non-linear and non-stationary time
series analysis. Proc. R. Soc. Lond. A 454, 903–995. (doi:10.1098/rspa.1998.0193)
4. Costa M, Priplata AA, Lipsitz LA, Wu Z, Huang NE, Goldberger AL, Peng CK. 2007 Noise and
poise: enhancement of postural complexity in the elderly with a stochastic-resonance-based
therapy. Europhys. Lett. EPL 77, 68008. (doi:10.1209/0295-5075/77/68008)
5. Cummings DA, Irizarry RA, Huang NE, Endy TP, Nisalak A, Ungchusak K, Burke DS. 2004
Travelling waves in the occurrence of dengue haemorrhagic fever in Thailand. Nature 427,
344–347. (doi:10.1038/nature02225)
6. Hu K, Peng CK, Huang NE, Wu Z, Lipsitz LA, Cavallerano J, Novaka V. 2008 Altered
phase interactions between spontaneous blood pressure and flow fluctuations in type 2
diabetes mellitus: nonlinear assessment of cerebral autoregulation. Phys. A 387, 2279–2292.
(doi:10.1016/[Link].2007.11.052)
7. Lo MT et al. 2013 A new method to estimate the amplitude spectrum analysis of
ventricular fibrillation during cardiopulmonary resuscitation. Resuscitation 84, 1505–1511.
(doi:10.1016/[Link].2013.07.004)
8. Huang NE, Wu Z. 2008 A review on Hilbert–Huang transform: method and its applications
to geophysical studies. Rev. Geophys. 46, RG2006. (doi:10.1029/2007RG000228)
9. Wu Z, Norden EH, Chen X. 2009 The multi-dimensional ensemble empirical mode
decomposition method. Adv. Adapt. Data Anal. 1, 339–372. (doi:10.1142/S1793536909000187)
10. Wu Z, Huang NE. 2009 Ensemble empirical mode decomposition: a noise-assisted data
analysis method. Adv. Adapt. Data Anal. 1, 1–41. (doi:10.1142/S1793536909000047)
11. Rehman N, Mandic DP. 2010 Multivariate empirical mode decomposition. Proc. R. Soc. A 466,
1291–1302. (doi:10.1098/rspa.2009.0502)
12. Chu PC, Fan C, Huang N. 2012 Compact empirical mode decomposition: an algorithm to
reduce mode mixing, end effect, and detrend uncertainty. Adv. Adapt. Data Anal. 4, 1250017.
(doi:10.1142/S1793536912500173)
13. Singh P, Srivastava PK, Patney RK, Joshi SD, Saha K. 2013 Nonpolynomial spline based
empirical mode decomposition. In 2013 Int. Conf. on Signal Processing and Communication,
pp. 435–440.
14. Singh P, Patney RK, Joshi SD, Saha K. 2014 Some studies on nonpolynomial interpolation and
error analysis. Appl. Math. Comput. 244, 809–821. (doi:10.1016/[Link].2014.07.049)
15. Singh P, Patney RK, Joshi SD, Saha K. 2015 The Hilbert spectrum and the energy preserving
empirical mode decomposition. ([Link] [[Link]])
16. Mandic DP, Rehman N, Wu Z, Huang NE. 2013 Empirical mode decomposition-based time-
frequency analysis of multivariate signals. IEEE Signal Process. Mag. 30, 74–86. (doi:10.1109/
MSP.2013.2267931)
17. Boashash B. 1992 Estimating and interpreting the instantaneous frequency of a signal. I.
Fundamentals. Proc. IEEE 80, 520–538. (doi:10.1109/5.135376)
18. Flandrin P, Rilling G, Goncalves P. 2004 Empirical mode decomposition as a filter bank. IEEE
Signal Process. Lett. 11, 112–114. (doi:10.1109/LSP.2003.821662)
19. Rehman N ur, Mandic DP. 2011 Filter bank property of multivariate empirical mode
decomposition. IEEE Trans. Signal Process. 59, 2421–2426. (doi:10.1109/TSP.2011.2106779)
 20. Daubechies I, Lu J, Wu H. 2011 Synchrosqueezed wavelet transforms: An empirical mode
27
decomposition-like tool. Appl. Comput. Harmonic Anal. 30, 243–261. (doi:10.1016/[Link].
2010.08.002)

[Link] Proc. R. Soc. A 473: 20160871

...................................................
21. Dragomiretskiy K, Zosso D. 2014 Variational mode decomposition. IEEE Trans. Signal Process.
62, 531–544. (doi:10.1109/TSP.2013.2288675)
22. Jain P, Pachori RB. 2015 An iterative approach for decomposition of multi-componentnon-
stationary signals based on eigenvalue decomposition of the Hankel matrix. J. Franklin Inst.
352, 4017–4044. (doi:10.1016/[Link].2015.05.038)
23. Hou TY, Shi Z. 2011 Adaptive data analysis via sparse time-frequency representation. Adv.
Adapt. Data Anal. 3, 1–28. (doi:10.1142/S1793536911000647)
24. Feldman M. 2006 Time-varying vibration decomposition and analysis based on the Hilbert
transform. J. Sound Vib. 295, 518–530. (doi:10.1016/[Link].2005.12.058)
25. Selesnick IW. 2011 Resonance-based signal decomposition: a new sparsity-enabled signal
Downloaded from [Link] on 11 August 2023

analysis method. Signal Process. 91, 2793–2809. (doi:10.1016/[Link].2010.10.018)

26. Gabor D. 1946 Theory of communication, part 1: the analysis of information. J. Inst. Electr.
Eng. III 93, 429–441. (doi:10.1049/ji-3-2.1946.0074)
27. Bedrosian E. 1963 A product theorem for Hilbert transform. Proc. IEEE 51, 868–869.
(doi:10.1109/PROC.1963.2308)
28. Nuttall AH, Bedrosian E. 1966 On the quadrature approximation to the Hilbert transform of
modulated signals. Proc. IEEE 54, 1458–1459. (doi:10.1109/PROC.1966.5138)
29. Boashash B. 1992 Estimating and interpreting the instantaneous frequency of a signal–part 2:
algorithms and applications. Proc. IEEE 80, 540–568. (doi:10.1109/5.135378)
30. Claasen TACM, Mecklenbrauker WFG. 1980 The Wigner distribution–a tool for time-
frequency signal analysis, Part II: Discrete time signals. Philips J. Res. 35, 276–300.
31. Boashash B. 1990 Time-frequency signal analysis. In Advances in spectrum analysis and array
process, vol. 1 (ed. S Haykin). New York, NY: Prentice Hall.
32. Huang NE, Wu Z. 2008 A review on Hilbert-Huang transform: method and its applications
to geophysical studies. Rev. Geophys. 46, 1–23. (doi:10.1029/2007RG000228)
33. Thomas YH, Zuoqiang S. 2011 Adaptive data analysis via sparse time-frequency
representation. Adv. Adapt. Data Anal. 3, 1–28. (doi:10.1142/S1793536911000647)
34. Huang NE, Chen X, Lo MT, Wu Z. 2011 On Hilbert spectral representation: a true time-
frequency representation for nonlinear and nonstationary data. Adv. Adapt. Data Anal. 63,
63–93. (doi:10.1142/S1793536911000659)