## 3 Materials and Methods
### Multivariate Empirical Mode Decomposition (MEMD)  
The fundamental approach of EMD involves separating the signal into natural IMFs [31]. Upon obtaining the IMFs, the application of the HHT makes it possible to determine the amplitude, phase, and frequency of the present signal, enabling a time-frequency analysis of the observed phenomenon. This approach has been successfully utilized in analyzing the dynamic behavior of power systems [32; 33; 34; 35]. In the end, the EMD process separates the signal into a specific number of IMFs and a residue. The total of all IMFs combined with the residual component enables the restoration of the original signal without any loss of data or alteration. IMFs possess distinct characteristics: a local equilibrium near zero, all peaks (or troughs) being higher (or lower) than zero, and the count of extreme points and zeros being either identical or differing by no more than one. The average of the envelopes formed by local maximums and those formed by local minimums is zero. Every IMF, regardless of being linear or nonlinear, embodies basic oscillations with an equal count of extreme points and crossings of the zero line.  
\[m_{1}(t)=\frac{x_{up}(t)+x_{low}(t)}{2} \tag{4}\]  
where \(x_{up}\) signifies the upper envelope and \(x_{low}\) indicates the lower envelope; the disparity between the original signal \(x(t)\) and the mean value of the signal envelope \(m_{1}(t)\) results in \(h_{1}\):  
\[h_{1}=x(t)-m_{1} \tag{5}\]  
By iterating \(k\) times, \(h_{1k}\) becomes an IMF:  
\[h_{1k}=h_{1(k-1)}-m_{1k} \tag{6}\]  
When the specified criteria are fulfilled [36], \(h_{1k}\) transforms into the initial IMF:  
\[c_{1}=h_{1k} \tag{7}\]The original signal \(x(t)\) is reduced by the acquired IMF function, and subsequent IMFs are derived from the residual signals.  
\[r_{n}=r_{n-1}-c_{n-1} \tag{8}\]  
The iterative process concludes when one of the conditions is met: either the component \(c_{n}\) or the residue \(r_{n}\) diminishes to a value below a predetermined threshold, or when the residue \(r_{n}\) turns into a monotonic function, making it impossible to extract any further IMF functions. Regardless of whether the original signal has a zero mean value, the final residue may not be zero. In cases where the signal shows a tendency towards growth or decline, this final residue serves as an indicator of the signal's overall trend. The original signal is obtained by summing all IMF functions and the residue [36]:  
\[x(t)=\sum_{j=1}^{n}c_{j}+r_{n} \tag{9}\]  
An extended version of standard EMD, known as multivariate EMD (MEMD), is developed for multivariate data. As highlighted by [37], MEMD exhibits properties such as coherence and a data fusion algorithm. In this approach, data are scrutinized within a p-dimensional domain, estimating the local p-dimensional mean by averaging multiple envelopes.  
In the MEMD procedure, the first step is to generate a series of directional vectors, \(K\), by utilizing a point set on a p-dimensional hypersphere generated with a low-discrepancy Hammersley sequence. The procedure starts by initializing the residue of the p-dimensional signal as \(r(t)=x(t)\) and \(i=1\). After that, the set of projections along each vector \(K\) is obtained:  
\[\left\{e_{\theta_{k}}(t)\right\}_{k=1}^{K} \tag{10}\]  
The next step is to identify the local maxima and minima within \(\left\{e_{\theta_{k}}(t)\right\}_{k=1}^{K}\) by interpolating the local maxima and minima with a cubic spline. Over then determines the average of the upper and lower envelopes by calculating their mean:  
\[E_{mean}(t)=\frac{1}{2K}\sum_{k=1}^{K}\left(E_{min}^{\theta_{k}}(t)+E_{max}^{ \theta_{k}}(t)\right) \tag{11}\]  
The subsequent phase is to update \(r(t)=r(t)-E_{mean}(t)\) and repeat steps until a stopping criterion is met [36]. The process involves extracting the \(i\)th IMF denoted as \(E_{mean}(t)\) and subsequently updating residual \(r(t)\) from the original signal \(x(t)\) as:  
\[r(t)=x(t)-E_{mean}(t) \tag{12}\]