Variability Analysis of Observational Time Series: An Overview of the Decomposition Methods for Non-Stationary and Noisy Signals

Abstract

The analysis of observational data sequences in Geophysics consists of characterizing the underlying dynamics. An important preliminary step aims to analyze the variability related to the observed dynamic. The specific objectives related to this step are to remove noise, to determine the overall trend of the observational time series and to identify the relevant components contributing significantly to the original time series variability knowing that their number determines the dimensionality of the observed dynamics. Most of the observational time series have characteristics of non-stationarity and present fluctuations at all-time scales. In this context, variability analysis consists in representing time series in the time-frequency space and requires the development of specific numerical signal decomposition methods. The most commonly used techniques are adaptive and data-driven and among the most cited in the literature are the empirical mode decomposition, the empirical wavelet transform, and singular spectrum analysis. In this work, we describe all of these techniques and evaluate their ability to remove noise and to identify components corresponding to the physical processes involved in the evolution of the observed system and deduce the dimensionality of the associated dynamics. Results obtained with all of these methods on experimental total ozone columns and rainfall time series will be discussed and compared.

1. Introduction

The notion of complexity is a characteristic present in most of systems of the real world. Intuitively, complexity constitutes a measure of the organization, the structure and the dynamics of a system i.e. its behavior over time. In the physical sciences, a complex system can be defined as composed of several entities interacting with each other on several scales of time and space. The complexity is created by the superimposition of the individual behavior of each entity with their interactions and makes the overall behavior of the system difficult to characterize and predict. Therefore, the characterization of the system dynamic is carried out from the global properties of the system which cannot be reduced to the sum of the characteristics of the entities of which it is composed. Many studies that aim to gain a better understanding of the complex system dynamics are based on the analysis of observational measurements data sequence. Most of the observational data sequences in geophysics derive from complex dynamics resulting from the superimposition of the individual behaviors of the physical processes involved in the evolution of the overall system with the interactions they have between them. As a consequence, time series are non-stationary and have fluctuations at all scales of time. Then, the analysis methodology consists of analyzing the variability of the time series. The method used is to decompose the observed time into the sum of noise, physically meaning full components and trends.

$$\mathit{TS}=\mathit{\text{noise}}+{\sum }_i{c}_i+\mathit{\text{Trend}},\mathit{TS}=\left\{x_k=x\left(\mathit{k\delta t}\right),k=1,N\right\}$$E1

Where TS is the original times series composed of N measurement values x_k of the signal x with a sample time δt and c_i are the components that contribute significantly to the variability of the original signal. Each of c_i can be identified with a physical process involved in the overall evolution of the system. To implement such decomposition, two approaches have been considered: The first approach consists in getting a representation of TS in the time-frequency space which amounts to classify the fluctuations according to the temporal resolution at which they occur. Among the most commonly used techniques enabling to obtain a representation of time series in time frequency space, the empirical mode decomposition (EMD) is very popular. EMD is a relatively new method proposed by Huang et al. in 1998 well suited to non-stationary signals that decompose a time series into a finite number of components called “Intrinsic Mode Function” (IMFs) [1]. IMF captures the repeating behavior of the signal at some particular time scales. Unlike the Fourier or wavelet transforms, EMD enables to decompose a time signal into a set of basis signals derived from the data itself. The biggest advantage of the EMD method is to be totally adaptive and data-driven without the need for a-priori basis function selection for signal decomposition. As EMD acts in the time domain, another advantage of this method is to be close to the observed dynamics. Moreover, the EMD acts as a bank of bandpass filter [2] in the frequency domain and as a result, the main limiting factor of EMD is the frequency resolution which when it becomes too small can induce the mode-mixing phenomenon where the spectral content of some IMFs overlaps each other [3]. Although several techniques exist to overcome this problem [4, 5, 6], Gilles proposed an alternative entitled empirical wavelet transform (EWT) [7], which operates in the frequency domain and consists of elaborating a segmentation of the original signal’s Fourier spectrum and building an appropriate wavelet filter bank on the segmented Fourier spectrum. The EWT allows a better frequency resolution and thus overcomes the mode-mixing problem by partitioning the spectrum of the original signal into separate spectral bands. However, although the EWT technique enables to detect the relevant frequencies involved in the original signal, it does not allow to associate to a specific mode of variability as EMD does. Because the variability modes provided by the EMD technique are closer to the observed dynamics than EWT, O. Delage [8] proposed an optimization of the EMD technique entitled Empirical Adaptive Wavelet Decomposition (EAWD) which combines the advantages of both EMD and EWT techniques. The heart of the EWT method lies in the segmentation of the original signal’s Fourier spectrum. The main idea is to use the local maxima of the IMFs’ spectral support returned by EMD to design the Fourier spectrum segmentation required to initiate the EWT method. IMFs involved in the Fourier spectrum segmentation are those selected to be relevant i.e. whose contribute more than 1% of the variability of the original signal. The final step consists in applying the EWT method on the segmented Fourier spectrum.

Another technique called Singular Spectrum Analysis (SSA) using a different approach is described in the scope of this chapter. Such a technique can decompose a non-stationary signal into a sum of independent components and also be able to reconstruct the underlying dynamics.

This document is structured around three sections: the first section is devoted to the description of all the signal decomposition techniques mentioned above. In the second section, three observational time series are analyzed by using EMD, EAWD and SSA techniques and the corresponding results are presented. In the last section, the results obtained respectively with the EMD, EAWD and SSA techniques are compared and discussed.

2. Review of the signal decomposition techniques

The methodology related to the signal variability analysis is composed of three main steps: (1) express the original signal as a sum of a finite number of components, (2) identify the components dominated by noise, (3) select the relevant components that contribute significantly to the signal variability with the objectives of determining the dimensionality of the underlying dynamics and identifying the relevant components to physical processes involved in the underlying dynamics.

2.1 The empirical mode decomposition (EMD)

2.1.1 EMD basics

In 1998, Huang et al. proposed an original method called the Empirical Mode Decomposition (EMD) whose purpose is to adaptively decompose any signals into oscillatory contributions. The EMD technique can be summarized as an iterative method where the signal can be decomposed into a local average m called trend and a strongly oscillating part called details. The trend is related to low frequencies while the details characterized by strong fluctuations are related to high frequencies. At each iteration, the high-frequency fluctuations part is separated from the low-frequencies trend and are reinjected as a new signal in the next iteration. During the iterations of the algorithm [1], the “details” related to the high frequencies are successively separated from the low-frequency part by using a procedure called “Sifting process”. The main steps of the sifting process consist in: identify the signal local extrema; the local maxima and minima are then interpolated by using cubic splines to form respectively the upper and lower envelopes of the signal; the mean envelope is then determined by calculating the half sum of the upper and lower envelopes; the mean envelope is subtracted from the initial signal. The same procedure is reiterated on the remainder until the mean envelope is close to zero everywhere and the resultant signal is designated as the first IMF. The criterion for the sifting process to stop is generally set as:

$$0.2\le \mathit{SD}=\sum _{\mathit{i}=1}^{\mathit{T}}\left[\frac{{\left|{\mathit{h}}_{\mathit{k}-1}\left(\mathit{t}\right)-{\mathit{h}}_{\mathit{k}}\left(\mathit{t}\right)\right|}^2}{{{\mathit{h}}_{\mathit{k}-1}}^2}\right]\le 0.3$$E2

Where h_k is the result after k iterations of the sifting process. The higher order IMFs are iteratively extracted by using the above sifting process procedure until the remaining signal cannot be assimilated to an IMF knowing that a signal can be called IMF if it satisfies the two following criteria: the numbers of local maxima and local minima must differ by at most one; the half sum of its upper and lower envelopes is locally close to zero. The original signal can be finally expressed as:

$$\mathit{TS}=\sum _{\mathit{m}=1}^{\mathit{M}}{\mathit{IMF}}_{\mathit{m}}+\mathit{R}$$E3

Where R is the residual trend. EMD decomposes the data into M fundamental components with distinct time scales where the first component has the smallest time scale. The interesting fact about this algorithm is that it is adaptive and data-driven and is able to extract non-stationary parts of the original signal at different time scales as multi-resolution analysis does. In this context, Flandrin et al. [9] described the EMD as behaving as a dyadic filter bank as those involved in the multi-resolution analysis. Consequently, the maximum number of IMFs that can be extracted from the original time series is Log(N)/Log (2) if N is the time series size.

2.1.2 Noise removing-Noisy components identification

Generally, the original signal is corrupted by noise. As a consequence, a few IMFs may be the oscillations of the noise-free signal and the others correspond to noise. To determine which IMFs are noise-related components, a robust threshold is required. In this chapter, a “Detrended Fluctuations analysis” (DFA) technique is used to obtain such a threshold. The basic principle of the DFA [10, 11] is to compute how the time series fluctuations around the local trend varies as a function of the time scale. When applied to a time series TS(t), the first step of the DFA technique is to compute the time series TSI(t) composed of the cumulative sums of TS after removing its mean:

$$\mathit{TSI}\left(\mathit{k}\right)=\sum _{\mathit{i}=1}^{\mathit{k}}\left[\mathit{x}\left(\mathit{i\delta t}\right)-\overline{\mathit{TS}}\right],1\le \mathit{k}\le \mathit{N}$$E4

Where $\overline{\mathit{TS}}$ states the average of TS time series over [1, N], TSI is then divided into time windows of size n samples each. For each time window, the estimated local trend TSI_n(k) is determined by using least-square linear fitting. Finally, the average root mean square (RMS) of the fluctuations at a specific time scale n, F(n) may be written as:

$$\mathit{F}\left(\mathit{n}\right)=\sqrt{\frac{1}{\mathit{N}}{\sum }_{\mathit{k}=1}^{\mathit{N}}{\left[\mathit{TSI}\left(\mathit{k}\right)-{\mathit{TSI}}_{\mathit{n}}\left(\mathit{k}\right)\right]}^2}$$E5

The method suggested to identify the noise-related IMFs is to use the slope α of the curve log[F(n)]/log(n) as a threshold. The method consists in excluding IMFs with α value less than a threshold θ. For DFA, an α value of 0.5 characterizes an uncorrelated white noise. The commonly threshold value θ taken in the literature is 0.5 with a 0.2 confidence interval that is θ = 0.7.

2.1.3 Relevant components selection

Although the EMD is a powerful tool for analyzing complicated datasets, many irrelevant IMF may appear in the decomposition. A relevant IMF may be defined as an IMF that retains most of the information content of the signal. As a consequence, relevant IMFs would have a good correlation with the original signal while irrelevant ones would have a poor correlation. For discriminating between relevant and irrelevant IMFs, we use the Pearson correlation coefficient between each IMF and the original time series, i.e.:

$$\mathit{\text{CORR}}\left({\mathit{IMF}}_{\mathit{m}},\mathit{TS}\right)=\frac{\mathit{cov}\left({\mathit{IMF}}_{\mathit{m}},\mathit{TS}\right)}{\mathit{sd}\left({\mathit{IMF}}_{\mathit{m}}\right)\ast \mathit{sd}\left(\mathit{TS}\right)}\mathrm{m}=1,\mathrm{M}$$E6

Where cov(IMFm,TS) is the covariance between the m^th IMF and the original time series, sd(IMF_m) is the standard deviation of the m^th IMF while sd(TS) is the standard deviation of the original time series. A threshold s is required for selecting relevant IMF: if CORR(IMF_m,TS) > s, keep the m^th IMF. Otherwise, eliminate the m^th IMF and add it to the residual. As the threshold is different for different signals, a suitable threshold must be selected according to the relativity between IMfs and the original signal. In general threshold s can be the ratio of the mean value between the maximum and the minimum of $\left\{\mathit{\text{CORR}}\left({\mathit{IMF}}_{\mathit{m}},\mathit{TS}\right)\mathit{m}\mathit{=}1,\mathit{M}\right\}$ that is:

$$\mathit{s}=\frac{\mathit{M}+\mathit{m}}{2\mathit{\eta }}$$E7

Where η is a coefficient greater than 1, M = max_m(CORR(IMF_m,TS), m = min_m(CORR(IMF_m,TS).

So that the correlation coefficients between each of the IMFs and the original signal would be in the interval [0,1], all IMF and original signals will be normalized at first. In the literature, two different thresholds are proposed ${\mathit{s}}_1\mathit{=}\frac{\mathit{M}\mathit{+}\mathit{m}}{20}$ [12] and ${\mathit{s}}_2\mathit{=}\frac{\mathit{M}}{10\mathit{M}\mathit{-}3}$ [13]. The threshold proposed in this work depends on the relative position de M/m and 10. If we suppose that 0.7 ≤ M ≤ 1, and M/m > 10, m could be less than 0.07 which characterizes a very weak correlation. In this case s has to be greater than m and from (Eq. (7)), η > 5.5 which is verified when, s = s₁. If on the contrary when M/m ≤ 10, m could be greater than 0.1 and s must be close to m which from (Eq. (7)) can be written:

$$1-\mathit{\epsilon }\le \frac{\mathit{s}}{\mathit{m}}=\frac{1+\frac{\mathit{M}}{\mathit{m}}}{2}.\frac{1}{\mathit{\eta }}\le 1+\mathit{\epsilon }$$E8

Given that 5 ≤ M/m ≤ 10, then 3 ≤ (1 + M/m)/2 ≤ 5.5 and (Eq. (8)) becomes

$$\frac{5.5}{1+\mathit{\epsilon }}\le \mathit{\eta }\le \frac{3}{1-\mathit{\epsilon }}$$E9

If we set ε = 0.052, (Eq. (9)) becomes 3 ≤ η ≤ 5.2 which is verified when s = s₂.

In summary, s is defined as follows:

$$\mathit{s}=\left\{\begin{array}{l}\frac{\mathit{M}+\mathit{m}}{20}\mathit{\text{if}}\frac{\mathit{M}}{\mathit{m}}>10\\ \frac{\mathit{M}}{10\mathit{M}-3}\mathit{\text{if}}\frac{\mathit{M}}{\mathit{m}}\le 10\end{array}\right.$$E10

2.1.4 EMD limiting factors-IMF oscillation cycle

As the EMD technique acts as a bank of bandpass filters [2], each IMF is associated with a specific bandpass and the oscillation cycle corresponds to the prevailing frequency in the bandpass. One way to get the IMFs oscillation cycle is to calculate the spectral density for each IMF and identify the frequency for which the spectral density is maximum. The local maxima occurring in the IMFs spectrum characterize the frequencies involved in the corresponding EMD mode of variability. One of the major limiting factors of the EMD is the frequency resolution which can, when it is not sufficient, induce the mode-mixing phenomenon where the spectral content of some IMFs overlaps each other. The source at the origin of mode-mixing can be categorized in three main groups: (1) presence of noise, (2) presence of intermittency, (3) presence of closely spaced spectral components. To overcome the inherent mode-mixing problem of the EMD, Wu and Huang [14] proposed the Ensemble Empirical Mode Decomposition (EEMD). In EEMD, white noise signals n_i(t) are added to the original signal x(t); Because the white noise spectrum is evenly distributed, the white noise signals will be automatically distributed to the appropriate reference scales. Moreover, because of its zero-mean characteristic, the white noise will cancel itself out after many rounds of averaging. The specific steps of the EEMD algorithm are as follows:

• Initialize the number of ensemble members M. - Compute M realizations of white noise with the different variances that are added to the original signal: ${\mathit{x}}_{\mathit{i}}\left(\mathit{t}\right)\mathit{=}{\mathit{n}}_{\mathit{i}}\left(\mathit{t},{\mathit{\sigma }}_{\mathit{i}}\right)\mathit{+}\mathit{x}\left(\mathit{t}\right)\mathit{i}\mathit{=}1\mathit{,}\mathit{M}$ where n_i is a white noise of variance σ_i. – Use the EMD algorithm to decompose x_i(t) into IMFs: ${\mathit{x}}_{\mathit{i}}\left(\mathit{t}\right)\mathit{=}{\sum }_{\mathit{j}\mathit{=}1}^{\mathit{J}}{\mathit{IMF}}_{\mathit{i}\mathit{,}\mathit{j}}\mathit{+}{\mathit{R}}_{\mathit{i}}$ i = 1,M - Calculate the ensemble means of the decomposed IMFs: ${\mathit{IMF}}_{\mathit{j}}\left(\mathit{t}\right)\mathit{=}\frac{1}{\mathit{M}}{\sum }_{\mathit{i}\mathit{=}1}^{\mathit{M}}{\mathit{IMF}}_{\mathit{i}\mathit{,}\mathit{j}}$

2.2 Wavelet approaches

Wavelets are commonly used to analyze the variability of a signal. In the temporal domain, a wavelet basis is defined as the mother wavelet Ψ of zero mean, dilated with a parameter s > 0 and translated by u∈ℜ:

$${\mathit{\psi }}_{\mathit{u},\mathit{s}}\left(\mathit{t}\right)=\frac{1}{\sqrt{\mathit{s}}}\mathit{\psi }\left(\frac{\mathit{t}-\mathit{u}}{\mathit{s}}\right)$$E11

For the wavelet decomposition of a time series x(t), the most widely used case is the dyadic one, s = 2^j. Then the wavelet decomposition of x is obtained by computing the inner product W_x(kj) as:

$${\mathit{W}}_{\mathit{x}}\left(\mathit{k},\mathit{j}\right)=⟨\mathit{x},{\mathit{\psi }}_{\mathit{k},\mathit{j}}⟩\mathit{\text{with}}{\mathit{\psi }}_{\mathit{k},\mathit{j}}\left(\mathit{t}\right)=\frac{1}{\sqrt{2^{\mathit{j}}}}\mathit{\psi }\left(\frac{\mathit{t}-\mathit{k}}{2^{\mathit{j}}}\right),\mathit{k}\in \mathit{Z}$$E12

where j represents the resolution level. The decomposition is then similar to a multi-resolution analysis carrying out successive projections of x on a sequence of nested subspaces V_j∈L²(ℜ) j = [0, n], which leads to increasingly coarse approximations of x as j increases. At each resolution level j, the information about the signal x, at resolution j-1 in the subspace V_j-1, is split into two parts named approximation and details corresponding, respectively to the low frequencies and the high frequencies contained in the signal x. The approximation of x results from the orthogonal projection of x onto the subspace V_j and the information of details results from the orthogonal projection of x onto the subspace W_j orthogonal to V_j such that ${\mathit{V}}_{\mathit{j}\mathit{-}1}\mathit{=}{\mathit{V}}_{\mathit{j}}\mathit{⨁}{\mathit{W}}_{\mathit{j}}$, where ⊕ denotes the direct sum of vector subspaces. Wavelets $\left\{{\mathit{\psi }}_{\mathit{k}\mathit{,}\mathit{j}}\left(\mathit{t}\right),\mathit{k}\mathit{\in }\mathit{Z}\right\}$ form the basis of W_j. According to the definition of multi-resolution analysis, there exists a function φ(t), called scaling function, such that $\left\{\mathit{\phi }\left(\mathit{t}\mathit{-}\mathit{k}\right),\mathit{k}\mathit{\in }\mathit{Z}\right\}$ form a basis of V₀ corresponding to the coarsest approximation of x. The reconstruction of x is obtained from:

$$\mathit{x}\left(\mathit{t}\right)=⟨\mathit{x}\left(\mathit{t}\right),\mathit{\phi }\left(\mathit{t}\right)⟩.\mathit{\phi }\left(\mathit{t}\right)+\sum _{\mathit{j}=1}^{\mathit{n}}⟨\mathit{x}\left(\mathit{t}\right),{\mathit{\psi }}_{\mathit{k},\mathit{j}}\left(\mathit{t}\right)⟩.{\mathit{\psi }}_{\mathit{k},\mathit{j}}\left(\mathit{t}\right)$$E13

where $⟨⟩$ represents the inner product. The approximation coefficients corresponding to the coarsest resolution level are given by $⟨\mathit{x},\mathit{\phi }⟩$ and the details coefficients corresponding to the successively decreasing resolution level $\mathit{∆}{\mathit{s}}_{\mathit{j}}\mathit{=}\frac{\mathit{∆}{\mathit{s}}_{\mathit{j}\mathit{-}1}}{2}$ are given by $⟨\mathit{x},{\mathit{\psi }}_{\mathit{k}\mathit{,}\mathit{j}}⟩$ as follows:

$$⟨\mathit{x}\left(\mathit{t}\right),\mathit{\psi }\left(\mathit{t}\right)⟩=\int \mathit{x}\left(\mathit{\tau }\right).\overline{\mathit{\phi }\left(\mathit{\tau }-\mathit{t}\right)}\mathit{d\tau },⟨\mathit{x}\left(\mathit{t}\right),{\mathit{\psi }}_{\mathit{k},\mathit{j}}\left(\mathit{t}\right)⟩=\int \mathit{x}\left(\mathit{\tau }\right).\overline{{\mathit{\psi }}_{\mathit{k},\mathit{j}}\left(\mathit{\tau }-\mathit{t}\right)}\mathit{d\tau }$$E14

2.2.1 Empirical wavelets-EWT

The essence of EMD is that the time domain functions into which a signal is decomposed have the same length as the original signal, allowing time-varying frequencies to be preserved. EMD is described [2, 9] as behaving as a dyadic filter bank as those involved in the multi-resolution analysis. In this context, the mode-mixing phenomenon specific to EMD can be interpreted as the presence of several filters of overlapping frequency content. As a result, the spectral content of some IMFs restituted by the EMD overlap each other. To overcome this problem, Gilles [7] proposed an alternative named the EWT. This method acts in the spectral domain and starts from the segmentation of the original signal’s Fourier spectrum. The Fourier support [0,π] is subdivided into N non-overlapping contiguous segments denoted Δ_n = [ω_n-1,ω_n]. An appropriate wavelet filter bank is then used to extract spectra relative to each Fourier segment. In the time domain, the components related to the original signal’s decomposition are obtained by performing inverse Fourier transform on each Fourier spectrum segments.

The filter bank [15,16] is defined by the empirical scaling function and the empirical wavelets on each Δ_n through the following equations:

$$\widehat{{\mathit{\phi }}_{\mathit{n}}\left(\mathit{\omega }\right)}=\left\{\begin{array}{l}\begin{array}{l}\begin{array}{l}1\mathit{\text{if}}⌈\mathit{\omega }⌉\le \left(1-\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}}\\ \mathit{cos}\left[\frac{\mathit{\pi }}{2}\mathit{\beta }\left(\frac{1}{2\mathit{\gamma }{\mathit{\omega }}_{\mathit{n}}}\left(\left|\mathit{\omega }\right|-\left(1-\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}}\right)\right)\right]\end{array}\\ \mathit{\text{if}}\left(1-\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}}\le \left|\mathit{\omega }\right|\le \left(1+\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}}\end{array}\\ 0\mathit{\text{otherwise}}\end{array}\right.$$E15

and

$$\widehat{{\mathit{\psi }}_{\mathit{n}}\left(\mathit{\omega }\right)}=\left\{\begin{array}{l}1\mathit{\text{if}}\left(1+\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}}\le \left|\mathit{\omega }\right|\le \left(1-\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}+1}\\ \mathit{cos}\left[\frac{\mathit{\pi }}{2}\mathit{\beta }\left(\frac{1}{2\mathit{\gamma }{\mathit{\omega }}_{\mathit{n}+1}}\left(\left|\mathit{\omega }\right|-\left(1-\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}+1}\right)\right)\right]\\ \mathit{\text{if}}\left(1-\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}+1}\le \left|\mathit{\omega }\right|\le \left(1+\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}+1}\\ \mathit{sin}\left[\frac{\mathit{\pi }}{2}\mathit{\beta }\left(\frac{1}{2\mathit{\gamma }{\mathit{\omega }}_{\mathit{n}}}\left(\left|\mathit{\omega }\right|-\left(1-\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}}\right)\right)\right]\\ \mathit{\text{if}}\left(1-\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}}\le \left|\mathit{\omega }\right|\le \left(1+\mathit{\gamma }\right){\mathit{\omega }}_{\mathit{n}}\\ 0\mathit{\text{otherwise}}\end{array}\right.$$E16

The function β(x) is an arbitrary C^k([0,1]) function defined as:

$$\mathit{\beta }\left(\mathit{x}\right)=\left\{\begin{array}{l}0\mathit{ifx}\le 0\\ \mathit{\beta }\left(\mathit{x}\right)+\mathit{\beta }\left(1-\mathit{x}\right)=1\mathit{x}\in \left[0,1\right]\\ 1\mathit{ifx}\ge 1\end{array}\right.$$E17

Many functions satisfy this property and the one the most used in the literature [17] is:

$$\mathit{\beta }\left(\mathit{x}\right)={\mathit{x}}^4\left(35-84\mathit{x}+70{\mathit{x}}^2-20{\mathit{x}}^3\right)$$E18

The parameter γ is chosen to satisfy the following criterion:

$$\mathit{\gamma }<{\mathit{Min}}_{\mathit{n}}\left(\frac{{\mathit{\omega }}_{\mathit{n}+1}-{\mathit{\omega }}_{\mathit{n}}}{{\mathit{\omega }}_{\mathit{n}+1}+{\mathit{\omega }}_{\mathit{n}}}\right)$$E19

The details and approximation coefficients are calculated by using (Eqs. (15) and (16)) and are respectively given by inner products with the empirical wavelets ψ_n and the scaling function φ₁:

$${\mathit{W}}_{\mathit{x}}\left(\mathit{n},\mathit{t}\right)=⟨\mathit{x}\left(\mathit{t}\right),{\mathit{\psi }}_{\mathit{n}}⟩=\mathit{\text{IFFT}}$$E20

$${\mathit{W}}_{\mathit{x}}\left(1,\mathit{t}\right)=⟨\mathit{x}\left(\mathit{t}\right),{\mathit{\phi }}_1⟩=\mathit{\text{IFFT}}$$E21

where X is the Fourier transform of the original signal x, $\overline{}$ represents the complex conjugate, IFFT represents the inverse Fourier transform, and ψ_n and φ₁ are the results of the inverse Fourier transforms of $\widehat{{\mathit{\psi }}_{\mathit{n}}}$ and $\widehat{{\mathit{\phi }}_1}$ respectively. The segmentation of the original signal’s Fourier spectrum is obtained from the detection of local maxima. Each segment is centered around one or a group of successive local maxima. The boundary between two contiguous segments is set as the nearest local minimum to the midpoint between two adjacent local maxima groups. Many of the detected local maxima are irrelevant as their contributions to the original signal variability are negligible. Selecting the relevant local maxima requires determining a threshold which is not always possible.

2.2.2 The empirical adaptive wavelet decomposition (EAWD)

The EMD technic enables an observational data sequence to be decomposed into multiple empirical modes of variability, each of them reflecting the observed dynamics at a specific timescale. As the spectral contents of the IMFs returned by the EMD are determined from the relative positions of the original signal’s maxima, one of the major drawbacks of the EMD is the mode-mixing phenomenon resulting in a few IMFs to have overlapping spectral supports. On the contrary, the EWT has a more solid theoretical context using wavelets and therefore provides a better frequency resolution. On the other hand, the EWT does not allow to associate the detected frequencies to a specific mode of variability as the EMD technique does. The main idea of the proposed EAWD method is to combine EMD and EWT techniques by setting non-overlapping groups of local maxima from the spectral contents of the IMFs returned by the EMD technique. Each IMF local maximum group will be associated with a segment of the original signal Fourier spectrum segmentation. The boundaries of each of these segments will be set as the local minima the closest to the midpoint between local maximum groups of two consecutive IMFs. As the EMD method acts as a bank of dyadic band-pass filters, the result of each of these filters in the frequency domain is composed of a set of local maxima relative to a specific time scale in which the resolution is divided by 2 in comparison with the timescale immediately above it. Considering that a time scale is characterized by a set of values in the range [2ⁿ,2^n + 1], to carry out a segmentation of the Fourier spectrum, it is necessary to distribute the local maxima groups relatively to a grid [2ⁱ,2^i + 1], i∈[2,J] with $\mathit{J}\mathit{=}\int \left(\frac{\mathit{log}\left(\mathit{N}\right)}{\mathit{log}\left(2\right)}\right)\mathit{-}1$, where N represents the size of the original time series and int. is the integer part of a number. The proposed Fourier spectrum segmentation algorithm has three main steps: (1) calculation of the spectral content of each IMF based on spectral density and selection of significant local maxima whose energy contribution is >1%. (2) Compute the cutoff frequency between the spectral supports of two consecutive IMFs. (3) EWT technic is run from the Fourier transform obtained from step 2. To get more information on the segmentation algorithm, it is described in detail in Ref. [8].

2.3 Signal decomposition based on singular spectrum analysis

In time series analysis, singular spectrum analysis (SSA) combines elements of classical time series analysis, multivariate statistics and dynamical systems. SSA is used both to decompose time series into components each having a meaningful interpretation and to reconstruct the underlying dynamics from a single time series based on the embedding theorem [18]. In this paragraph, as the reconstruction of the underlying dynamics is beyond the scope of this work, we will describe the SSA technique and its use in signal decomposition.

2.3.1 Singular spectrum decomposition

2.3.1.1 The basic singular spectrum analysis method

As EMD and EWT, SSA is a data-adaptive and non-parametric method for time series decomposition which is suitable to non-stationary time series. The efficiency of such a technic has been recognized for its ability to provide meaningful results in a wide range of application fields, without making any assumption on the processed data. Generally, the components extracted by SSA can be identified as trends, periodic (possibly amplitude-modulated) components or noise components. Furthermore, the characteristics of the components provided by SSA which are to be independent make it particularly suitable for blind source separation. SSA consists of four stages: embedding, decomposition, grouping, and reconstruction [19, 20, 21].

Embedding: The starting point of SSA is to embed a time series TS = {x_k, k = [1,N]} of size N in a vector space of dimension L. L is a strictly positive integer named the embedding dimension or window length, 1 < L < N. The embedding procedure forms K=N-L + 1 lagged vectors X_i = (x_i, ….x_i + L-1)^T with i = 1….N-L + 1. The trajectory matrix of the time series TS is then given by:

$$\mathit{MT}=\left[{\mathit{X}}_1,\dots {\mathit{X}}_{\mathit{K}}\right]$$E22

MT is a Hankel matrix (constant cross-diagonals) of size (LxK). Hence, the embedding procedure builds a sequence of L-dimensional vectors from the original time series TS, by using a sliding window of size L.

Decomposition: The singular value decomposition (SVD) of the trajectory matrix MT is then computed, providing MT = U.D.V^T with U = (LxL) and V = (KxK) being orthogonal matrices containing respectively the left and right singular vectors. D = (LxK) is a matrix containing the singular values σ_I on the main diagonal and zero elsewhere (where ${\mathit{\sigma }}_{\mathit{i}}\mathit{=}\sqrt{\mathit{K}{\mathit{\lambda }}_{\mathit{i}}}$) with λ_I being the eigenvalues of C the covariance matrix of MT which can be expressed $\mathit{C}\mathit{=}\frac{1}{\mathit{K}}\mathit{MT}\mathit{.}{\mathit{MT}}^{\mathit{T}}$. The matrix MT is thus decomposed in a sum of rank-one matrices MT_i such that:

$$\mathit{MT}=\sum _{\mathit{i}=1}^{\mathit{L}}{\mathit{MT}}_{\mathit{i}}=\sum _{\mathit{i}=1}^{\mathit{L}}{\mathit{\sigma }}_{\mathit{i}}{\mathit{u}}_{\mathit{i}}{\mathit{v}}_{\mathit{i}}^{\mathit{T}}$$E23

where u_i is i^th column of matrix U, v_i is the i^th column of matrix V. In other word, the SVD of the trajectory matrix corresponds to an orthogonal matrix whose column vectors form an eigenbasis of the multidimensional space created by the embedding step. Relevant components of the signal can be obtained by projecting the data onto that eigenbasis and omitting noise-related components and finally reconstruct an improved (denoised) version of the original time series.

Grouping step: This step consists of splitting the set of elementary matrices MT_i into r disjoint groups and summing the matrices within each group. The result may be written as:

$$\mathit{MT}=\sum _{\mathit{k}=1}^{\mathit{r}}{\mathit{MT}}_{{\mathit{I}}_{\mathit{k}}}\mathit{\text{with}}{\mathit{MT}}_{{\mathit{I}}_{\mathit{k}}}={\sum }_{\mathit{i}\in {\mathit{I}}_{\mathit{k}}}{\mathit{MT}}_{\mathit{i}}$$E24

Thus ${\mathit{MT}}_{{\mathit{I}}_{\mathit{k}}}$ is the resulting matrix of the group I_k, k = 1, …., r.

Diagonal averaging: The averaging along cross-diagonals of the matrix ${\mathit{MT}}_{{\mathit{I}}_{\mathit{k}}}$ aims at solving the problem of finding the time series TS^k for which the trajectory matrix is the closest to ${\mathit{MT}}_{{\mathit{I}}_{\mathit{k}}}$ in a least-squares sense. In other words, the diagonal averaging of ${\mathit{MT}}_{{\mathit{I}}_{\mathit{k}}}$={x_i,j, i = 1..L,j = 1…K} provide the elements of a time series $\left\{{\mathit{TS}}_{\mathit{n}}^{\mathit{k}},\mathit{n}\mathit{=}1,\dots \mathit{N}\right\}$ as:

$${\mathit{TS}}_{\mathit{n}}^{\mathit{k}}=\left\{\begin{array}{l}\frac{1}{\mathit{n}}\sum _{\mathit{m}=1}^{\mathit{n}}{\mathit{x}}_{\mathit{m},\mathit{n}-\mathit{m}+1}\mathit{\text{for}}1\le \mathit{n}<\mathit{L}\\ \frac{1}{\mathit{L}}\sum _{\mathit{m}=1}^{\mathit{L}}{\mathit{x}}_{\mathit{m},\mathit{n}-\mathit{m}+1}\mathit{\text{forL}}\le \mathit{n}\le \mathit{K}\\ \frac{1}{\mathit{N}-\mathit{n}+1}\sum _{\mathit{m}=\mathit{n}-\mathit{K}+1}^{\mathit{L}}{\mathit{x}}_{\mathit{m},\mathit{n}-\mathit{m}+1}\mathit{\text{forK}}+1\le \mathit{n}\le \mathit{N}\end{array}\right.$$E25

This cross-averaging process can also be applied to each MT_i matrix. The resulting time series can be assimilated to elementary components. This process finally provides an exact expansion of the original time series TS into r components that satisfies: $\mathit{TS}\mathit{=}{\sum }_{\mathit{k}\mathit{=}1}^{\mathit{r}}{\mathit{TS}}_{\mathit{n}}^{\mathit{k}}$.

2.3.1.2 Separability and choice of parameters

The very important question is how to choose parameters to build the proper decomposition of the observed time series. We need first to study the concept of separability which characterizes how well different components can be separated from each other. To do so, the so-called w-correlation matrix is computed. This is the matrix consisting of a weighted correlation between the reconstructed time series components. The weight reflects the number of entries of the time series terms into its trajectory matrix. Well separated components have a weak correlation while badly separated components have a large correlation. Therefore, looking at the w-correlation matrix, one can find groups of correlated elementary reconstructed series and use this information for the consequent grouping procedure. MATWCORR, the w-correlation matrix is expressed as follows:

$${\mathit{\text{MATWCORR}}}_{{\mathit{k}}_1,{\mathit{k}}_2}=\frac{\left({\mathit{TS}}_{{\mathit{k}}_1},{\mathit{TS}}_{{\mathit{k}}_2}\right)}{{‖{\mathit{TS}}_{{\mathit{k}}_1}‖}_{\mathit{w}}.{‖{\mathit{TS}}_{{\mathit{k}}_2}‖}_{\mathit{w}}}$$E26

where

$${‖{\mathit{TS}}_{{\mathit{k}}_{\mathit{i}}}‖}_{\mathit{w}}=\sqrt{{\left({\mathit{TS}}_{{\mathit{k}}_{\mathit{i}}},{\mathit{TS}}_{{\mathit{k}}_{\mathit{i}}}\right)}_{\mathit{w}}},{\left({\mathit{TS}}_{{\mathit{k}}_{\mathit{i}}},{\mathit{TS}}_{{\mathit{k}}_{\mathit{j}}}\right)}_{\mathit{w}}={\sum }_{\mathit{l}=1}^{\mathit{N}}{\mathit{w}}_{\mathit{l}}{\mathit{x}}_{\mathit{l}}^{{\mathit{k}}_{\mathit{i}}}{\mathit{x}}_{\mathit{l}}^{{\mathit{k}}_{\mathit{j}}}$$E27

with w_l = min (l, L, N-l) where L ≤ N/2.

An important issue is the selection of window length L which is a main parameter to determine the rank of the trajectory matrix. The window length is in the range [2, N/2]. The construction of trajectory matrix is similar to the phase space reconstruction of a nonlinear time series. The most common phase space reconstruction method is the method of delays (MOD). Many techniques have been suggested to estimate the parameters of MOD [22, 23, 24], i.e. the time delay τ and the embedding dimension m. The embedding dimension refers to the samples (separated by a fixed time delay τ) from a time window length L = (m-1)τ. The difference between the MOD and the trajectory matrix of SSA is that in MOD the m coordinates are samples separated by a fixed τ which can be greater than 1, while in the standard SSA, all the available samples in window length L are samples separated by a constant τ = 1. For the selection of window length L in SSA, there are some typical methods available, e.g. Kugiumtzis [25] suggested to take the window length L proportional to the mean orbital period” (MP), which can be determined from the spectral density of the original time series:

$$\mathit{MP}=\frac{{\sum }_{\mathit{i}=1}^{\mathit{M}}\mathit{DS}\left(\mathit{f}\right)}{{\sum }_{\mathit{i}=1}^{\mathit{M}}\mathit{f}.\mathit{DS}\left(\mathit{f}\right)}+1$$E28

where M is the number of frequencies f contained in the spectral support of the original time series, DS(f) is the spectral density of the f frequency. Several approaches to determine an optimal window length L are proposed in the literature [26, 27]. Unfortunately, those proposed methods are all case sensitive and based on the assumption that the measured time series is stationary. So far there does not exist any universal method to find an optimal window length L for an arbitrary time series.

3. Time series analysis and results

Numerical techniques mentioned above as EMD, EAWD, SSA have been applied to three monthly experimental time series of observation. The first time series is 22 years total ozone columns (in Dobson units) from January 1998 recorded in Natal (Brazil). The second time series is 42 years total ozone columns from 1978 recorded in Argentina. The third one is 60 years rainfall time series recorded in Conakry (Guinea).

3.1 Natal monthly total columns ozone time series analysis

The original time series is displayed in Figure 1, its size is N = 264:

Figure 1.

Natal monthly total ozone columns.

An EEMD technique has been applied to the time series presented above. After having selected relevant IMFs and add irrelevant ones to the residue, 5 IMFs and the trend has been determined. Results are shown in the Figure 2.

Figure 2.

EEMD decomposition of the Natal total ozone columns time series.

The IMF cycle is determined by the dominant frequency of its spectrum. In Figure 3, the residual trend R returned by the EEMD is compared with the trend of the original signal obtained from a moving window with a size set at the maximum of the relevant IMF cycles, i.e. 135 months named Tmb.

Figure 3.

Trend obtained by the EEMD technique (black curve). The trend of the original signal obtained from a moving window (blue curve) with a 135 months size.

The accuracy of the Trend returned by the EEMD is estimated by using the following expression:

$${\mathit{AC}}_{\mathit{\text{EEMD}}}=\sqrt{\frac{{\sum }_{\mathit{i}=1}^{\mathit{N}}{\left(\mathit{R}\left(\mathit{i}\right)-\mathit{Tmb}\left(\mathit{i}\right)\right)}^2/N}{\overline{\mathit{Tmb}}}}$$E29

where N is the size of the original time series and $\stackrel{-}{\dots }$ represents the mean operator.

No mode-mixing has been detected in the results returned by the EEMD. However, because the EEMD suffers from a lack of theory, the spectral supports of IMFs are not totally disjoint. The EAWD method allows, while relying on the spectral content of the IMFs provided by the EEMD to restitute components whose spectral contents are disjoint. In this context, the EAWD technique has been applied to the Natal time series and the obtained results are displayed on Figure 4.

Figure 4.

EAWD decomposition of the Natal total ozone columns time series.

The accuracy of the trend returned by the EAWD has been estimated by using (Eq. (29)). AC_EAWD = 0.2%.

The trend restituted by respectively EEMD and EAWD techniques have been superimposed with the trend obtained from the original signal by using a moving window of size 135 months as shown in Figure 5.

Figure 5.

Superimposition of the EAWD trend (red curve), EEMD trend (black curve) and the trend of the original signal obtained by using a moving window (blue curve).

The SSA technique has been applied to the Natal total ozone columns. The mean orbit period MP has been estimated by using (Eq. (28)), MP = 11. The embedding window size L has been fixed to L = 110. After having removing noisy components and applying the SSA basic procedure, 37 principal components have been found including trend and relevant periodic components. As the EEMD, EWT, EAWD acts as a bank of dyadic filter, to compare the results provided by SSA with those returned by EEMD, EAWD techniques, cycles of the 37 components extracted from the original signal by SSA are distributed along a grid $\left[2^{\frac{\mathit{log}\left(\mathit{\text{Cymax}}\right)}{\mathit{log}\left(2\right)}},2^2\right]$. Thus, the SSA components whose cycles are in the interval [2ⁱ,2^i + 1] are summed to form a component that can be compared to those returned by EEMD, EAWD methods. The components thus obtained are displayed in the Figure 6.

Figure 6.

SSA results on Natal ozone total columns time series.

Comparison between EAWD and SSA results are displayed on Figure 7.

Figure 7.

Comparison between SSA results (black curve) and EAWD results (red curve).

Picture (a), (c) and (d) represents the comparison between SSA and EAWD for respectively the cycles 135 months, 12 months and 6 months. In picture (c) it’s the superimposition of summed EAWD components 2 and 3 (cycles 33.75 months and 27 months) with SSA component 4 (cycle 27 months). Picture (c) is interesting from a climate forcing point of view. SSA2 component can be identified to the ENSO climate forcing (45 months) and on the graph of SSA2 (Figure 6), we can see that SSA2 energy contribution is very weak. SSA3 component can be identified to the QBO climate forcing (27 months). In Figure 6(c), the red curve is the sum the EAWD components 2 and 3 of 33.75 months and 27 months respectively (Figure 4). The difference the two curves in picture (c) is due to the presence of a very weak contribution of ENSO climate forcing in red curve (EAWD results).

3.2 Argentina monthly total columns ozone time series analysis

The original time series is displayed in Figure 8, its size is N = 507.

Figure 8.

Argentina monthly total ozone columns.

EEMD results are displayed in Figure 9.

Figure 9.

EEMD results of Argentina total ozone columns.

Mode-mixings are detected between IMFs 2 and 3 and between IMFs 5 and 6 as their dominant frequency are the same. The accuracy of the Trend returned by the EEMD is estimated to be 0.6%. An EAWD technique has been applied then to overcome mode-mixings encountered when EEMD has been applied. Results are displayed in Figure 10.

Figure 10.

EAWD results of Argentina total ozone columns.

Mode-mixings occurring in EEMD results have been removed. In EAWD results three new components have been appeared with cycles of respectively 28.5, 21.3 and 6 months. The IMF with cycle 2.4 months returned by EEMD disappeared. The trend restituted by respectively EEMD and EAWD techniques have been superimposed with the trend obtained from the original signal by using a moving window of size 170 months (Figure 11).

Figure 11.

Superimposition of the EAWD trend (red curve), EEMD trend (black curve) and the trend of the original signal obtained by using a moving window (blue curve).

The accuracy of the trend returned by the EAWD (red curve) has been estimated to be equal to 0.2%. The SSA technique has been applied to the Argentina total ozone columns. The mean orbit period MP has been estimated by using (Eq. (28)), MP = 10. The embedding window size L has been fixed to L = 110. After having removing noisy components and applying the SSA basic procedure 28 principal components have been found including trend and relevant periodic components. SSA principal components have been distributed along a grid $\left[2^{\frac{\mathit{log}\left(\mathit{\text{Cymax}}\right)}{\mathit{log}\left(2\right)}},2^2\right]$. The components obtained are displayed in Figure 12.

Figure 12.

SSA results on Argentina ozone total columns time series.

Figure 13 displays the superimposition of the EAWD results (red curve) and SSA results (black curve). In pictures (a), (b) (d), (e) and (f) SSA results in black and EAWD results in red are shown for respectively the cycles 6, 12, 64, 170 months and trends. In picture (c) it’s the superimposition of summed EAWD components 3 and 4 (cycles 28 months and 21 months) with SSA component 4 (cycle 28 months). Curves in picture (c) can be identified to QBO climate forcing. Curves in pictures (d) and (e) can be respectively identified to ENSO and solar cycle.

Figure 13.

Superimposition of SSA results (black curve) and EAWD results (red curve).

3.3 Conakry monthly rainfall time series analysis

The original time series is displayed in Figure 14, its size is N = 720:

Figure 14.

Conakry monthly rainfall in millimeters.

EEMD results are displayed in Figure 15.

Figure 15.

EEMD results of Conakry rainfall.

Mode-mixings are detected between IMFs 2 and 3 and between IMFs 5 and 6 as their cycles are the same (12 months). An EAWD technique has been applied then to overcome mode-mixings encountered when EEMD has been applied. Results are displayed in Figure 16.

Figure 16.

EAWD results of Conakry rainfall time series.

Mode-mixing occurring between IMFs 2 and 3 in EEMD results has been removed. The trend restituted by respectively EEMD and EAWD techniques have been superimposed with the trend obtained from the original signal by using a moving window of size 144 months (Figure 17).

Figure 17.

Superimposition of the EAWD trend (red curve), EEMD trend (black curve) and the trend of the original signal obtained by using a moving window (blue curve).

The accuracy of the trend returned by the EAWD (red curve) has been estimated to be equal to 0.97%. The SSA technique has been applied to the Conakry rainfall time series. The mean orbit period MP has been estimated by using (Eq. (28)), MP = 9. The embedding window size L has been fixed to L = 108. After having removing noisy components and applying the SSA basic procedure 42 principal components have been found including trend and relevant periodic components. SSA principal components have been distributed along a grid $\left[2^{\frac{\mathit{log}\left(\mathit{\text{Cymax}}\right)}{\mathit{log}\left(2\right)}},2^2\right]$. The components obtained are displayed in Figure 18.

Figure 18.

SSA results on Conakry rainfall time series.

Figure 19 displays the superimposition of the EAWD results (red curve) and SSA results (black curve). In pictures (a), (b), (c), (d) (e), and (f) SSA results in black and EAWD results in red are shown for respectively the cycles 6, 12, 32.7, 65.5, 144 months, and trend. As the segmentation of the original signal spectrum is not identical in EAWD and SSA, a 21 months cycle component appears in the EAWD decomposition and a 9 months cycle component appears in the SSA decomposition.

Figure 19.

Superimposition of SSA results (black curve) and EAWD results (red curve).

4. Conclusion

A comparative analysis of signal decomposition methods has been performed. The methods include empirical mode decomposition (EMD), ensemble empirical mode decomposition (EEMD) which can be seen as an optimization of EMD against noise, empirical wavelet transform and empirical adaptive wavelet decomposition which are methods that allows the mode-mixing problem to be avoid, and finally the singular spectrum decomposition which acts as a bank of narrow band filters and requires grouping techniques to be compared with results obtained with EAWD. Assuming that the variability components returned by the EMD are close to the physics implicitly contained in the observation time series. The EAWD technique seems to be an advantageous technique as it combines the rigor of wavelets while taking into account the spectral content of the components returned by the EMD. Thus the mode-mixing is avoid and in addition if an EEMD technique is used instead EMD, this technique becomes more robust against noise. Moreover, such a technique does not require a priori specification of the number of decomposed components. However, the SSA technique requires setting parameters like the width of the embedding window L and selecting an appropriate grouping technique to obtain components similar to those returned by the EAWD. SSA technique was found to be sensitive to parameter changes and L determines how well the components are separated. Despite the fact that SSA requires an appropriate embedding window size L, SSA is close to the underlying dynamics to the time series. Moreover, as the components extracted by SSA are independent, such a technique can be very useful for blind source separation purpose.

Acknowledgments

This study was led in the scope of the CNRS-LEFE-MANU SOLSTYCE project. Rainfall data were collected and qualified in the framework of the French South-African International Research Group IRP-ARSAIO (International Research Project - Atmospheric Research in Southern Africa and the Indian Ocean) supported by the NRF and CNRS, while the ozone time series were constructed within the French-Brazilian COFECUB-CAPES program.