3.1 Sifting

The algorithm will be outlined here as it was introduced in Huang et al. (Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu1998) before the various extensions are discussed. To this end, some notation is introduced here. Let $ |\{\cdot\}| $ be the cardinality of the set, $ \gamma_k(t) $ be the $ k{\rm th} $ IMF defined over $ t \in [0, T] $ , $ \tilde{\gamma}^M_k(t) $ be the maxima envelope fitted through the maxima, $ \tilde{\gamma}^m_k(t) $ be the minima envelope fitted through the minima, and let $ \tilde{\gamma}^{\mu}_k(t) $ be the mean envelope being the average of $ \tilde{\gamma}^M_k(t) $ and $ \tilde{\gamma}^m_k(t) $ , with $ \frac{d\gamma_k(t)}{dt} $ and $ \frac{d^2\gamma_k(t)}{dt^2} $ being the first and second derivatives, respectively, of the IMF, $ \gamma_k(t) $ . With this notation, the two defining conditions of an IMF (first introduced by Campi et al. Reference Campi, Peters, Azzaoui and Matsui2021) are expressed mathematically as below:

This definition, specifically Condition 2, while precise in its formulation of one of the characteristics of an IMF that are to be extracted by the EMD sifting procedure, is not practical for implementation. Satisfying this condition is difficult to attain in practice and may lead to an excessive computational effort which can also often lead to the propagation of errors and, potentially, the removal of meaningful structures from the signal basis representation. As such, one often needs to apply an approximation that is more computationally efficient and iteratively favourable. Such a definition can be seen below:

for some chosen $ \epsilon $ with $\tilde{\gamma}^{\mu}_k(t)$ defined discretely ( $t \in \{t_0, \dots, t_N\}$ ) during implementation. If Condition 1 and Modified Condition 2 are met, then the IMF is stored and removed from the remainder or residual ( $ r_k(t) = r_{k-1}(t) - \gamma_k(t) $ ), if either are not met then the process is repeated (or an additional optional third condition, known as a stopping criterion, is checked). The process is repeated until an IMF candidate either meets the required conditions or a stopping criterion is met. The end result of the sifting procedure will be an IMF representation of the original time series, x(t), with the total number of IMFs being denoted by K and the remainder or trend being denoted by $ r_K(t) $ :

(1) \begin{equation} x_{IMF}(t) = \sum_{k=1}^{K}\gamma_k(t) + r_K(t).\end{equation}

3.2 Hilbert transform and instantaneous frequency

Once the basis of the time series has been extracted, they need to be analysed using the HT to provide a time-frequency characterisation. The IMF basis represents the time-domain signal decomposition, and the objective is to transform this basis decomposition into a frequency-based characterisation. The joint process of decomposing the time series using EMD and Hilbert transforming the result is known as the HHT. The HT of a time series only makes physical sense if the time series satisfies the IMF conditions. Given an IMF, $ \gamma_k(t) $ , the corresponding HT, $ \check{\gamma}_k(t) = \mathcal{HT}[\gamma_k(t)] $ with $ \mathcal{HT}[\cdot] $ being the HT, is calculated as follows:

(2) \begin{equation}\check{\gamma}_k(t) = \frac{1}{\pi}PV\int_{-\infty}^{\infty}\frac{\gamma_k(t^*)}{t - t^*}dt^*,\end{equation}

with PV being the Cauchy Principle Value integral. It is an improper integral that assigns values to otherwise undefined integrals. The Cauchy principal value integral is formally defined as:

(3) \begin{equation} \check{\gamma}_k(t) = \frac{1}{\pi}\lim_{\epsilon\rightarrow{0}^+}\bigg[\int_{t - \frac{1}{\epsilon}}^{t - \epsilon}\frac{\gamma_k(t^*)}{t - t^*}dt^* + \int_{t^* + \epsilon}^{t + \frac{1}{\epsilon}}\frac{\gamma_k(t^*)}{t - t^*}dt^*\bigg].\end{equation}

The analytical signal, $ \gamma^a_k(t) $ , can then be defined as:

(4) \begin{equation}\gamma^a_k(t) = \gamma_k(t) + i\check{\gamma}_k(t) = a_k(t)e^{i\theta_k(t)} = a_k(t)e^{i\int\omega_k(t)dt},\end{equation}

with $ a_k(t) $ and $ \theta_k(t) $ being defined as:

(5) \begin{equation}a_k(t) = \sqrt{\gamma_k(t)^2 + \check{\gamma}_k(t)^2},\end{equation}

and

(6) \begin{equation}\theta_k(t) = \mbox{tan}^{-1}\left(\frac{\check{\gamma}_k(t)}{\gamma_k(t)}\right).\end{equation}

The IF, $ \omega_k(t) $ , for the k-th IMF, is then defined as:

(7) \begin{equation} \omega_k(t) = \frac{d\theta_k(t)}{dt}.\end{equation}

The Hilbert spectrum of IMF $\gamma_k(t)$ is calculated as:

(8) \begin{equation}H_k(t)\;:\!=\; H(t;\;\omega_k) =\left\{\begin{array}{l@{\quad}l}a_k(t), & \omega_{k} \leq \omega(t) \lt \omega_{k+1}\\[5pt] 0, & {}\mbox{otherwise}\end{array},\right.\end{equation}

with $\omega_{k}$ and $\omega_{k+1}$ being the discretisation of the frequency domain for analysis and plotting purposes – higher resolutions require higher computational expense and without any smoothing (Gaussian smoothing recommended – see Figures 13 and 14 for context) would be imperceptible to most in plots. The Hilbert spectrum allows for a far higher resolution than other frequency analysis techniques such as the short-time Fourier transform and Morlet wavelet transform.

3.3 Cubic B-spline empirical mode decomposition

One will immediately notice that the conditions on an IMF are characterising conditions but they are not constructive; in other words, one must still define a functional form for parameterisation of the representation of the IMFs. This is distinct from methods such as Fourier decomposition where a fixed basis, in this case, cosines, is determined as fundamentally part of the basis decomposition. This is both advantageous and challenging for EMD methods. On the one hand, it provides a lot of flexibility to the modeller to achieve the conditions for the IMF in the EMD sifting procedure; on the other hand, it opens up a model specification and selection problem. Fortunately, to some extent, this question has been addressed in the literature, at least to the extent that numerous choices have been proposed.

This section follows the work done in Chen et al. (Reference Chen, Huang, Riemenschneider and Xu2006) on applying cubic B-splines to EMD. In this case, one expresses the EMD representation as follows:

(9) \begin{equation} x_{IMF}(t) = \sum_{k=1}^{K}\gamma_k(t) + r_K(t),\end{equation}

with each $ \gamma_k(t) $ and the trend, $ r_K(t) $ , being the sum of compact cubic B-splines. The optimality of cubic splines concerning MSE minimisation has been shown in Craven & Wahba (Reference Craven and Wahba1978), Wahba (Reference Wahba1990), and Bowman & Evers (Reference Bowman and Evers2017), and in addition, the cubic B-splines have the further practical advantage that the HT of B-splines can be obtained in closed-form in a recursive representation. B-splines and the recursive relationship that exists to define the higher-order splines are shown in de Boor (Reference de Boor1978). The B-spline of order 1 is shown below:

(10) \begin{equation} B_{i,1,\boldsymbol{\tau}}(t) = \mathbb{1}_{[\tau_{i},\tau_{i+1})}.\end{equation}

This is the base case and as such the higher-order splines can be defined recursively. $ B_{i,j,\boldsymbol{\tau}}(t) $ is the $ i{\rm th} $ basis function of order j defined over knot sequence $ \boldsymbol{\tau} $ with $ t \in \mathbb{R} $ . With this framework, $ B_{i,j,\boldsymbol{\tau}}(t) $ is defined recursively as:

(11) \begin{equation} B_{i,j,\boldsymbol{\tau}}(t) = \dfrac{t - \tau_i}{\tau_{i+j-1} - \tau_i} B_{i,j-1,\boldsymbol{\tau}}(t) + \dfrac{\tau_{i+j} - t}{\tau_{i+j} - \tau_{i+1}}B_{i+1,j-1,\boldsymbol{\tau}}(t).\end{equation}

Equation (11) is proven in de Boor (Reference de Boor1978). The HT of Equation (10), with $ \check{B}_{i,1,\boldsymbol{\tau}}(t) = \mathcal{HT}[B_{i,1,\boldsymbol{\tau}}(t)], $ can be shown to be:

(12) \begin{equation} \check{B}_{i,1,\boldsymbol{\tau}}(t) = \frac{1}{\pi}ln\bigg|\dfrac{t-\tau_i}{t-\tau_{i+1}}\bigg|,\end{equation}

with singularities at $ \tau_i $ and $ \tau_{i+1} $ . It is for this reason that the knot sequence is advised to be a proper subset of the time series. The corresponding HT recursive relationship is proven in Chen et al. (Reference Chen, Huang, Riemenschneider and Xu2006) and is stated below:

(13) \begin{equation} \check{B}_{i,j,\boldsymbol{\tau}}(t) = \dfrac{t - \tau_i}{\tau_{i+j-1} - \tau_i}\check{B}_{i,j-1,\boldsymbol{\tau}}(t) + \dfrac{\tau_{i+j} - t}{\tau_{i+j} - \tau_{i+1}}\check{B}_{i+1,j-1,\boldsymbol{\tau}}(t).\end{equation}

With this framework in place, the cubic B-spline fitting problem can be stated. The time series is defined at discrete points, $ \textbf{t} = \{t_0, \dots, t_N\} $ , so that the time series can be represented by a vector, $ \boldsymbol{{s}} $ , below:

(14) \begin{equation}\boldsymbol{{s}} = \begin{bmatrix}s(t_0)\\[5pt] s(t_1)\\[5pt] \vdots\\[5pt] s(t_{N})\end{bmatrix}.\end{equation}

The knot sequence is implied by expressing the problem objective function in matrix form, and as such, the knot sequence subscript is dropped below:

(15) \begin{equation}\boldsymbol{{B}} = \left[\begin{array}{ccc}B_{0,4}(t_0) {}& \cdots & {} B_{(M-4),4}(t_0)\\[5pt] B_{0,4}(t_1) {}& \cdots & {} B_{(M-4),4}(t_1)\\[5pt] \vdots {}& \ddots {}& \vdots\\[5pt] B_{0,4}(t_{N}) {}& \cdots & {} B_{(M-4),4}(t_{N}) \end{array}\right].\end{equation}

Finally the coefficient vector, $ \boldsymbol{{c}} $ , is defined below:

(16) \begin{equation}\boldsymbol{{c}} = \begin{bmatrix}c_0\\[5pt] c_1\\[5pt] \vdots\\[5pt] c_{(M-4)}\end{bmatrix}.\end{equation}

Using the above, the objective function can be defined as follows:

(17) \begin{equation} MSE(\boldsymbol{{c}}|\boldsymbol{{s}}) = (\boldsymbol{{s}} - \boldsymbol{{B}}\boldsymbol{{c}})^T(\boldsymbol{{s}} - \boldsymbol{{B}}\boldsymbol{{c}}).\end{equation}

The result of this optimisation is a cubic B-spline interpolation spline, p(t), such that:

(18) \begin{equation} p(t) = \sum_{i=0}^{M-4}c_iB_{i,4,\boldsymbol{\tau}}(t).\end{equation}

The power of using cubic B-splines is that the HT of this spline, $ \check{p}(t) = \mathcal{HT}[p(t)] $ , is simply:

(19) \begin{equation} \check{p}(t) = \sum_{i=0}^{M-4}c_i\check{B}_{i,4,\boldsymbol{\tau}}(t).\end{equation}

Several variations of this general framework will be explored. An explicit smoothing term can be incorporated in Equation (17) to create explicitly smoothed cubic B-splines as in Equation (36). These discretely penalised B-splines are referred to as P-splines. The number and distribution of knots can be changed to implicitly smooth the spline either manually or using knot point optimisation as in section 7.8, and the knot points are extended beyond the edge of the time series to create a non-natural spline to better fit the edges of the time series as shown in Figure 2.

Assuming a consistent set of knots throughout the sifting procedure (either optimise_knots=0 or optimise_knots=1, but more about this in sections 7.1 and 7.8) and by letting remaining trend structure K, $ r_K(t) $ , be denoted as structure $ K+1 $ , Equation (9) can be expressed as:

(20) \begin{align} x_{IMF}(t) {}&= \sum_{k=1}^{K}\sum_{i=0}^{M-4}c_{i,k}B_{i,4,\boldsymbol{\tau}}(t) + r_K(t)\nonumber \\[5pt] {}&= \sum_{k=1}^{K+1}\sum_{i=0}^{M-4}c_{i,k}B_{i,4,\boldsymbol{\tau}}(t),\end{align}

with $ c_{i,k} $ denoting coefficient i for structure k.

4.1 Python – emd 0.3.3

In this package, Quinn (Reference Quinn2020), there are numerous variations of EMD available. The standard univariate EMD does not have many algorithmic variations. The stopping criteria used are the Cauchy-type convergence (section 10.5.4 of the supplement), mean fluctuation threshold (section 10.5.7 of the supplement), and fixed iteration (section 10.5.2 of the supplement). The splines used in the packages are cubic B-splines, piece-wise cubic Hermite splines, and monotonic piece-wise cubic Hermite splines. The edges are reflected to prevent extreme extrapolation at the edges causing errors to propagate. Ensemble EMD (EEMD) (Wu & Huang, Reference Wu and Huang2009) (section 7.9), complete EEMD (CEEMD) (Torres et al., Reference Torres, Colominas, Schlotthauer and Flandrin2011), masking signal-assisted EMD (Deering & Kaiser, Reference Deering and Kaiser2005), and Holo-Hilbert spectral analysis (Huang et al., Reference Huang, Hu, Yang, Chang, Jia, Liang, Yeh, Kao, Juan, Peng, Meijer, Wang, Long and Wu2016) are also available. CEEMD adds noise to every IMF candidate, rather than once at the start of each iteration of the algorithm – this does show promising results but deserves further study with an appropriate sorting mechanism. Adding a masking signal is very time series specific – it depends on the frequency and amplitude composition of the IMF candidates and therefore usually requires an initial unmasked sifting. Holo-Hilbert spectral analysis attempts to handle more complex time series combinations.

4.2 Python – PyEMD (EMD-signal) 0.2.10

In this package, Laszuk (Reference Laszuk2020), some variations of EMD are available. Univariate EMD does have limited algorithmic variations. The splines available are cubic, Akima (section 10.6.2 of the supplement), and linear. The stopping criteria used are the standard Cauchy-type convergence (section 10.5.4 of the supplement), fixed iteration (section 10.5.2 of the supplement), and S Stoppage (section 10.5.3 of the supplement). Extrema are calculated using standard finite difference methods or by fitting a parabola and inferring extrema from the parabola. Extrema at the edges are reflected to deal with the edge effect. Different algorithm variations are also available: ensemble EMD (EEMD) (Wu & Huang, Reference Wu and Huang2009) (section 7.9), complete EEMD (CEEMD) (Torres et al., Reference Torres, Colominas, Schlotthauer and Flandrin2011), and two-dimensional EMD grouped into Bi-dimensional EMD (BEMD) and two-dimensional EMD (EMD2D) – both have been classified as in-progress or abandoned projects.

4.3 R – EMD 1.5.8

This package, Kim & Oh (Reference Kim and Oh2018), is a more complete and refined EMD library. The stopping criteria used are fixed iteration (section 10.5.2 of the supplement), S Stoppage (section 10.5.3 of the supplement), modified mean (section 10.5.1 of the supplement), Cauchy-type converge (section 10.5.4 of the supplement), Cauchy-type 11a (section 10.5.5 of the supplement), and Cauchy-type 11b (section 10.5.6 of the supplement). The optional edge effects are one of the symmetric edge effects discussed in section 7.3.1, the anti-symmetric edge effect (section 7.3.1), the Huang characteristic wave (section 10.3.3 of the supplement), and the Coughlin characteristic wave (section 10.3.3 of the supplement). Cubic spline interpolation is done. Furthermore, several smoothing techniques are used such as spline smoothing, kernel smoothing, and local polynomial smoothing. As a result, various forms of statistical EMD (SEMD) are available. Plotting of individual IMF candidates and final IMFs is available for debugging purposes, and certain periods may be excluded in a different method to prevent mode mixing. Confidence limits for the IMFs may also be calculated using vector auto-regressive techniques. Bi-dimensional EMD is available as well as a de-noising version of EMD. A limit can be imposed on the maximum number of IMFs.

4.4 MATLAB – Empirical Mode Decomposition 1.0.0.0

The package, Ortigueira (Reference Ortigueira2021), is older and more tested with few, if any, recent additions. Cauchy-type convergence (section 10.5.4 of the supplement) is the stopping criterion used in this package. Fixed iterations (section 10.5.2 of the supplement) is another criterion used. An energy ratio stopping criterion is used that may be directly translated to a rephrasing of a relative modified mean criterion (section 10.5.1 of the supplement). A limit is imposed on the allowed number of IMFs. Spline methods used are cubic splines and piece-wise cubic Hermite splines. This package, while having few algorithmic variations or extensions, is well-tested – much like Kim & Oh (Reference Kim and Oh2018).

7.1 Base implementation of AdvEMDpy package

Before discussing the nuances of the algorithmic variations, one would benefit from some choice remarks concerning non-essential, but helpful outputs of the AdvEMDpy package. These play no direct part in the implementation of the algorithm but assist in iterative understanding and assessment of EMD’s performance. The base implementation of the algorithm may be seen below. More detailed descriptions may be found in the associated scripts as well as a base implementation in the associated README.md file.

# instantiate EMD class with time (optional keyword argument) and time series

emd = EMD(time=time, time_series=time_series)

# base implementation with optional helpful inputs shown

# knots and knot_time are optional keyword arguments

imfs, hts, ifs = \

emd.empirical_mode_decomposition(knots=knots,

knot_time=knot_time,

debug=False,

dft=‘envelopes’,

verbose=True,

stopping_criterion=‘sd’,

output_coefficients=False,

output_knots=False)[:3]

The details of each field in the input to the function are detailed in the following subsections. If time is not specified as a keyword argument, then the following code is implemented where time is simply an index set of the same length as the original time series.

self.time = np.arange(len(self.time_series))

If the knot sequence or the associated knot time are also unspecified keyword arguments, then the following code is implemented that sets knot time equal to time and spaces the knots 10 times further apart than the time points. This is problematic if the default knot sequence is insufficient to capture the highest frequency structures which is demonstrated in new_user.ipynb.

knots = np.linspace(0, self.time[-1], int(len(self.time)/10 + 1))

knot_time = self.time

The only required argument is the time series. This is detailed further in new_user.ipynb which is intended for new users to get comfortable with the package and the available options.

7.1.3 Output coefficients flag in AdvEMDpy package

If output_coefficients=True, cubic B-spline coefficients corresponding to each IMF output (including initial smoothed time series and trend) are output. One must remember the six additional coefficients corresponding to the three additional knots on either edge of the time series as demonstrated in Figure 2. These additional bases result in implicit smoothing through a non-natural spline which does not tend to zero at the edges in both position and the higher-order derivatives. It would be unnecessarily restrictive to impose a natural spline on the sifting algorithm and the resulting IMFs.

Figure 2 Figure demonstrating incomplete bases at the boundary of the time series to create non-natural cubic B-splines that can accommodate non-zero edge values, derivatives, and curvatures.

7.1.4 Output knots flag in AdvEMDpy package

If output_knots=True, knot points are output. If knot points are not optimised, that is if optimise_knots=0, then the original knot points are repeated at each iteration of the sifting procedure when extracting recursively each IMF. In other words, all IMFs will have the same knot sequence. If the user provides a uniform knot sequence and chooses not to optimise the knot locations, then a uniform knot sequence will be used throughout as demonstrated in Figure 3.

Figure 3 Figure demonstrating predefined uniform knot placement and the resulting IMFs.

Alternatively, one can optimise the knots on an original spline (optimise_knots=1) representation of the input time series signal using optimal knot placement methodology as outlined in section 7.8 where two options are available using simple bisection and serial bisection, but more about this will follow in the dedicated section. In this way, the knot points will still be universal for all IMFs; however, they will be optimised for the given signal input and will not dynamically adjust as the residuals become increasingly less complex. For this reason, this can be referred to as the statically optimised knots, and an example demonstrating this can be seen in Figure 4 where the optimised non-uniform knots can be seen to be used throughout the algorithm despite the decreasing complexity.

Figure 4 Figure demonstrating statically optimised knot placement, which is optimised once at outset and used throughout the sifting, and the resulting IMFs.

In addition, one could also optimise the knot points per iteration (optimise_knots=2) of the EMD sifting procedure; once an IMF is isolated, the next iteration of EMD on the residual signal could first have the knot point optimisation performed when fitting the spline to the residual to proceed with the next rounds of sifting to extract the next IMF. This results in different numbers of knot points per IMF and different placements; in general, one would expect to require far fewer knots for later IMFs extracted which have lower frequency content compared to those first few IMF bases extracted. An example demonstrating this can be seen in Figure 5 where there are far fewer non-uniform knots for the second IMF than was required for the first IMF.

Figure 5 Figure demonstrating dynamically optimised knot placement, which is optimised at the beginning of each internal sifting routine, and the resulting IMFs.

Performing the third option may increase the speed of the algorithm if a large number of siftings need to take place and a large number of IMFs need to be extracted, despite the added time required to find the optimised knots at each stage. Progressively fewer knots will be required throughout the algorithm with the added advantage of having a more parsimonious representation with fewer parameters required for higher-order IMFs which have guaranteed reducing oscillation in their representation.

7.2 Pre-processing flag in AdvEMDpy package

The raw time series signals that can be passed to the EMD package may contain a wide variety of structures. This could include different non-linear structures, non-stationarity features in trend, volatility, etc., as well as corruption by observation noise that may be present in the collection of the data. Therefore, in the AdvEMDpy package there is the option to pre-process the time series data input to reduce the effect of these features if they may be adverse to the user’s analysis. This is particularly the case in the context of signals observed in various noisy environments.

Therefore, before the sifting algorithm is implemented, there are several choices available for the pre-processing of the raw time series. This was a necessary inclusion in the algorithm owing to the wide variety of time series and their potentially vastly different statistical characteristics that would need to be accommodated in the EMD basis decomposition. These pre-processing techniques may be broadly grouped into filtering and smoothing. The filtering methods developed out of a necessity to mitigate the corruption of the IMFs that would otherwise arise as a result of the permeation of error that would occur due to the presence in the input time series signal of noise corruptions such as heavy-tailed noise, mixed noise (Gaussian noise with different variances), and Poisson noise. Smoothing developed out of the broader field of trend extraction among cyclical components where within EMD defining an individual cyclical component among others and within a non-stationary setting becomes increasingly challenging. The base implementation is as follows:

imfs, hts, ifs = \

emd.empirical_mode_decomposition(knots=knots, knot_time=knot_time,

initial_smoothing=True,

preprocess=‘none’,

preprocess_window_length=51,

preprocess_quantile=0.9,

preprocess_penalty=1,

preprocess_order=13,

preprocess_norm_1=2,

preprocess_norm_2=1,

downsample_window=‘hamming’,

downsample_decimation_factor=20,

downsample_window_factor=20)[:3].

The details of the implementation of filtering techniques are available such as a mean filter, median filter, Winsorization, and Winsorization and interpolation; one can see section 10.2.2 in “Supplement to: Package AdvEMDpy: Algorithmic Variations of Empirical Mode Decomposition in $\textsf{Python}$ ”. The smoothing techniques available such as generalised Hodrick-Prescott smoothing, Henderson-Whittaker graduation, downsampling, and downsampling followed by decimation can be found in section 10.2.3 of the supplementary materials.

7.3 Edge effects in AdvEMDpy package

The most prevalent challenge in EMD is the propagation of errors throughout the IMFs and the permeation of the errors throughout the sifting process. Owing to the iterative nature of the algorithm, incorrectly estimated edges lead to errors being ubiquitous throughout all the IMFs. There are several techniques used to estimate the extrema beyond the edges of the signal – this AdvEMDpy package attempts to give the reader as many reasonable and researched options as possible. The base implementation of the edge effects code is as follows:

imfs, hts, ifs = \

emd.empirical_mode_decomposition(knots=knots,

knot_time=knot_time,

edge_effect=‘symmetric’,

sym_alpha=0.1,

nn_m=200,

nn_k=100,

nn_method=‘grad_descent’,

nn_learning_rate=0.01,

nn_iter=100)[:3].

7.3.1 Symmetric methods

Without loss of generality, the procedure will be explained for the right edge maxima of the signal. Figure 6 demonstrates examples of the three possible symmetric edge effects. In no particular order, the technique demonstrated in Rilling et al. (Reference Rilling, Flandrin and Goncalves2003) and Wu & Qu (Reference Wu and Qu2008) is referred to in this paper as the symmetric discard technique owing to the discarding of the end of the time series in the new extrema approximation.

Figure 6 Example time series demonstrating four (five if conditional symmetric anchor is included) different symmetric edge effect techniques with axes of symmetry included.

Symmetric Discard – In Figure 6, an implementation of the Symmetric Discard technique (edge_effect=“symmetric_discard”) maximum is shown with a purple dot. The end of the time series between the last blue dot and the orange dot is disregarded when approximating the next extreme. Taking the last maximum value, $ X(t^{\max}_M) $ , the associated time point, $ t^{\max}_M $ , the last minimum value, $ X(t^{\min}_m) $ , and the associated time point, $ t^{\min}_m $ , the next extreme is calculated as:

(21) \begin{equation} X(t^{\max}_{M+1}) = X(t^{\max}_M),\end{equation}

with the associated time point being calculated as:

(22) \begin{equation} t^{\max}_{M+1} = t^{\min}_{m} + (t^{\min}_{m} - t^{\max}_{M}).\end{equation}

Symmetric Anchor – In Figure 6, an implementation of the Symmetric Anchor technique (edge_effect=“symmetric_anchor”) maximum is shown with an orange dot. In Zhao & Huang (Reference Zhao and Huang2001) and Zeng & He (Reference Zeng and He2004), the technique creates an extreme at the endpoint – this is why this technique is referred to in this paper as the Symmetric Anchor technique. In this paper, this technique has been generalised to conditionally create an extreme depending on the difference in vertical displacement between the last two extrema and the difference in vertical displacement between the last extremum and the end of the signal – this can be referred to as the conditional symmetric anchor technique. The conditional symmetric anchor is calculated as follows – if $ \beta{L} \geq (1 - \alpha)L $ where $ \beta $ is the ratio of the vertical height between the last extremum and the end of the time series to the vertical height between the last two extrema, L is the vertical height between the last two extrema, and $ \alpha $ (sym_alpha=0.1) is the significance level input, then:

(23) \begin{equation} X(t^{\max}_{M+1}) = X(t_N),\end{equation}

with the associated time point being calculated as:

(24) \begin{equation} t^{\max}_{M+1} = t_N.\end{equation}

The above is the method followed (without conditions) in Zhao & Huang (Reference Zhao and Huang2001) and Zeng & He (Reference Zeng and He2004). If, however, $ \beta{L} \lt (1 - \alpha)L $ , then:

(25) \begin{equation}X(t^{\max}_{M+1}) = X(t^{\max}_M),\end{equation}

with the associated time point being calculated as:

(26) \begin{equation}t^{\max}_{M+1} = t_N + (t_N - t^{\max}_{M}).\end{equation}

The conditional symmetric anchor technique collapses to the symmetric anchor technique when $ \alpha = 1 $ . The other extreme where $ \alpha = -\infty $ leads to the following method.

Symmetric – The symmetric technique (edge_effect=“symmetric”) does not anchor the extrema envelope to the ends of the signal under any condition and is equivalent to the conditional symmetric anchor technique where $ \alpha = -\infty $ . The values are calculated as follows:

(27) \begin{equation}X(t^{\max}_{M+1}) = X(t^{\max}_M),\end{equation}

with the associated time point being calculated as:

(28) \begin{equation}t^{\max}_{M+1} = t_N + (t_N - t^{\max}_{M}).\end{equation}

This point is denoted with a grey dot in Figure 6.

Anti-Symmetric – The anti-symmetric (edge_effect=“anti-symmetric”) approach reflects the univariate signal about both axes – the approximated maximum will be the reflected minimum. It is formally calculated as:

(29) \begin{equation}X(t^{\max}_{M+1}) = X(t_N) + (X(t_N) - X(t^{\min}_m)),\end{equation}

with the associated time point being calculated as:

(30) \begin{equation}t^{\max}_{M+1} = t_N + (t_N - t^{\min}_{m}).\end{equation}

This point is denoted with a green dot in Figure 6. This technique is a more practical variation of that proposed in Zeng & He (Reference Zeng and He2004), where the points are reflected about the axis rather than about the endpoints.

Descriptive figures, like Figure 6, detailing the other techniques can be found in “Supplement to: Package AdvEMDpy: Algorithmic Variations of Empirical Mode Decomposition in $\textsf{Python}$ ”. The slope-based family of techniques approach the problem using finite differences by approximating the forthcoming extrema by estimating the forthcoming slopes – the details of which can be found in section 10.3.2 of supplement. Characteristic waves are another family of techniques that approximate the edges using sinusoids with the details in section 10.3.3 of the supplement. Explicit techniques, like the single-neuron neural network (without an activation function) which can be found in section 10.3.4 of the supplement, formally forecast the time series and forecast extrema are taken from the forecast time series.

7.4 Detrended fluctuation analysis in AdvEMDpy package

At EMD’s core, the algorithm iterates by progressively removing the local mean from the IMF candidate. Owing to the numerous ways of defining and estimating a local mean, local mean estimation should be referred to as detrended fluctuation analysis. Detrended fluctuation analysis was originally introduced in Peng et al. (Reference Peng, Buldyrev, Havlin, Simons, Stanley and Goldberger1994) and intended for a different purpose where data were partitioned into equal sets and detrended using discontinuous linear approximations of the local trends. The average of the variances of the detrended data in these partitioned sets was calculated before the subsets were increased in size and the process was repeated. This was done to determine the presence of long-memory processes or rather long-term correlation structures. For our purposes, each set of data (data points between two consecutive knots) is detrended using continuous cubic B-splines to approximate the local trend or local mean. The trends are extracted using local windows defined using the knot sequences. The time series can therefore be estimated as a series of bases and thought of as a sequence of locally defined segments.

In Figure 7, detrended fluctuation analysis is demonstrated using the following time series:

(31) \begin{equation} f(t) = \text{cos}(2t) + \text{cos}(4t) + \text{cos}(8t) + \epsilon(t),\end{equation}

where $ \epsilon(t) \in \mathcal{N}(0,1) $ for $ t \in \{t_0, \dots, t_{1,000}\} $ . One can observe in Figure 7 that as the distance between the uniformly placed knots is increased, the fitted splines cannot detect the higher-frequency structures.

Figure 7 Example time series demonstrating detrended fluctuation analysis of time series with different knot sequences resulting in different trend estimation.

In Figure 7, one can observe the relationship between detrended fluctuation analysis and the iterative IMF extraction method known as EMD. In the top image of Figure 7, one can observe the trend as the sum of all the IMFs extracted with a sufficient knot sequence. In the middle image, the trend is approximated with a knot sequence that is insufficient to capture the highest frequency structure. One can also see how closely this trend resembles the sum of IMF 2 and IMF 3 for EMD when the knot sequence is sufficient. Finally, in the bottom image, one can see how the trend is estimated with a wholly insufficient knot sequence compared with the lowest order IMF of each previous sifting procedure.

In Figure 8, the different frequency structures are more easily visible. It can be shown that:

(32) \begin{equation} \text{lim}_{n\rightarrow{\infty}}\sqrt{\frac{1}{n-1}\sum_{k=1}^{n}\text{cos}^2\bigg(\frac{2\pi{fk}}{n}\bigg)} = \frac{1}{\sqrt{2}} \approx 0.707,\end{equation}

with $ f \in \mathbb{N} $ and $ f <\lt n $ . With this in mind, and noting the structure of Equation (31), one can understand the following results:

(33) \begin{align} SD\bigg(f(t) - \sum_{i=1}^3{IMF^{51}_i(t)}\bigg) {}&= 1.005 \approx 1 \nonumber \\[5pt] SD\bigg(f(t) - \sum_{i=1}^2{IMF^{31}_i(t)}\bigg) {}&= 1.558 \\[5pt] SD\bigg(f(t) - IMF^{11}_1(t)\bigg) {}&= 2.072, \nonumber\end{align}

with $ IMF^{j}_i $ being the $ i{\rm th} $ IMF with knot sequence j. The first equality in Equation (33) most closely calculates the true underlying noise present in the system as a result of the random fluctuations caused by the standard normally distributed Gaussian noise. In the other equalities, the standard deviation calculations are confounded by undetected underlying high-frequency structures with their individual “standard deviations” being approximated by Equation (32). The parallels between detrended fluctuation analysis and EMD should be noted by the user. By noting Figures 7, 8, and 9, one can observe that as a result of the defined knot sequence that fluctuation is measured relative to some frequency band. Each IMF exists in some frequency range, and the fluctuations calculated in Equation (33) measure the local fluctuation relative to some implicit frequency boundary as a result of the knot sequence defined.

Figure 8 Example time series demonstrating detrended fluctuation analysis of time series with different knot sequences resulting in different trend estimation.

Figure 9 Hilbert spectrum of example time series demonstrating the frequencies of the three IMFs present when sufficient knots are used.

In Figure 9, the frequencies of the three constituent structures (apart from the added standard normal Gaussian noise) are visible. The highest frequency structure is perturbed by the noise as expected. EMD can be viewed as a generalisation of detrended fluctuation analysis whereby the trend is decomposed into separate frequency structures in descending order of IF using a defined knot sequence (in general).

The local mean can be estimated using the framework in section 3.3 without any explicit or implicit smoothing. Explicit smoothing that deals simultaneously with smoothing and the edge effects is referred to in the literature as Statistical EMD (SEMD). The base implementation of the detrended fluctuation analysis is as follows:

imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,

knot_time=knot_time,

smooth=True,

smoothing_penalty=0.1,

dft=‘envelopes’,

order=15,

increment=10)[:3].

7.4.1 Statistical EMD (P-splines)

This method (smooth=True) is introduced in Kim et al. (Reference Kim, Kim and Oh2012). SEMD, in the cubic B-spline envelope setting, can be implemented by introducing a smoothing parameter into the objective function, Equation (17). The specific introduction of a penalty into B-splines is discussed in Eilers & Marx (Reference Eilers and Marx1996) and is referred to as P-splines. The matricification of the P-spline objective function can be seen below. Second-order smoothing is done on the coefficients, but this can be generalised to higher-order smoothing. The second-order smoothing of the coefficients is incorporated using the $ \boldsymbol{{D}} $ matrix seen below:

(34) \begin{equation}\boldsymbol{{D}} = \begin{bmatrix}1\;\;\; & {} -2\;\;\; & {} 1\;\;\;& {} 0\;\;\;& {} 0\;\;\;& {} \cdots\;\;\; & {} 0\\[5pt] 0\;\;\; & {} 1\;\;\; & {} -2\;\;\; & {} 1\;\;\; & {} 0\;\;\; & {} \cdots\;\;\; {}& 0\\[5pt] \vdots\;\;\; {} & & {} \ddots\;\;\; {}& \ddots\;\;\; & {} \ddots\;\;\; {}& & {} \vdots \\[5pt] 0\;\;\; & {} \cdots\;\;\; & {} 0\;\;\; & {} 1\;\;\; {}& -2\;\;\; {}& 1\;\;\; & {} 0\\[5pt] 0\;\;\; {}& \cdots\;\;\; {}& 0\;\;\; {}& 0\;\;\; & {} 1\;\;\; {}& -2\;\;\; & {} 1 \end{bmatrix},\end{equation}

and with $ \boldsymbol{{P}} $ defined as below,

(35) \begin{equation} \boldsymbol{{P}} = \boldsymbol{{D}}^T\boldsymbol{{D}},\end{equation}

Equation (17), with $ \boldsymbol{{s}} $ , $ \boldsymbol{{B}} $ , and $ \boldsymbol{{c}} $ defined as before and with discrete penalty term $ \lambda $ , becomes:

(36) \begin{equation} DPMSE(\boldsymbol{{c}}|\boldsymbol{{s}}) = (\boldsymbol{{s}} - \boldsymbol{{B}}\boldsymbol{{c}})^T(\boldsymbol{{s}} - \boldsymbol{{B}}\boldsymbol{{c}}) + \lambda\boldsymbol{{c}}^T\boldsymbol{{P}}\boldsymbol{{c}}.\end{equation}

The magnitude of the penalty term, $ \lambda $ (smoothing_penalty=0.1), determines the amount of smoothing. SEMD is implemented separately to the other detrended fluctuation techniques as smoothing can be done irrespective of the detrended fluctuation technique used.

It is highly recommended that smooth=True and the smoothing penalty is non-zero when dft=“envelopes” as the extrema are not guaranteed to satisfy the Schoenberg–Whitney Conditions (SWC) for an arbitrary knot sequence. This will allow envelopes to be fitted even when there are no extrema between successive knot points; otherwise, nonsensical envelopes may result. For an arbitrary set of extrema $ \textbf{y} = \{y_1, y_2, \dots, y_j\} $ and a cubic B-spline knot sequence $ \boldsymbol{\tau} = \{\tau_1, \tau_2, \dots, \tau_{j+4}\} $ , the SWC can be stated as:

(37) \begin{equation} \tau_i \leq y_i, \leq \tau_{i+1}\text{ }\forall\text{ }i \in \{1, 2, \dots, j-1\}.\end{equation}

In Figure 10, the following time series is plotted to demonstrate the necessity for either the SWC to be satisfied or for the envelopes to be smoothed:

(38) \begin{equation} g(t) = \text{cos}(t) + \text{cos}(5t),\end{equation}

with $ t \in [0, 5\pi] $ . With $ \lambda = 0 $ in Equation (36) and with the SWC not being satisfied by either the maxima or the minima, the envelopes are stretched towards zero. This is, unfortunately, unavoidable without either the SWC being satisfied or some form of smoothing.

Figure 10 Example time series demonstrating unsmoothed extrema envelopes being fitted when SWC are not satisfied resulting in nonsensical envelopes.

The other DFA techniques available in this package are demonstrated in Figure 11 and again in “Supplement to: Package AdvEMDpy: Algorithmic Variations of Empirical Mode Decomposition in $\textsf{Python}$ ” for easy reference. The implementation of Enhanced EMD (not called EEMD or E-EMD to avoid confusion with Ensemble EMD (EEMD)) can be found in section 10.4.2 of the supplement – this technique relies on derivatives to optimise, in some sense, the positioning of the “extrema” for the splines. Inflection point interpolation and binomial average interpolation can be found in sections 10.4.3 and 10.4.4, respectively – these two methods are more prone to discontinuities and rely on well-behaved time series, in some sense.

Figure 11 Example time series demonstrating five different local mean estimation techniques through detrended fluctuation analysis.

7.5 Stopping criteria in AdvEMDpy package

To prevent over-sifting resulting in physically meaningless IMFs and little discernible information about the process under observation, several stopping criteria are provided. The validity of the various stopping criteria warrants further study as some are more related to algorithmic steps than others. The base implementation of the stopping criteria is as follows:

imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,

knot_time=knot_time,

stop_crit=‘S_stoppage’,

stop_crit_threshold=10,

mft_theta_1=0.05,

mft_theta_2=0.5,

mft_alpha=0.05,

mean_threshold=10,

max_internal_iter=30,

max_imfs=10)[:3].

The numerous stopping criteria in this package can be found in the various subsections of section 10.5 in “Supplement to: Package AdvEMDpy: Algorithmic Variations of Empirical Mode Decomposition in $\textsf{Python}$ ”. They are not listed here as they are not the major contributions of this work besides the contribution of section 10.5 being the most complete review of available EMD stopping criteria. A final contribution is the recommendation in the following section which is the most robust stopping criterion available which manages the stopping condition with a logical stopping process (using extrema count variables), as well as managing the total number of external siftings (number of IMFs), the total number of internal siftings, and the modified mean threshold in Modified Condition 2.

7.6 Spline methods in AdvEMDpy package

All the above techniques are listed with the cubic B-spline implementation of EMD in mind. Other spline techniques are effective when using EMD. The other splines used have a different basis, but because they are both cubic bases, they may be easily mapped onto one another. The base implementation of cubic B-spline EMD is:

imfs, hts, ifs = \

emd.empirical_mode_decomposition(knots=knots, knot_time=knot_time,

spline_method=‘b_spline’)[:3].

Cubic Hermite spline interpolation and Akima spline interpolation are also available in this package and are discussed, along with the differences in their respective construction when compared against cubic B-splines, in sections 10.6.1 and 10.6.2, respectively, of “Supplement to: Package AdvEMDpy: Algorithmic Variations of Empirical Mode Decomposition in $\textsf{Python}$ ”.

7.7 Discrete-time Hilbert transforms in AdvEMDpy package

Despite B-splines having closed-form solutions to the HT, discrete solutions to the HT may be needed for a number of reasons. Here are two widely known and used methods that are included. The base implementation of the DTHT is as follows:

imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,

knot_time=knot_time,

dtht=False,

dtht_method=‘fft’)[:3].

The basic DTHT (title from cited literature) and the fast-Fourier transform DTHT (FFT DTHT) can be found in more detail in sections 10.7.1 and 10.7.2, respectively, of “Supplement to: Package AdvEMDpy: Algorithmic Variations of Empirical Mode Decomposition in $\textsf{Python}$ ”. As discussed in the following section, it is often advantageous to compare this output with the closed-form solution and should be output during initial debugging and exploratory investigation of the algorithm.

7.8 Knot point optimisations in AdvEMDpy package

Two methods are available for the optimisation of knot point allocation. Both methods are in the bisection family of knot optimisation techniques. The first technique simply bisects the domain iteratively until some error bound is met. The next method is a variation on this that extends the domain or diminishes the domain based on some error bound. The base implementation of the knot point optimisation is as follows:

imfs, hts, ifs = \

emd.empirical_mode_decomposition(knots=knots,

knot_time=knot_time,

optimise_knots=0,

knot_method=‘ser_bisect’,

knot_error=10,

knot_lamda=1,

knot_epsilon=0.5)[:3].

In section 10.8.1 of “Supplement to: Package AdvEMDpy: Algorithmic Variations of Empirical Mode Decomposition in $\textsf{Python}$ ”, classical bisection is discussed, before the more computationally expensive serial bisection technique is discussed in section 10.8.2 of the supplementary materials. As well as including the formal mathematics, these sections also include descriptive figures to assist with explainability and readability.

7.9 Ensemble empirical mode decomposition in AdvEMDpy package

This is a noise-assisted data analysis (NADA) technique introduced in Wu & Huang (Reference Wu and Huang2009) that utilises white noise. This technique relies on the ability of the algorithm to discern consistent structures in the presence of different sets of randomised white noise. The base implementation of the EEMD is implemented as follows:

imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,

knot_time=knot_time,

ensemble=False,

ensemble_sd=0.5,

ensemble_iter=10)[:3].

The original time series is median filtered to approximate the original level of noise in the system. The median is removed from the time series before the standard deviation is calculated, $ \sigma_s $ . The noise added to the system has a mean of zero and a standard deviation of $ \sigma_s \times \sigma_e $ , where $ \sigma_e $ (ensemble_sd=0.5) is the chosen level of standard deviation in the ensemble system as a fraction of the original standard deviation in the system. The IMFs are calculated for this modified time series and stored before a new set of noise is added to the original time series and the sifting process is repeated – the procedure is repeated a predetermined number of iterations (ensemble_iter=10).

Further applications of this technique will benefit from random additions of other colours of noise such as violet noise, blue noise, pink noise, and red noise (also known as Brownian noise), which all have non-constant power spectral densities, compared with white noise that has a constant power spectral density. This technique, with the random addition of other colours of noise, is called full-spectrum ensemble empirical mode decomposition (FSEEMD). This technique is experimental and is also included with the package in the script entitled emd_experimental.py for completeness sake.

This method performs a similar task, utilising a different methodology, to ICA-EMD, introduced in van Jaarsveldt et al. (Reference van Jaarsveldt, Peters, Ames and Chantler2021), which is finding the most consistent structures among noisy data. EEMD proceeds by introducing noise to the time series and isolating IMFs from the noisy time series. This task is repeated several times and IMFs are simply averaged. A more sophisticated sorting technique such as the minimum description length principle (introduced in Fayyad & Irani (Reference Fayyad and Irani1993)) should be used to ensure structures of similar frequency are grouped as in van Jaarsveldt et al. (Reference van Jaarsveldt, Peters, Ames and Chantler2021). ICA-EMD proceeds by applying independent component analysis (ICA) directly to all the noisy IMFs to isolate the most consistent structure, before EMD is performed on the ICA component to isolate IMFs.

8.1 United Kingdom of Great Britain and Northern Ireland birth rates

Trends in birth rates occur for a variety of reasons such as global events (World War 2), national governmental intervention (such as the NHS providing the contraceptive pill and the Abortion Act), and societal trends (such as people delaying childbearing years to later ages). These structures can be seen in the plotted time series shown in Figure 12 along with numerous other structures. An annual cycle is faintly visible in Figure 12 with the annual structure representing the first IMF in the EMD of the data in Figure 12 and available here: Mortality Database (2022).

Figure 12 Monthly births in the United Kingdom from January 1938 until December 2018 inclusive, with some notable features and causality.

The IF of the first IMF can be seen in Figure 13. Little analysis can be performed on this IF as noise is still present. The highest intensities can be seen every 9 months at a few points in the previous century. It is, however, clear that some structure is present. By using EMD-X11, first proposed in van Jaarsveldt et al. (Reference van Jaarsveldt, Peters, Ames and Chantler2021), one can improve the resolution of the IF – this can be seen in Figure 14.

Figure 13 IF of first IMF extracted using EMD of monthly births in the United Kingdom from January 1938 until December 2018 inclusive, with some features and limits included.

Figure 14 IF of first IMF extracted using EMD and refined using X11 of monthly births in the United Kingdom from January 1938 until December 2018 inclusive, with some features and limits included.

Better analysis can be performed on Figure 14. After World War 2 and before the implementation of the Abortion Act (1968) and the NHS providing of contraceptive pills to unmarried women, there is a clear increasing in frequency and intensity structure. After the contraceptive pill became available to the general public, this structure began to depreciate until it finally stabilised at the beginning of the 1990s. After this period, a clear annual structure is observed with only a minor break in the 2000s during the most recent surge in birth rates. The increased resolution of EMD-X11 has allowed improved analysis of the birth rates.

8.2 United Kingdom of Great Britain and Northern Ireland mortality trends

In Figures 15 and 16, the total number of deaths in the United Kingdom can be seen for women and men, respectively. One of the most prevalent structures in both figures is the exponential decrease in the infant (aged 0–5 non-inclusive) mortality rate over the past century. This was part of a major campaign by many European countries and others worldwide to decrease infant mortality. Another, less obvious, almost meta-structure can be seen as the increasing of year and deaths in the maxima as a result of the significant increase in births after World War 1. This meta-structure would be more difficult to analyse with most modern techniques, but EMD is able to iteratively decompose the structures and associate these maxima.

Figure 15 Annual deaths among females in the United Kingdom from 1922 until 2020 inclusive, in a 5-year stratification.

Figure 16 Annual deaths among males in the United Kingdom from 1922 until 2020 inclusive, in a 5-year stratification.

In Figures 17, 18, 19, and 20, the first IMF for each of the five year stratifications over the age of 50 are plotted. In Figure 17, the cyclical nature of this mortality data is clear, but the true underlying driving-force behind this structure is obscured. By plotting the deaths based on when people were born, as in Figure 18, the fewer people dying who were born between 1915–1920 is as a result of the World War 1 killing a statistically significant number of people in the general population. This is followed by the post-World War 1 baby boom in the 1920s. The analogous data for men in the United Kingdom are plotted in Figures 19 and 20 – the same structures are observable in these figures as can be expected, but the amplitude of the World War 1 event structure among males is more pronounced owing to the statistically more men that died during World War 1.

Figure 17 IMF 1 of deaths among females in the United Kingdom from 1922 until 2020 inclusive, in a 5-year stratification.

Figure 18 IMF 1 of deaths among females in the United Kingdom from 1922 until 2020 inclusive, in a 5-year stratification, sorted by birth year.

Figure 19 IMF 1 of deaths among males in the United Kingdom from 1922 until 2020 inclusive, in a 5-year stratification.

Figure 20 IMF 1 of deaths among males in the United Kingdom from 1922 until 2020 inclusive, in a 5-year stratification, sorted by birth.