Pseudo-bi-dimensional empirical mode decomposition^[3]

The “pseudo-BEMD” method is not limited to one-spatial dimension; rather, it can be applied to data of any number of spatial-temporal dimensions. Since the spatial structure is essentially determined by timescales of the variability of a physical quantity at each location and the decomposition is completely based on the characteristics of individual time series at each spatial location, there is no assumption of spatial coherent structures of this physical quantity. When a coherent spatial structure emerges, it better reflects the physical processes that drive the evolution of the physical quantity on the timescale of each component. Therefore, we expect this method to have significant applications in spatial-temporal data analysis.

To design a pseudo-BEMD algorithm the key step is to translate the algorithm of the 1D EMD into a Bi-dimensional Empirical Mode Decomposition (BEMD) and further extend the algorithm to three or more dimensions which is similar to the BEMD by extending the procedure on successive dimensions. For a 3D data cube of ${\displaystyle i\times j\times k}$ elements, the pseudo-BEMD will yield detailed 3D components of ${\displaystyle m\times n\times q}$ where ${\displaystyle m}$, ${\displaystyle n}$ and ${\displaystyle q}$ are the number of the IMFs decomposed from each dimension having ${\displaystyle i}$, ${\displaystyle j}$, and ${\displaystyle k}$ elements, respectively.

Mathematically let us represent a 2D signal in the form of ${\displaystyle i\times j}$ matrix form with a finite number of elements.

${\displaystyle X(i,j)={\begin{bmatrix}x(1,1)&x(1,2)&\cdots &x(1,j)\\x(2,1)&x(2,2)&\cdots &x(1,j)\\\vdots &\vdots &&\vdots \\x(i,1)&x(i,2)&\cdots &x(i,j)\end{bmatrix}}}$^[3]

At first we perform EMD in one direction of X(i, j), Row wise for instance, to decompose the data of each row into m components, then to collect the components of the same level of m from the result of each row decomposition to make a 2D decomposed signal at that level of m. Therefore, m set of 2D spatial data are obtained

${\displaystyle {\begin{aligned}RX(1,i,j)&={\begin{bmatrix}x(1,1,1)&x(1,1,2)&\cdots &x(1,1,j)\\x(1,2,1)&x(1,2,2)&\cdots &x(1,2,j)\\\vdots &\vdots &&\vdots \\x(1,i,1)&x(1,i,2)&\cdots &x(1,i,j)\end{bmatrix}}\\[6pt]RX(2,i,j)&={\begin{bmatrix}x(2,1,1)&x(2,1,2)&\cdots &x(2,1,j)\\x(2,2,1)&x(2,2,2)&\cdots &x(2,2,j)\\\vdots &\vdots &&\vdots \\x(2,i,1)&x(2,i,2)&\cdots &x(2,i,j)\end{bmatrix}}\\\vdots \qquad &\,\,\,\vdots \qquad \vdots \\[6pt]RX(m,i,j)&={\begin{bmatrix}x(m,1,1)&x(m,1,2)&\cdots &x(m,1,j)\\x(m,2,1)&x(m,2,2)&\cdots &x(m,2,j)\\\vdots &\vdots &&\vdots \\x(m,i,1)&x(m,i,2)&\cdots &x(m,i,j)\end{bmatrix}}\end{aligned}}}$^[3]

where RX (1, i, j), RX (2, i, j), and RX (m, i, j) are the m sets of signal as stated (also here we use R to indicate row decomposing). The relation between these m 2D decomposed signals and the original signal is given as ${\displaystyle X(i,j)=\sum _{k\mathop {=} 1}^{m}RX(k,i,j)}$.^[3]

The first row of the matrix RX (m, i, j) is the mth EMD component decomposed from the first row of the matrix X (i, j). The second row of the matrix RX (m, i, j) is the mth EMD component decomposed from the second row of the matrix X (i, j), and so on.

Suppose that the previous decomposition is along the horizontal direction, the next step is to decompose each one of the previously row decomposed components RX(m, i, j), in the vertical direction into n components. This step will generate n components from each RX component.

For example, the component

1. RX(1,i,j) will be decomposed into CRX(1,1,i,j), CRX(1,2,i,j),...,CRX(1,n,i,j)
2. RX(2,i,j) will be decomposed into CRX(2,1,i,j), CRX(2,2,i,j),..., CRX(2,n,i,j)
3. RX(m,i,j) will be decomposed into CRX(m,1,i,j), CRX(m,2,i,j),..., CRX(m,n,i,j)

where C means column decomposing. Finally, the 2D decomposition will result into m× n matrices which are the 2D EMD components of the original data X(i,j). The matrix expression for the result of the 2D decomposition is

${\displaystyle CRX(m,n,i,j)={\begin{bmatrix}crx(1,1,i,j)&crx(2,1,i,j)&\cdots &crx(m,1,i,j)\\crx(1,2,i,j)&crx(2,2,i,j)&\cdots &crx(m,2,i,j)\\\vdots &\vdots &&\vdots \\crx(1,n,i,j)&crx(2,n,i,j)&\cdots &crx(m,n,i,j)\end{bmatrix}}}$^[3]

where each element in the matrix CRX is an i × j sub-matrix representing a 2D EMD decomposed component. We use the arguments (or suffixes) m and n to represent the component number of row decomposition and column decomposition, respectively rather than the subscripts indicating the row and the column of a matrix. Notice that the m and n indicate the number of components resulting from row (horizontal) decomposition and then column (vertical) decomposition, respectively.

By combining the components of the same scale or the comparable scales with minimal difference will yield a 2D feature with best physical significance. The components of the first row and the first column are approximately the same or comparable scale although their scales are increasing gradually along the row or column. Therefore, combining the components of the first row and the first column will obtain the first complete 2D component (C2D1). The subsequent process is to perform the same combination technique to the rest of the components, the contribution of the noises are distributed to the separate component according to their scales. As a result, the coherent structures of the components emerge, In this way, the pseudo-BEMD method can be applied to reveal the evolution of spatial structures of data.

${\displaystyle C2D_{L}=\sum _{k=1}^{m}crx_{k,l}+\sum _{k=l+1}^{n}crx_{l,k}}$^[3]

Following the convention of 1D EMD, the last component of the complete 2D components is called residue.

The decomposition scheme proposed here could be extended to data of any dimensions such as data of a solid with different density or other measurable properties

given as ${\displaystyle I=f(x_{1},x_{2},x_{3},x_{4},\ldots ,x_{n})}$

In which the subscription, n, indicated the number of dimensions. The procedure is identical as stated above: the decomposition starts with the first dimension, and proceeds to the second and third until all the dimensions are exhausted. The decomposition is still implemented by slices. This new approach is based on separating the original data into one-dimensional slices, then applying ensemble EMD to each one-dimensional slice. The key part of the method is in the construction of the IMF according to the principle of combination of the comparable minimal scale components.

For example, the matrix expression for the result of a 3D decomposition is TCRX(m,n,q,i,j,k) where T denotes the depth (or time) decomposition. Based on the comparable minimal scale combination principle as applied in the 2D case, the number of complete 3D components will be the smallest value of m, n, and q. The general equation for deriving 3D components is

${\displaystyle C3D_{L}=\sum _{m=l}^{m}\sum _{n=\ell }^{n}tcrx_{m,n,\ell }+\sum _{m=\ell +1}^{m}\sum _{q=\ell +1}^{n}tcrx_{m,\ell ,q}+\sum _{n=\ell +1}^{m}\sum _{q=\ell +1}^{m}tcrx_{\ell ,n,q}}$  ^[3]

where ℓ denotes the level of the C3D i.e.

${\displaystyle {\begin{cases}\ell =m=n&{\text{if }}m=n;\\\ell =m&{\text{if }}m<n;\\\ell =n&{\text{if }}m>n.\end{cases}}}$

The pseudo-BEMD method has several advantages. For instance, the sifting procedure of the pseudo-BEMD is a combination of one dimensional sifting. It employs 1D curve fitting in the sifting process of each dimension, and has no difficulty as encountered in the 2D EMD algorithms using surface fitting, which has the problem of determining the saddle point as a local maximum or minimum. Sifting is the process which separates the IMF and repeats the process until the residue is obtained. The first step of performing sifting is to determine the upper and lower envelopes encompassing all the data by using the spline method. Sifting scheme for pseudo-BEMD is like the 1D sifting where the local mean of the standard EMD is replaced by the mean of multivariate envelope curves.

The major disadvantage of this method is that although we could extend this algorithm to any dimensional data we only use it for Two dimension applications. Because the computation time of higher dimensional data would be proportional to the number of IMF's of the succeeding dimensions. Hence, it could exceed the computation capacity for a Geo-Physical data processing system when the number of EMD in the algorithm is large. Hence, we have mentioned below faster and better techniques to tackle this disadvantage.

Multi-dimensional ensemble empirical mode decomposition.^[4]

A Fast and efficient data analysis is very important for large sequences hence the MDEEMD focuses on two important things

1. Data compression which involves decomposing data into simpler forms.
2. EEMD on the compressed data; this is the most challenging since on decomposing the compressed data there is a high probability to lose key information. A data compression method that uses principal component analysis (PCA)/empirical orthogonal function (EOF) analysis or principal oscillation pattern analysis is used to compress data.

Fast multidimensional ensemble empirical mode decomposition^[4]

For a temporal signal of length M, the complexity of cubic spline sifting through its local extrema is about the order of M, and so is that of the EEMD as it only repeats the spline fitting operation with a number that is not dependent on M. However, as the sifting number (often selected as 10) and the ensemble number (often a few hundred) multiply to the spline sifting operations, hence the EEMD is time-consuming compared with many other time series analysis methods such as Fourier transforms and wavelet transforms. The MEEMD employs EEMD decomposition of the time series at each division grids of the initial temporal signal, the EEMD operation is repeated by the number of total grid points of the domain. The idea of the fast MEEMD is very simple. As PCA/EOF-based compression expressed the original data in terms of pairs of PCs and EOFs, through decomposing PCs, instead of time series of each grid, and using the corresponding spatial structure depicted by the corresponding EOFs, the computational burden can be significantly reduced.

The fast MEEMD includes the following steps:

1. All pairs of EOF's, V_i, and their corresponding PCs, Y_i, of spatio-temporal data over a compressed sub-domain are computed.
2. The number of pairs of PC/EOF that are retained in the compressed data is determined by the calculation of the accumulated total variance of leading EOF/PC pairs.
3. Each PC Y_i is decomposed using EEMD, i.e.

${\displaystyle Y_{i}=\sum _{j=1}^{n}c_{j,i}+r_{n,i}}$^[4]

where c_j,i represents simple oscillatory modes of certain frequencies and r_n,i is the residual of the data Y_i. The result of the ith MEEMD component C_j is obtained as

${\displaystyle C_{j}=\sum _{j=1}^{40}c_{j,i}V_{i}.}$ ^[4]

In this compressed computation, we have used the approximate dyadic filter bank properties of EMD/EEMD.

Note that a detailed knowledge of the intrinsic mode functions of a noise corrupted signal can help in estimating the significance of that mode. It is usually assumed that the first IMF captures most of the noise and hence from this IMF we could estimate the Noise level and estimate the noise corrupted signal eliminating the effects of noise approximately. This method is known as denoising and detrending. Another advantage of using the MEEMD is that the mode mixing is reduced significantly due to the function of the EEMD.
The denoising and detrending strategy can be used for image processing to enhance an image and similarly it could be applied to Audio Signals to remove corrupted data in speech. The MDEEMD could be used to break down images and audio signals into IMF and based on the knowledge of the IMF perform necessary operations. The decomposition of an image is very advantageous for radar-based application the decomposition of an image could reveal land mines etc.

OpenMP implementation^[5]

The EEMDs comprising MEEMD are assigned to independent threads for parallel execution, relying on the OpenMP runtime to resolve any load imbalance issues. Stride memory accesses of high-dimensional data are eliminated by transposing these data to lower dimensions, resulting in better utilization of cache lines. The partial results of each EEMD are made thread-private for correct functionality. Memory requirements depend on the number of OpenMP threads and are managed by OpenMP runtime.

CUDA implementation^[5]

In the GPU CUDA implementation, each EMD, is mapped to a thread. The memory layout, especially of high-dimensional data, is rearranged to meet memory coalescing requirements and fit into the 128-byte cache lines. The data is first loaded along the lowest dimension and then consumed along a higher dimension. This step is performed when the Gaussian noise is added to form the ensemble data. In the new memory layout, the ensemble dimension is added to the lowest dimension to reduce possible branch divergence. The impact of the unavoidable branch divergence from data irregularity, caused by the noise, is minimized via a regularization technique using the on-chip memory. Moreover, the cache memory is utilized to amortize unavoidable uncoalesced memory accesses.^[5]

Concept

The fast and adaptive bidimensional empirical mode decomposition (FABEMD) is an improved version of traditional BEMD.^[6] The FABEMD can be used in many areas, including medical image analysis, texture analysis and so on. The order statistics filter can help in solving the problems of efficiency and restriction of size in BEMD.

Based on the algorithm of BEMD, the implementation method of FABEMD is really similar to BEMD, but the FABEMD approach just changes the interpolation step into a direct envelope estimation method and restricts the number of iterations for every BIMF to one. As a result, two order statistics, including MAX and MIN, will be used for approximating the upper and lower envelope. The size of the filter will depend on the maxima and minima maps obtained from the input. The steps of the FABEMD algorithm are listed below.

FABEMD algorithm^[6]

Step 1 – Determine and detect local maximum and minimum

As the traditional BEMD approach, we can find the jth ITS-BIMF ${\displaystyle F_{Tj}}$ of any source of input ${\displaystyle S_{i}}$ by neighboring window method. For FABEMD approach, we choose a different implementation approach.

From the input data, we can get an 2D matrix represent

${\displaystyle A=\left({\begin{array}{*{20}{c}}a_{11}&\cdots &a_{1N}\\\vdots &\ddots &\vdots \\a_{M1}&\cdots &a_{MN}\end{array}}\right)}$^[6]

where ${\displaystyle a_{MN}}$ is the element location in the matrix A, and we can define the window size to be ${\displaystyle w_{ex}\times w_{ex}}$. Thus, we can obtain the maximum and minimum value from the matrix as follows:

${\displaystyle a_{mn}\triangleq {\begin{cases}{\text{local max}}&{\text{if }}a_{mn}>a_{k\ell },\\{\text{local min}}&{\text{if }}a_{mn}<a_{k\ell },\end{cases}}}$^[6]

where

${\displaystyle k=m-{\frac {w_{ex}-1}{2}}:m+{\frac {w_{ex}-1}{2}},(k\neq m)}$^[6]

${\displaystyle \ell =n-{\frac {w_{ex}-1}{2}}:n+{\frac {w_{ex}-1}{2}},(\ell \neq n)}$^[6]

Flow chart for FABEMD algorithm^[7]

Step 2 – Obtain the size of window for order-statistic filter

At first, we define ${\displaystyle d_{\mathrm {adj} -\max }}$ and ${\displaystyle d_{\mathrm {adj} -\min }}$ to be the maximum and minimum distance in the array which is calculated from each local maximum or minimum point to the nearest nonzero element. Also, ${\displaystyle d_{\mathrm {adj} -\max }}$ and ${\displaystyle d_{\mathrm {adj} -\min }}$ will be sorted in descending order in the array according to the convenient selection. Otherwise, we will only consider square window. Thus, the gross window width will be as follows:

${\displaystyle w_{\max {\rm {\_en-g}}}=\operatorname {minimum} \{d_{\mathrm {adj} -\max }\}}$^[6]

${\displaystyle w_{\max {\rm {\_en-g}}}=\operatorname {maximum} \{d_{\mathrm {adj} -\max }\}}$^[6]

${\displaystyle w_{\min {\rm {\_en-g}}}=\operatorname {minimum} \{d_{\mathrm {adj} -\min }\}}$^[6]

${\displaystyle w_{\min {\rm {\_en-g}}}=\operatorname {maximum} \{d_{\mathrm {adj-\min } }\}}$^[6]

Step 3 – Apply order statistics and smoothing filters to obtain the MAX and MIN filter output

To obtain the upper and lower envelopes, there should be defined two parameter ${\displaystyle U_{Ej}(x,y)}$ and ${\displaystyle L_{Ej}(x,y)}$, and the equation will be as follows:

${\displaystyle U_{Ej}(x,y)=\max\{F_{Tj}(s,t)\}={\frac {1}{w_{sm}\times w_{sm}}}\sum _{(s,t)\in Z_{xy}}{U_{Ej}}(s,t)}$^[6]

${\displaystyle L_{Ej}(x,y)=\min\{F_{Tj}(s,t)\}={\frac {1}{w_{sm}\times w_{sm}}}\sum _{(s,t)\in Z_{xy}}L_{Ej}(s,t)}$^[6]

where ${\displaystyle Z_{xy}}$ is defined as the square region of window size, and ${\displaystyle w_{sm}}$ is the window width of the smoothing filter which ${\displaystyle w_{sm}}$ equals to ${\displaystyle w_{en}}$. Therefore, the MAX and MIN filter will form a new 2-D matrix for envelope surface which will not change the original 2-D input data.^[8]

Step 4 – Set up an estimation from upper and lower envelopes

This step is to make sure that the envelope estimation in FABEMD is nearly closed to the result from BEMD by using interpolation. In order to do comparison, we need to form corresponding matrices for upper envelope, lower envelope and mean envelope by using thin-plate spline surface interpolation to the max and min maps.

Advantages

This method (FABEMD) provides a way to use less computation to obtain the result rapidly, and it allows us to ensure more accurate estimation of the BIMFs. Even more, the FABEMD is more adaptive to handle the large size input than the traditional BEMD. Otherwise, the FABEMD is an efficient method that we do not need to consider the boundary effects and overshoot-undershoot problems.

Limitations

There is one particular problem that we will face in this method. Sometimes, there will be only one local maxima or minima element in the input data, so it will cause the distance array to be empty.

Concept

The Partial Differential Equation-Based Multidimensional Empirical Mode Decomposition (PDE-based MEMD) approach is a way to improve and overcome the difficulties of mean-envelope estimation of a signal from the traditional EMD. The PDE-based MEMD focus on modifying the original algorithm for MEMD. Thus, the result will provide an analytical formulation which can facilitate theoretical analysis and performance observation. In order to perform multidimensional EMD, we need to extend the 1-D PDE-based sifting process^[9] to 2-D space as shown by the steps below.

Here, we take 2-D PDE-based EMD as an example.

PDE-based BEMD algorithm^[9]

Step 1 – Extend super diffusion model from 1-D to 2-D

Considered a super diffusion matrix function

${\displaystyle G_{q}(x)=\left({\begin{array}{*{20}{c}}g_{q,1}(x)&0\\0&g_{q,2}(x)\end{array}}\right)}$^[9]

where ${\displaystyle {g_{q,i}}(x)}$ represent the qth-order stopping function in direction i.

Then, based on the Navier–Stokes equations, diffusion equation will be:

${\displaystyle u_{t}(x,t)=\operatorname {div} (\alpha {G_{1}}\nabla u(x,t)-(1-\alpha ){G_{2}}\nabla \Delta u(x,t))}$^[9]

where ${\displaystyle \alpha }$ is the tension parameter, and we assumed that ${\displaystyle q=2}$ .

Step 2 – Connect the relationship between diffusion model and PDEs on implicit surface

In order to relate to PDEs, the given equation will be

${\displaystyle u_{t}(x,t)=-(-1)^{q}\nabla _{S}^{2q}u(x,t)}$  ^[9]

where ${\displaystyle \nabla _{S}^{2q}}$ is the 2qth-order differential operator on u intrinsic to surface S, and the initial condition for the equation will be ${\displaystyle u(y,t)=f}$ for any y on surface S.

Step 3 – Consider all the numerical resolutions

To obtain the theoretical and analysis result from the previous equation, we need to make an assumption.

Assumption:

The numerical resolution schemes are assumed to be 4th-order PDE with no tension, and the equation for 4th-order PDE will be

${\displaystyle u_{t}=-\sum _{i,j=1}^{2}\partial _{x_{i}}^{1}(g_{i}\partial _{x_{i}}^{1}\partial _{x_{j}}^{2}u)}$  ^[9]

First of all, we will explicit scheme by approximating the PDE-based sifting process.

${\displaystyle U^{k+1}=\left(I-\Delta t\sum _{i,j=1}^{2}L_{ij}\right)U^{k}}$  ^[9]

where ${\displaystyle U}$ is a vector which consists of the value of each pixel, ${\displaystyle L_{ij}}$ is a matrix which is a difference approximation to the operator, and ${\displaystyle \Delta t}$ is a small time step.

Secondly, we can use additive operator splitting (AOS)^[10] scheme to improve the property of stability, because the small time step ${\displaystyle \Delta t}$ will be unstable when it comes to a large time step.

Finally, we can use the alternating direction implicit (ADI) scheme. By using ADI-type schemes, it is suggested that to mix a derivative term to overcome the problem that ADI-type schemes can only be used in second-order diffusion equation. The numerically solved equation will be :

${\displaystyle U^{k+1}=\left(\prod _{n=1}^{2}(I-\Delta tA_{nn})\right)^{-1}(I+\Delta t\sum _{i=1}^{2}\sum _{j\neq i}A_{ij})U^{k}}$^[9]

where ${\displaystyle A_{ij}}$ is a matrix which is the central difference approximation to the operator ${\displaystyle a_{ij}\partial _{x_{i}}^{1}\partial _{x_{j}}^{2}}$

Advantages

Based on the Navier–Stokes equations directly, this approach provides a good way to obtain and develop theoretical and numerical results. In particular, the PDE-based BEMD can work well for image decomposition fields. This approach can be applied to extract transient signal and avoiding the indeterminacy characterization in some signals.

Concept

There are some problems in BEMD and boundary extending implementation in the iterative sifting process, including time-consuming, shape and continuity of the edges, decomposition results comparison and so on. In order to fix these problems, the Boundary Processing in Bidimensional Empirical Decomposition (BPBEMD) method was created. The main points of the new method algorithm will be described next.

BPBEMD algorithm^[11]

The few core steps for BPBEMD algorithm are:

Step 1

Assuming the size of original input data and resultant data to be ${\displaystyle N\times N}$ and ${\displaystyle (N+2M)\times (N+2M)}$, respectively, we can also define that original input data matrix to be in the middle of resultant data matrix.

Step 2

Divide both original input data matrix and resultant data matrix into blocks of ${\displaystyle M\times M}$ size.

Step 3

Find the block which is the most similar to its neighbor block in the original input data matrix, and put it into the corresponding resultant data matrix.

Step 4

Form a distance matrix which the matrix elements are weighted by different distances between each block from those boundaries.

Step 5

Implement iterative extension when resultant data matrix faces a huge boundary extension, we can see that the block in original input data matrix is corresponding to the block in resultant data matrix.

Advantages

This method can process larger number of elements than traditional BEMD method. Also, it can shorten the time-consuming for the process. Depended on using nonparametric sampling based texture synthesis, the BPBEMD could obtain better result after decomposing and extracting.

Limitations

Because most of image inputs are non-stationary which don't exist boundary problems, the BPBEMD method is still lack of enough evidence that it is adaptive to all kinds of input data. Also, this method is narrowly restricted to be use on texture analysis and image processing.