Content deleted Content added
ce |
→BPBEMD algorithm: move ref out of heading |
||
| (38 intermediate revisions by 26 users not shown) | |||
Line 1:
In [[signal processing]],
==Motivation==
Multidimensional empirical mode decomposition is a popular method because of its applications in many fields, such as texture analysis, financial applications, [[Digital image processing|image processing]], [[Ocean Engineering|ocean engineering]], [[Seismology|seismic]] research,
==Introduction to empirical mode decomposition (EMD)==
[[File:Flow chart of EMD algorithm.jpg|thumb|400x400px|Flow chart of basic EMD algorithm<ref>{{Cite journal|url=http://www.ripublication.com/irph/ijeee_spl/ijeeev7n8_14.pdf|author=Sonam Maheshwari |author2=Ankur Kumar |title=Empirical Mode Decomposition: Theory & Applications |journal=International Journal of Electronic and Electrical Engineering |issn=0974-2174 |volume=7 |issue=8 |year=2014 |pages=873–878}}</ref>{{Predatory open access publisher}}]]
The
The EMD method was developed so that
The EMD method decomposes the input signal into
: <math>I(n)=\sum_{m=1}^M \operatorname{IMF}_m(n)+\operatorname{Res}_M(n)</math>
where <math>I(n)</math> is the multi-component signal. <math>\operatorname{IMF}_m(n)</math> is the <math>M^\text{th}</math> intrinsic mode function, and <math>\operatorname{Res}_M(n)</math> represents the residue corresponding to <math>M</math> intrinsic modes.
==Ensemble empirical mode decomposition==
In the case of only one observation, one of the multiple-observation ensembles is mimicked by adding
The EEMD consists of the following steps:
Line 26 ⟶ 25:
# Adding a white noise series to the original data.
# Decomposing the data with added white noise into oscillatory components.
# Repeating step 1 and step 2 again and again, but with a different white noise series added each time.
# Obtaining the
In these steps, EEMD uses two properties of white noise:
# The added white noise leads to a relatively even distribution of extrema distribution on all timescales.
# The [[Dyadic transformation|dyadic]] filter bank property provides a control on the periods of oscillations contained in an oscillatory component, significantly reducing the chance of scale mixing in a component. Through ensemble average, the added noise is averaged out.<ref name=":9" />
=== Pseudo-bi-dimensional empirical mode decomposition ===
Source:<ref name=":5" /> It should be pointed out here that the “pseudo-BEMD” method is not limited to only 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 reflects better 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 [[Hilbert huang transform|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 i × j × k elements, the pseudo-BEMD will yield detailed 3D components of m × n × q where m, n and q are the number of the IMFs decomposed from each dimension having i, j, and k elements, respectively.▼
▲To design a pseudo-BEMD algorithm the key step is to translate the algorithm of the 1D [[Hilbert huang transform|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 <math>i
Mathematically let us represent a 2D signal in the form of ixj matrix form with a finite number of elements.▼
▲Mathematically let us represent a 2D signal in the form of
:<math>
Line 46 ⟶ 47:
\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}
</math><ref name=":5" />
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
: <math>
\begin{align}
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,
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,
\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,
\end{align}
</math><ref name=":5" />
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 <math>X(i,j)=\sum_{ k \mathop =1}^mRX(k,i,j)</math>
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.
Line 64 ⟶ 65:
For example, the component
# RX(1,i,j) will be decomposed into CRX(1,1,i,j), CRX(1,2,i,j),
# RX(2,i,j) will be decomposed into CRX(2,1,i,j), CRX(2,2,i,j),
# RX(m,i,j) will be decomposed into CRX(m,1,i,j), CRX(m,2,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
: <math>
Line 76 ⟶ 77:
</math><ref name=":5" />
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.
Line 88 ⟶ 89:
given as <math>I=f(x_1,x_2,x_3,x_4,\ldots,x_n)</math>
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
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
Line 102 ⟶ 103:
</math>
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.
Source:<ref name=":7" />
▲=== Multi-dimensional ensemble empirical mode decomposition.<ref name=":7" /> ===
A Fast and efficient data analysis is very important for large sequences hence the MDEEMD focuses on two important things
Line 112 ⟶ 115:
# 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.
==== Principal component analysis (PCA) or empirical orthogonal function analysis (EOF)
The [[principal component analysis]]/[[Empirical orthogonal functions|empirical orthogonal function]] analysis (PCA/EOF) has been widely used in data analysis and image compression, its main objective is to reduce a data set containing a large number of variables to a data set containing fewer variables, but that still represents a large fraction of the variability contained in the original data set. In climate studies, EOF analysis is often used to study possible spatial modes (i.e., patterns) of variability and how they change with time . In statistics, EOF analysis is known as [[principal component analysis]] (PCA).
Typically, the EOFs are found by computing the eigenvalues and eigen vectors of a spatially weighted anomaly covariance matrix of a field. Most commonly, the spatial weights are the cos(latitude) or, better for EOF analysis, the sqrt(cos(latitude)). The derived eigenvalues provide a measure of the percent variance explained by each mode. Unfortunately, the eigenvalues are not necessarily distinct due to sampling issues. North et al. (Mon. Wea. Rev., 1982, eqns
Atmospheric and oceanographic processes are typically 'red' which means that most of the variance (power) is contained within the first few modes. The time series of each mode (aka, principle components) are determined by projecting the derived eigen vectors onto the spatially weighted anomalies. This will result in the amplitude of each mode over the period of record.
By construction, the EOF patterns and the principal components are independent. Two factors inhibit physical interpretation of EOFs: (i) The orthogonality constraint and (ii) the derived patterns may be ___domain dependent. Physical systems are not necessarily orthogonal and if the patterns depend on the region used they may not exist if the ___domain changes.
==== Spatial-temporal signal using multi-dimensional ensemble empirical mode decomposition
Assume, we have a spatio-temporal data
Using PCA/EOF, one can express
where
If the data subjected to PCA/EOF analysis is all white noise, all eigenvalues are theoretically equal and there is no preferred vector direction for the principal component in PCA/EOF space. To retain most of the information of the data, one needs to retain almost all the PC's and EOF's, making the size of PCA/EOF expression even larger than that of the original but If the original data contain only one spatial structure and oscillate with time, then the original data can be expressed as the product of one PC and one EOF, implying that the original data of large size can be expressed by small size data without losing information, i.e. highly compressible.
The variability of a smaller region tends to be more spatio-temporally coherent than that of a bigger region containing that smaller region, and, therefore, it is expected that fewer PC/EOF components are required to account for a threshold level of variance hence one way to improve the efficiency of the representation of data in terms of PC/EOF component is to divide the global spatial ___domain into a set of sub-regions. If we divide the original global spatial ___domain into n sub-regions containing N1, N2, . . . , Nn spatial grids, respectively, with all Ni, where i=1, . . . , n, greater than M, where M denotes the number of temporal locations, we anticipate that the numbers of the retained PC/EOF pairs for all sub-regions K1, K2, . . . , Kn are all smaller than K, the total number of data values in PCA/EOF representation of the original data of the global spatial ___domain by the equation given up is K×(N+M). For the new approach of using spatial division, the total number of values in PCA/EOF
: <math> \sum_{i=1}^nK_i(M+N_i)=K'_iN+M\sum_{i=1}^nK_i</math>
where
: <math>K'_i=\frac{\sum_{i=1}^nK_iN_i}{N}. </math>
Therefore, the compression rate of the spatial ___domain is as follows
Line 143 ⟶ 146:
The advantage of this algorithm is that an optimized division and an optimized selection of PC/EOF pairs for each region would lead to a higher rate of compression and result into significantly lower computation as compared to a Pseudo BEMD extended to higher dimensions.
=== Fast multidimensional ensemble empirical mode decomposition ===
Source:<ref name=":7" /> 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
The fast MEEMD includes the following steps:
Line 159 ⟶ 164:
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.<br />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.
== Parallel implementation of multi-dimensional ensemble empirical mode decomposition.
In MEEMD, although ample parallelism potentially exists in the ensemble dimensions and/or the non-operating dimensions, several challenges still face a high performance MEEMD implementation.<ref name=":8" />
[[File:Sample_BEMD_Simulation_results_for_a_noisy_signal.jpg|thumb|Bi-Dimensional EMD corrupted with Noise]]
Line 165 ⟶ 170:
# Dynamic data variations: In EEMD, white noises change the number of extrema causing some irregularity and load imbalance, and thus slowing down the parallel execution.
# Stride memory accesses of high-dimensional data: High dimensional data are stored in non-continuous memory locations. Accesses along high dimensions are thus strided and uncoalesced, wasting available memory bandwidth.
# Limited resources to harness parallelism: While the independent EMDs and/or EEMDs comprising an MEEMD provide high parallelism, the computational capacities of multi-core and many-core processors may not be sufficient to fully exploit the inherent parallelism of MEEMD. Moreover, increased parallelism may increase memory requirements beyond the memory capacities of these processors.
[[File:Sample_BEMD_Simulation_results_for_a_noisy_signal_with_imf.jpg|thumb|Bi-Dimensional EMD Intrinsic mode function along with the residue eliminating the noise level.]] In MEEMD, when a high degree of parallelism is given by the ensemble dimension and/or the non-operating dimensions, the benefits of using a thread-level parallel algorithm are threefold.<ref name=":8" /> # It can exploit more parallelism than a block-level parallel algorithm.
Line 171 ⟶ 178:
# Its implementation is like the sequential one, which makes it more straightforward.
===
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.
=== CUDA implementation
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.<ref name=":8" />
Line 182 ⟶ 189:
'''The fast and adaptive bidimensional empirical mode decomposition (FABEMD)''' is an improved version of traditional BEMD.<ref name=":0" /> 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
===FABEMD algorithm===
Source:<ref name=":0" /> ;Step 1 – Determine and detect local maximum and minimum
Line 213 ⟶ 221:
: <math>\ell = n - \frac{w_{ex} - 1}{2}:n + \frac{w_{ex} - 1}{2},(\ell \ne n)</math><ref name=":0" />
;[[File:Flow chart for FABEMD algorithm.jpg|thumb|upright|Flow chart for FABEMD algorithm<ref>{{Cite book |doi=10.1109/ICASSP.2008.4517859|chapter=A novel approach of fast and adaptive bidimensional empirical mode decomposition|title=2008 IEEE International Conference on Acoustics, Speech and Signal Processing|pages=1313–1316|year=2008|last1=Bhuiyan|first1=Sharif M. A.|last2=Adhami|first2=Reza R.|last3=Khan|first3=Jesmin F.|isbn=978-1-4244-1483-3|s2cid=18226941 }}</ref>]]Step 2 – Obtain the size of window for order-statistic filter
At first, we define <math>d_{\mathrm{adj}-\max}</math> and <math>d_{\mathrm{adj}-\min}</math> 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, <math>d_{\mathrm{adj}-\max}</math> and <math>d_{\mathrm{adj}-\min}</math> 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:
Line 240 ⟶ 248:
===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===
Line 248 ⟶ 256:
===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<ref name=":2" /> 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===
Source:<ref name=":2" /> ;Step 1 – Extend super diffusion model from 1-D to 2-D
Line 265 ⟶ 274:
where <math>{g_{q,i}}(x)</math> represent the qth-order stopping function in direction i.
Then, based on the [[Navier–Stokes equations/Derivation|
: <math>u_t(x,t) = \operatorname{div}(\alpha {G_1}\nabla u(x,t) - (1 - \alpha ){G_2} \nabla \Delta u(x,t))</math><ref name=":2" />
Line 309 ⟶ 318:
===Concept===
There are some problems in BEMD and boundary extending implementation in the iterative sifting process, including time
===BPBEMD algorithm
The few core steps for BPBEMD algorithm are:<ref name=":3" />
;Step 1
Line 335 ⟶ 344:
===Advantages===
This method can process larger number of elements than traditional BEMD method. Also, it can shorten the time
===Limitations===
Because most of image inputs are non-stationary which
==Applications==
In the first part, these MEEMD techniques can be used on Geophysical data sets such as climate, magnetic, seismic data variability which takes advantage of the fast algorithm of MEEMD. The MEEMD is often used for nonlinear geophysical data filtering due to its fast algorithms and its ability to handle large amount of data sets with the use of compression without losing key information. The IMF can also be used as a signal enhancement of Ground Penetrating Radar for nonlinear data processing; it is very effective to detect geological boundaries from the analysis of field anomalies.<ref name=":6" />
In the second part, the PDE-based MEMD and FAMEMD can be implemented on audio processing, image processing and texture analysis. Because of its several properties, including stability, less time
==References==
Line 349 ⟶ 358:
<ref name=":9">N.E. Huang, Z. Shen, et al., "The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non- Stationary Time Series Analysis," Proceedings: Mathematical, Physical and Engineering Sciences, vol. 454, pp. 903–995, 1998.</ref>
<ref name=":5">Chih-Sung Chen, Yih Jeng,"Two-dimensional nonlinear geophysical data filtering using the multidimensional EEMD method", Department of Earth Sciences, National Taiwan Normal University, 88, Sec. 4, Ting-Chou Road, Taipei, 116, Taiwan, ROC</ref>
<ref name=":7">Wu Z, Feng J, Qiao F, Tan Z-M, "2016 Fast multidimensional ensemble empirical mode decomposition for the analysis of big spatio-temporal datasets.", Phil. Trans. R. Soc. A 374: 20150197.</ref>
| |||