# Jump Plus AM-FM Mode Decomposition

Mojtaba Nazari  Anders Rosendal Korshøj  Naveed ur Rehman  M. Nazari and N.
Rehman are with the Department of Electrical and Computer Engineering, Aarhus
University, 8200 Aarhus N, Denmark e-mails<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>Rosendal Korshøj is with 1: Aarhus University
Hospital, Dept. of Neurosurgery, Aarhus, Denmark, 2: Aarhus University, Dept.
of Clinical Medicine, Aarhus, Denmark e-mail<EMAIL_ADDRESS>
###### Abstract

A novel method for decomposing a nonstationary signal into amplitude- and
frequency-modulated (AM-FM) oscillations and discontinuous (jump) components
is proposed. Current nonstationary signal decomposition methods are designed
to either obtain constituent AM-FM oscillatory modes or the discontinuous and
residual components from the data, separately. Yet, many real-world signals of
interest simultaneously exhibit both behaviors i.e., jumps and oscillations.
Currently, no available method can extract jumps and AM-FM oscillatory
components directly from the data. In our novel approach, we design and solve
a variational optimization problem to accomplish this task. The optimization
formulation includes a regularization term to minimize the bandwidth of all
signal modes for effective oscillation modeling, and a prior for extracting
the jump component. Our method addresses the limitations of conventional AM-FM
signal decomposition methods in extracting jumps, as well as the limitations
of existing jump extraction methods in decomposing multiscale oscillations. By
employing an optimization framework that accounts for both multiscale
oscillatory components and discontinuities, our methods show superior
performance compared to existing decomposition techniques. We demonstrate the
effectiveness of our approaches on synthetic, real-world, single-channel, and
multivariate data, highlighting their utility in three specific applications:
Earth’s electric field signals, electrocardiograms (ECG), and
electroencephalograms (EEG).

###### Index Terms:

Variational mode decomposition, multivariate data, empirical mode
decomposition, biomedical applications, jump extraction

##  I Introduction

Signal decomposition techniques are crucial for extracting and separating
individual components from complex signals. Traditionally, nonstationary
signal decomposition methods have focused on the extraction of amplitude-
modulated and frequency-modulated (AM-FM) components [1, 2, 3]. The notable
algorithms include Empirical Mode Decomposition (EMD) [1], Variational Mode
Decomposition (VMD) [2], and Synchrosqueezed Transform (SST) [3]. These
techniques decompose non-stationary signals into their constituent intrinsic
mode functions (IMFs) and have found numerous cross-disciplinary applications
[4, 5, 6]. Extensions for multivariate signals have also been developed [7,
8].

Despite their effectiveness, EMD and VMD primarily focus on the extraction of
oscillatory components from the data. However, in many practical applications,
data often displays a combination of oscillations, abrupt discontinuities (or
jumps), and concurrent trends. Fig. 1 (a) and (g), demonstrate two practical
cases where jumps and oscillations are part of the signal and it is of
interest to decompose such signals into their constituent components. In such
signals, the discontinuities in the data significantly degrade the performance
of traditional decomposition algorithms [9]. This is demonstrated in Fig. 1
(b), where the output of the multivariate EMD (MEMD) algorithm mistakenly
models a jump as an oscillation. That happens because, in the frequency
domain, a sudden jump has a wide-band spectrum, affecting all frequency
components. Consequently, filtering methods that use cutoff frequencies to
distinguish jumps struggle to remove such artifacts effectively. These
anomalies can affect all frequency components, making it difficult to separate
them from the meaningful parts of the signal.

To address these challenges, a class of algorithms has been developed to
handle jumps and noise more effectively [10, 11, 12, 13, 14]. Particularly,
[10] presents a method that estimates both discontinuities and a corresponding
piecewise-smooth signal. Another study [11] combines Tikhonov and total
variation regularizations to reconstruct piecewise-smooth signals by
considering the signal as a sum of piecewise constant and smooth components.
Researchers in [12] introduce the $\ell_{1}$ trend filtering method, which is
a variant of Hodrick–Prescott filtering for estimating piecewise linear trends
in time series data. Further advancements include the inverse Potts energy
functionals [13], which use dynamic programming and the alternating direction
method of multipliers (ADMM) to recover jump-sparse and sparse signals from
noisy, incomplete data. This technique overcomes the limitations associated
with the $\ell_{1}$ norm by employing non-convex sparse regularizers (i.e.,
$\ell_{0}$ pseudonorm-based penalty), demonstrating superior performance in
signal reconstruction compared to classical and recent approaches such as TV
minimization [15], orthogonal matching pursuit [16], and iterative hard
thresholding [17].

Figure 1: Two examples of real-world signals with jump components:
multichannel electroencephalography (EEG) signal (a) and the single channel
ECG signal (g). The decomposition of these two real-world signals into their
AM-FM components using MEMD (A) and JOT (B). The MEMD wrongly models step-like
jumps as the oscillatory components as illustrated in row (b). JOT also fails
to decompose the AM-FM modes, as they are all mixed up, as shown in row (j).

The authors of [18] proposed the Jumps, Oscillations, and Trends (JOT) model
to decompose signals into three components. This model was inspired by the
image decomposition work in [19], which proposed decomposing an image into its
geometric part and its oscillatory texture or noisy part. The JOT method
effectively separates signal components without artifacts, making it suitable
for preprocessing. However, JOT has limitations in obtaining AM-FM oscillatory
components since it considers all potential oscillations as a single
component. Consequently, mid-low oscillatory components must be modeled either
as high-frequency oscillations or as trends. Additionally, being a two-stage
method, JOT has many input parameters that are challenging to tune. This
tuning becomes even more complex when JOT is combined with other decomposition
methods, each with their own set of parameters.

Here, we develop and present a novel method that jointly decomposes a
nonstationary signal into jump and oscillatory (AM-FM) modes. Mathematically,
the input signal is assumed to follow this model:

$\displaystyle f(t)=\sum_{k=1}^{K}u_{k}(t)+v(t)+n(t),$ (1)

where $u_{k}(t)$ denotes the $K$ number of oscillatory or AM-FM components,
$v(t)$ represents the jump component in the input data, and $n(t)$ refers to
noise. The main goal of this work is to obtain $u_{k}(t)$ and $v(t)$, given
the input signal $f(t)$. We achieve this by developing an optimization
formulation that includes two different priors: one for obtaining the signal’s
AM-FM modes and another for its discontinuities. One penalty term enforces the
AM-FM modes to be band-limited, while the other minimizes a function of the
derivative of the jump component to precisely compute the amplitude of the
piecewise-constant component. The proposed method is built on a single-stage
algorithm rather than a two-stage process, offering a robust solution with
fewer input parameters and a higher likelihood of convergence. The resulting
optimization problem is solved using the ADMM approach. In the sequel, the
proposed method is termed as Jump Plus AM-FM Mode Decomposition (JMD).

The rest of the paper is organized as follows: Section II reviews the related
works and methods and their limitations in acquiring AM-FM modes and jumps. In
section III we present the steps to build our formulation and solution to the
problem, followed by several experiments and applications of our methods to
both synthetic and real-life signals in section IV. Finally, discussions and
concluding remarks are given in sections V and VI.

##  II Related Works and Limitations

###  II-A AM-FM Mode Decomposition Methods

Mode decomposition methods aim to decompose a signal into a finite number of
amplitude- and frequency-modulated (AM-FM) modes. Let $f(t)$ denote an input
signal, then the AM-FM decomposition methods obtain the modes
$u_{k}(t)=a_{k}(t)$cos$(\phi_{k}(t))$, assuming the following model of the
signal $f(t)$:

$\displaystyle f(t)=\sum_{k=1}^{K}u_{k}(t)$ (2)

While several methods have been proposed to extract AM-FM modes from a
nonstationary signal, the most widely used are EMD [1], VMD [2], and SST [3].
Invariably, all mode decomposition methods have been designed to process data
containing smooth, oscillatory, and noise components only and their
performance suffers significantly when discontinuities or jumps are part of
the input signal. Specifically, the mode decomposition methods attempt to
approximate jumps through oscillatory modes, resulting in spurious components
(see row (b) in Fig. 1).

In their original formulations, AM-FM mode decomposition methods can only
handle univariate or single-channel data sets. However, many real-world
signals contain multiple channels exhibiting complex relationships e.g., the
EEG signal shown in Fig. 1 (a). To effectively process such data, several
multivariate extensions of the AM-FM mode decomposition methods have emerged
[7, 8]. Unlike their univariate counterparts, these methods can model cross-
channel dependencies in multivariate data leading to improved performances in
wide-ranging applications [20, 6, 21]. However, similar to their univariate
counterparts, the multivariate VMD (MVMD) and multivariate EMD (MEMD) are not
able to effectively extract jump components from data. For example, Fig. 1 (b)
- (f) demonstrate the weak performance of MEMD in extracting the jump
component from the multichannel signal (a), where the jump component (b) is
inaccurately modeled as an oscillation.

###  II-B Jump Decomposition Methods

Jump decomposition methods aim to accurately separate jump from signals [10,
11, 12, 13, 14, 15, 16, 17]. These methods are crucial for identifying and
removing anomalies, such as jumps or spikes, which can obscure the meaningful
analysis of the underlying signal trends and frequencies. Notably, Potts
focuses on partitioning the signal into segments with distinct levels,
effectively capturing jumps [13]. This approach uses an $\ell_{0}$ pseudonorm-
based penalty, making it robust against noise and capable of preserving sharp
transitions. More recently, the two-stage approach JOT has been proposed that
aims to decompose signals into three distinct components: Jumps ($v^{*}$),
Oscillations ($n^{*}$), and Trends ($w^{*}$) [18] as follows:

$\displaystyle f(t)=v^{*}(t)+n^{*}(t)+w^{*}(t).$ (3)

The method first identifies these components by solving a non-convex
minimization problem to capture jumps and oscillations, then refines these
components by redistributing trends and oscillatory residuals. Although JOT
and Potts are effective in separating signal components without artifacts,
they struggle with extracting AM-FM oscillatory modes due to their signal
models, which treat all oscillations as a single component. For instance, Fig.
1 (h) - (j) demonstrates the limitation of JOT in extracting the AM-FM modes
from a single-channel ECG signal with jump artifacts (g). In addition to the
mode mixing problem, all oscillatory components are inaccurately modeled as a
single oscillation (as illustrated in row (j)).

##  III Jump Plus AM-FM Mode Decomposition

In this section, we present our proposed variational JMD model to extract AM-
FM modes $u_{k}(t)$ and jump component $v(t)$, from an input signal $f(t)$, in
accordance with the signal model given in (1). Our approach aims to achieve
the following two goals: i) extracting the jumps by incorporating a sparsity-
inducing regularizer, and ii) decomposing the signal into a finite number of
oscillatory components (sub-signals or modes) by minimizing a term related to
the bandwidth of the modes. This combination leads to a minimization problem
for effectively separating the discontinuity features from oscillating parts
of the original signal.

In this framework, the assumption is that each oscillatory mode is
predominantly localized around a central frequency [2] while the
discontinuities are wide-band components on the spectrum. Considering these
two opposite features, incorporating two different terms in our formulation is
required. We will discuss this further in the following subsections.

###  III-A Oscillatory Components

In order to obtain the AM-FM oscillatory components from the data, we use the
following optimization formulation of the VMD algorithm that aims to minimize
the bandwidth of the AM-FM components:

$\displaystyle
J_{1}=\sum_{k}\Big{\|}\partial_{t}\Big{[}\tilde{u}_{k+}(t)e^{-j\omega_{k}t}\Big{]}\Big{\|}^{2}_{2},$
(4)

where the symbol $\partial_{t}$ represents the partial derivative with respect
to time, $u_{k+}(t)=u(t)+j\mathcal{H}u(t)=a(t)e^{j\phi(t)}$ (and
$\tilde{u}_{+}(\omega)=(1+\text{sgn}(\omega))\tilde{u}(\omega)$) refers to the
analytic signal in which $\mathcal{H}$ denotes the Hilbert transform, and
${\omega_{k}}$ is the center frequency of each mode. In essence, the squared
$l^{2}$-norm term estimates the bandwidths of the $k$-th oscillatory modes of
the signal. The minimization of the bandwidth of a mode is tantamount to
approximating it to a sinusoidal signal, which by definition has the minimum
possible bandwidth.

###  III-B Jump Component

To extract the jump component of a signal, we use a reparameterized and
rescaled minimax concave (MC) penalty term, originally proposed in [22]. This
penalty term penalizes the derivative of the jump component to promote a
piecewise-constant component. It aims to accurately preserve the amplitude of
the piecewise-constant signal component while improving the representation of
sharp discontinuities in the jump component:

$\displaystyle J_{2}=\int_{0}^{\infty}\phi(\|\partial_{t}v(t)\|_{2};b)\quad
dt,$ (5)

where $\partial_{t}v(t):=v^{\prime}(t)$ represents the first derivative of the
jump component. The non-convex sparsity-promoting penalty function
$\phi(\cdot;b):[0,+\infty)\to[0,1]$ is defined as a piecewise-quadratic
function [18], as follows:

$\displaystyle\centering\phi(x;b)\@add@centering$
$\displaystyle=\begin{cases}-\frac{b}{2}x^{2}+\sqrt{2b}x&x\in[0,\frac{2}{b}),\\\
1&x\in[\frac{2}{b},+\infty).\end{cases}$ (6)

In the above function, parameter $b$ changes the degree of non-convexity, such
that for $b\xrightarrow{}0$, the MC penalty is defined as $\phi(x;b)=|x|$ and
$\phi(\cdot;b)$ tends to $l_{0}$-pseudonorm for $b\xrightarrow{}\infty$.
Utilizing this term is motivated by an image decomposition study [19], which
proposed decomposing the image into three parts: geometric parts, and
oscillatory texture or noisy part by means of TV denoising [23]. This
formulation has also been used successfully to extract jump components from
input data [18, 14].

###  III-C Optimization Formulation

In our proposed solution to obtain $v(t)$ and $u_{k}(t)$ using input data
$f(t)$ in (1), we combine (4) and (5) within a single optimization formulation
as follows:

$\begin{split}\begin{aligned}
\underset{\\{{\tilde{u}}_{k}\\},\\{\omega_{k}\\},\\{v\\}}{\text{minimize}}\footnotesize\Big{\\{}\alpha{J}_{1}+\beta{J}_{2}\Big{\\}},\end{aligned}\end{split}$
(7) $\small\centering s.t.\quad v+\sum_{k}u_{k}=f.\@add@centering$ (8)

Next, we perform a variable splitting technique by defining an auxiliary
variable $x:=\partial_{t}v$ to solve the problem of non-differentiability of
the term $\phi(\cdot;b)$. Then we can write our problem in this form:

$\begin{split}\begin{aligned}
\underset{\\{{{u}}_{k}\\},\\{\omega_{k}\\},\\{v\\}}{\text{minimize}}\footnotesize\Big{\\{}\alpha\sum_{k}\Big{\|}\partial_{t}\Big{[}{u}_{k+}(t)e^{-j\omega_{k}t}\Big{]}\Big{\|}^{2}_{2}\\\
+\beta\int_{0}^{\infty}\phi(\|\partial_{t}v(t)\|_{2};b)\hskip
2.84526ptdt\Big{\\}},\end{aligned}\end{split}$ (9) $\small\centering
s.t.\begin{cases}v+\sum_{k}u_{k}=f,\\\
x=\partial_{t}v\end{cases}\@add@centering$ (10)

The above optimization formulation is constrained and cumbersome to solve. We
convert it into unconstrained formulation using Lagrangian multipliers,
resulting in the following augmented Lagrangian function:

$\displaystyle\mathcal{L}\\{{{u}}_{k},\omega_{k},v,x,\rho,\lambda\\}=\alpha\sum_{k}\Big{\|}\partial_{t}\Big{[}{u}_{k+}(t)e^{-j\omega_{k}t}\Big{]}\Big{\|}^{2}_{2}+\Big{\|}f(t)-$
(11)
$\displaystyle\Big{(}v(t)+\sum_{k}u_{k}(t)\Big{)}\Big{\|}^{2}_{2}+\Big{\langle}\lambda(t),f(t)-\Big{(}v(t)+\sum_{k}u_{k}(t)\Big{)}\Big{\rangle}$
$\displaystyle+\beta\int_{0}^{\infty}\phi(\|x(t)\|_{2};b)dt-\Big{\langle}\rho(t),x(t)-\partial_{t}v(t)\Big{\rangle}+\frac{\gamma}{2}\Big{\|}x(t)-$
$\displaystyle\partial_{t}v(t)\Big{\|}^{2}_{2},$

where $\gamma$ is the penalty scalar parameter, and $\lambda(t)$ and $\rho(t)$
are the vectors of Lagrangian multipliers associated with constraints
$\sum_{k}u_{k}+v=f$ and $x=\partial_{t}v$, respectively.

###  III-D Optimization Solution

The above formulation is solved using the alternate direction method of
multipliers (ADMM). To accomplish that, we divide the augmented Lagrangian
formulation into several sub-optimization problems and solve those separately
and iteratively.

####  III-D1 Minimization w.r.t $u_{k}$, $\omega_{k}$, and $\lambda$

The sub-optimization problem to solve for $u_{k}$, $\omega_{k}$, and $\lambda$
can be obtained from (11) by including only those terms that depend on them:

$\displaystyle\mathcal{L}\\{{{u}}_{k},\omega_{k},\lambda\\}=\alpha\sum_{k}\Big{\|}\partial_{t}\Big{[}{u}_{k+}(t)e^{-j\omega_{k}t}\Big{]}\Big{\|}^{2}_{2}+\Big{\|}f(t)-$
(12)
$\displaystyle\Big{(}v(t)+\sum_{k}u_{k}(t)\Big{)}+\frac{\lambda(t)}{2}\Big{\|}^{2}_{2}.$

As proposed in [2], by solving (12) in the Fourier domain, taking the partial
derivative, and some algebraic manipulations, the associated update equations
will be obtained as follows:

$\displaystyle\hat{u}_{k}^{(i+1)}(\omega)=\frac{\hat{f}(\omega)-\sum_{l\neq
k}u_{l}^{(i)}(\omega)-\hat{v}^{(i)}(\omega)+\frac{\hat{\lambda}^{(i)}(\omega)}{2}}{1+2\alpha(\omega-\omega_{k}^{(i)})^{2}},$
(13)

$\small\omega_{k}^{i+1}=\frac{\int^{\infty}_{0}\omega\Big{|}\hat{u}_{k}^{(i+1)}(\omega)\Big{|}^{2}d\omega}{\int^{\infty}_{0}\Big{|}\hat{u}_{k}^{(i+1)}(\omega)\Big{|}^{2}d\omega},$
(14)

$\displaystyle\lambda^{(i+1)}=\lambda^{(i)}(\omega)+\tau_{1}\Bigg{(}{\hat{f}}(\omega)-\Bigg{[}{\hat{v}}^{(i)}(\omega)+\sum_{k}u_{k}^{(i+1)}(\omega)\Big{)}\Bigg{]}\Bigg{)}.$
(15)

####  III-D2 Minimization w.r.t $v$, $x$, and $\rho$

The suboptimization problem to solve for $v$, $x$, and $\rho$ can be obtained
from (11) and is given as follows:

$\displaystyle\mathcal{L}\\{v,x,\rho\\}=\beta\int_{0}^{\infty}\phi(\|x(t)\|_{2};b)dt-\Big{\langle}\rho(t),x(t)-\partial_{t}v(t)\Big{\rangle}$
(16)
$\displaystyle+\frac{\gamma}{2}\Big{\|}x(t)-\partial_{t}v(t)\Big{\|}^{2}_{2}+\Big{\|}f(t)-\Big{(}v(t)+\sum_{k}u_{k}(t)\Big{)}\Big{\|}^{2}_{2}$
$\displaystyle+\Big{\langle}\lambda(t),f(t)-\Big{(}v(t)+\sum_{k}u_{k}(t)\Big{)}\Big{\rangle}.$

This sub-problem is cumbersome to solve directly in the continuous domain
owing to the first term involving the infinite integral. Fortunately, the term
involving $\phi$ can be converted to a discrete form and has been shown to be
strongly convex under necessary conditions in [24]. In the following, we
describe the steps to convert (16) in the discrete form and solve it to obtain
the update equations for $v$, $x$, and $\rho$.

Assume that we deal with a finite signal111Although we make this assumption to
simplify the sub-function including $\phi(\cdot;b)$ and proceed in the
discrete domain; practically this is also numerically established for (12) as
explained in [2].. Then, we can rewrite the optimization function (11) to
minimize w.r.t $v$, $x$ and $\rho$ as follows:

$\displaystyle\mathcal{L}\\{v,x,\rho\\}=\beta\sum_{j=1}^{N}\phi(\|x_{j}\|_{2};b)-\Big{\langle}\rho,x-Dv\Big{\rangle}+\frac{\gamma}{2}\Big{\|}x-Dv\Big{\|}^{2}_{2}$
(17)
$\displaystyle+\Big{\|}f-\Big{(}v+\sum_{k}u_{k}\Big{)}\Big{\|}^{2}_{2}+\Big{\langle}\lambda,f-\Big{(}v+\sum_{k}u_{k}\Big{)}\Big{\rangle},$

where $D$ is defined as the first-order derivative matrix as follows:

$\displaystyle D=\begin{bmatrix}-1&1\\\ &-1&1\\\ &&\ddots&\ddots\\\
&&&-1&1\end{bmatrix}\in\mathbb{R}^{(N-1)\times N}.$ (18)

Next, by solving (17), taking the partial derivative with respect to $v$, and
zero-initializing the vectors $x^{(0)}$ and $\rho^{(0)}$, we can obtain the
$i$-th iteration of $v$ as follows:

$\displaystyle v^{(i+1)}=(\gamma D^{T}D+2I)^{-1}\Big{(}D^{T}\rho^{(i)}+\gamma
D^{T}x^{(i)}+$ (19) $\displaystyle
2(f-\sum_{k}u_{k}^{(i+1)})-\lambda^{(i+1)}\Big{)}.$

Moreover, we can write the subproblem for obtaining $x$ by omitting the
constant terms as:

$\displaystyle
x^{(i+1)}=\arg\min_{x\in\mathbb{R}^{N}}\Big{\\{}\beta\sum_{n=1}^{N}\phi(|x_{j}|;b)-\langle\rho,x\rangle+\frac{\gamma}{2}\lVert
x-Dv\rVert_{2}^{2}\Big{\\}}.$ (20)

The above optimization problem, when expressed in terms of its individual
components, can be seen as comprising $N$ separate one-dimensional problems:

$\displaystyle
x_{j}^{(i+1)}=\arg\min_{x\in\mathbb{R}}\Big{\\{}\mu\phi(|x|;b)+\frac{1}{2}\lVert
x-(Dv^{(i)})_{j}+\frac{\rho^{(i)}_{j}}{\gamma}\rVert_{2}^{2}\Big{\\}},$ (21)
$\displaystyle j=1,\ldots,N,$

where $\mu=\frac{\beta}{\gamma}$. As demonstrated in [24], the problems in
(21) are strongly convex if and only if sufficient conditions for strong
convexity are met:

$\displaystyle
b<\frac{1}{\mu}\Rightarrow\gamma>b\beta\Rightarrow\gamma=\tau_{2}b\beta,\text{
for }\tau_{2}\in\mathbb{R},\tau_{2}>1.$ (22)

Given assumption (22), closed-form solutions for the problems stated in (21)
can be derived as:

$\displaystyle x_{j}^{(i+1)}=\min\left(\max\left(\frac{1}{1-\mu
b}-\frac{\frac{\mu\sqrt{2b}}{1-b\mu}}{\lvert
h_{j}\rvert},0\right),1\right)h_{j},$ (23)

where $h_{j}=(Dv^{(i)})_{j}+\frac{\rho^{(i)}_{j}}{\gamma}$. Finally, the
update for $\rho$ will be governed by the following equation:

$\displaystyle\rho^{(i+1)}=\rho^{(i)}-\gamma\left(x^{(i+1)}-Dv^{(i+1)}\right)$
(24)

With all relevant update equations at our disposal, we list the steps for
obtaining the AM-FM modes $u_{k}(t)$ and jump component $v(t)$, given an input
signal $f(t)$, in Algorithm 1.

Algorithm 1 JMD

Input and user-defined parameters: $f,K,\alpha,\beta,\overline{b}$
Initialize: $\epsilon\leftarrow
10^{-7},\hat{u}_{k}^{(1)},\omega_{k}^{(1)},x^{(1)},\lambda^{(1)},\rho^{(1)}\leftarrow
0$
Generate: Discrete operator $D$ in (18)

1:$b=2/{\overline{b}^{2}},\quad\gamma=\tau_{2}b\beta$

2:repeat

3: $i\leftarrow i+1$

4: for $k=1:K$ do

5: Update mode $\hat{u}_{k}$:

6: $\hat{u}_{k}^{(i+1)}(\omega)\leftarrow\text{Using Equation \eqref{uk}}$

7: Update $\omega_{k}$:

8: $\omega_{k}^{(i+1)}\leftarrow\text{Using Equation \eqref{omega}}$

9: end for

10: $\sum_{k}u_{k}\xleftarrow{IFT}\sum_{k}\hat{u}_{k}(\omega)$ (to compute
$v$)

11: Update $\lambda,v,x,$ and $\rho$:

12: $\lambda^{(i+1)}(\omega)\leftarrow\text{Using Equation \eqref{lambda}}$

13: $v^{(i+1)}\leftarrow\text{Using Equation \eqref{vjump}}$

14: $x^{(i+1)}\leftarrow\text{Using Equation \eqref{xupdate}}$

15: $\rho^{(i+1)}\leftarrow\text{Using Equation \eqref{rho}}$

16: $r^{(i+1)}=\sum_{k}{u}_{k}^{(i+1)}+{v}^{(i+1)}$

17: $\hat{v}^{(i+1)}(\omega)\xleftarrow{FT}v^{(i+1)}$ (to update
$\hat{u}_{k}(\omega)$)

18:until Convergence:
$\|r^{(i+1)}-r^{(i)}\|_{2}^{2}/\|r^{(i)}\|_{2}^{2}<\epsilon$

###  III-E Multivariate Extension

Here we present the multivariate extension of the JMD method. We start with
the multivariate extension of the jump and AM-FM mode decomposition model
given in (1). The model is given as follows:

$\mathbb{f}(t)=\sum_{k=1}^{K}\mathbf{u}_{k}(t)+\mathbb{v}(t)+\mathbb{n}(t),$
(25)

where $\mathbb{f}(t)$ denotes a multivariate signal consisting of C number of
data channels $\mathbb{f}(t)=[f_{1}(t),...f_{C}(t)]$. The signal can also be
represented by a matrix $\mathbb{F}$ in the discrete
$\mathbf{F}\in\mathbb{R}^{C\times N}$. Further, $\mathbb{u}_{k}(t)$ are the
set of AM-FM decomposed modes, all of which are multivariate with $C$ number
of data channels. Finally, $\mathbb{v}(t)$ and $\mathbb{n}(t)$ refer to the
multivariate jump component in the data and noise, respectively.

The goal of the multivariate extension of the JMD method is to extract the
jump components $\mathbb{v}(t)$ in all channels and the AM-FM modes
$\mathbb{u}_{k}(t)$ from input signal $\mathbb{f}(t)$. Unlike the univariate
algorithm, the AM-FM modes will be aligned in frequency in the multivariate
algorithm, enabling the direct comparison between the modes from different
channels [8]. To achieve this, we utilize the cost function of the MVMD method
in our problem such that the bandwidth of $\mathbb{u}_{k}(t)$ can be
determined by calculating the $\ell_{2}$ norm of the gradient of the
harmonically shifted $\mathbb{u}_{k}^{+}(t)$ [8]:

$\displaystyle
J_{1}=\alpha{\sum_{k}\sum_{c}}\Bigg{\|}\partial_{t}\Big{[}{u}_{k,c}^{+}(t)e^{-j\omega_{k}t}\Big{]}\Bigg{\|}^{2}_{2},$
(26)

where $c=1,2,\ldots,C$ denote number of channels. Said that, the term for
obtaining jumps will not be affected; thus, it remains almost the same:

$\displaystyle
J_{2}=\beta\int_{0}^{\infty}\phi(\|\partial_{t}v_{c}(t)\|_{2};b)\hskip
2.84526ptdt.$ (27)

Utilizing the same procedure of our univariate JMD optimization formulation
i.e. in (9) and (10), we can formulate the problem as follows:

$\begin{split}\begin{aligned}
\underset{\\{u_{k,c}\\},\\{\omega_{k}\\},\\{v_{c}\\}}{\text{minimize}}\footnotesize\Big{\\{}\alpha\sum_{k}\sum_{c}\Big{\|}\partial_{t}\Big{[}{u}_{k,c}^{+}(t)e^{-j\omega_{k}t}\Big{]}\Big{\|}^{2}_{2}\\\
+\beta\int_{0}^{\infty}\phi(\|\partial_{t}v_{c}(t)\|_{2};b)\hskip
2.84526ptdt\Big{\\}}\end{aligned}\end{split}$ (28)

$\small\centering s.t.\begin{cases}v_{c}+\sum_{k}u_{k,c}=f_{c}\\\
x_{c}=\partial_{t}v_{c}\end{cases}.\@add@centering$ (29)

Similarly, we convert the above optimization formulation into an unconstrained
formulation using Lagrangian multipliers, resulting in the following augmented
Lagrangian function:

$\displaystyle\mathcal{L}\\{{{u}}_{k,c},\omega_{k},v_{c},x_{c},\rho_{c},\lambda_{c}\\}=\alpha\sum_{k}\sum_{c}\Big{\|}\partial_{t}\Big{[}{u}_{k,c}^{+}(t)e^{-j\omega_{k}t}\Big{]}\Big{\|}^{2}_{2}$
(30)
$\displaystyle+\sum_{c}\Big{\|}f_{c}(t)-\Big{(}v_{c}(t)+\sum_{k}u_{k,c}(t)\Big{)}\Big{\|}^{2}_{2}+\sum_{c}\Big{\langle}\lambda_{c}(t),f_{c}(t)-$
$\displaystyle\Big{(}v_{c}(t)+\sum_{k}u_{k,c}(t)\Big{)}\Big{\rangle}+\beta\int_{0}^{\infty}\phi(\|x_{c}(t)\|_{2};b)dt$
$\displaystyle-\Big{\langle}\rho_{c}(t),x_{c}(t)-\partial_{t}v_{c}(t)\Big{\rangle}+\frac{\gamma}{2}\Big{\|}x_{c}(t)-\partial_{t}v_{c}(t)\Big{\|}^{2}_{2}.$

The above Lagrangian function can be similarly solved using ADMM as our
presented procedure for the univariate JMD method in previous subsections
(i.e., in equations (12) to (24)). Therefore, the relevant update equations
for $u_{k,c}$, $\omega_{k}$, $\lambda_{c}$, $v_{c}$, $x_{c}$, and $\rho_{c}$
will be given as follows:

$\displaystyle\hat{u}_{k,c}^{(i+1)}(\omega)=\frac{\hat{f_{c}}(\omega)-\sum_{l\neq
k}u_{l,c}^{(i)}(\omega)-\hat{v_{c}}^{(i)}(\omega)+\frac{\hat{\lambda_{c}}^{(i)}(\omega)}{2}}{1+2\alpha(\omega-\omega_{k}^{(i)})^{2}},$
(31)

$\small\omega_{k}^{i+1}=\frac{{\sum_{c}}\int^{\infty}_{0}\omega\Big{|}\hat{u}_{k,c}^{(i+1)}(\omega)\Big{|}^{2}d\omega}{{\sum_{c}}\int^{\infty}_{0}\Big{|}\hat{u}_{k,c}^{(i+1)}(\omega)\Big{|}^{2}d\omega},$
(32)

$\displaystyle\lambda_{c}^{(i+1)}=\lambda_{c}^{(i)}(\omega)+\tau_{1}\Bigg{(}{\hat{f_{c}}}(\omega)-\Bigg{[}{\hat{v_{c}}}^{(i)}(\omega)+\sum_{k}u_{k,c}^{(i+1)}(\omega)\Big{)}\Bigg{]}\Bigg{)},$
(33) $\displaystyle v_{c}^{(i+1)}=(\gamma
D^{T}D+2I)^{-1}\Big{(}D^{T}\rho_{c}^{(i)}+\gamma D^{T}x_{c}^{(i)}+$ (34)
$\displaystyle 2(f_{c}-\sum_{k}u_{k,c}^{(i+1)})-\lambda_{c}^{(i+1)}\Big{)},$
$\displaystyle x_{c,(j)}^{(i+1)}=\min\left(\max\left(\frac{1}{1-\mu
b}-\frac{\frac{\mu\sqrt{2b}}{1-b\mu}}{\lvert
h_{c,(j)}\rvert},0\right),1\right)h_{c,(j)},$ (35)

and

$\displaystyle\rho_{c}^{(i+1)}=\rho_{c}^{(i)}-\gamma\left(x_{c}^{(i+1)}-Dv_{c}^{(i+1)}\right).$
(36)

Having all these update equations, the multivariate JMD algorithm can be
summarized in Algorithm 2.

Algorithm 2 MJMD

Input and user-defined parameters: $f,K,\alpha,\beta,\overline{b}$
Initialize: $\epsilon\leftarrow
10^{-7},\hat{u}_{k,c}^{(1)},\omega_{k}^{(1)},x_{c}^{(1)},\lambda_{c}^{(1)},\rho_{c}^{(1)}\leftarrow
0$
Generate: Discrete operator $D$ in (18)

1:$b=2/{\overline{b}^{2}},\quad\gamma=\tau_{2}b\beta$

2:repeat

3: $i\leftarrow i+1$

4: for $k=1:K$ do

5: for $c=1:C$ do

6: Update mode $\hat{u}_{k,c}$:

7: $\hat{u}_{k,c}^{(i+1)}(\omega)\leftarrow\text{Using Equation \eqref{muk}}$

8: end for

9: Update $\omega_{k}$:

10: $\omega_{k}^{(i+1)}\leftarrow\text{Using Equation \eqref{momega}}$

11: end for

12: for $k=1:K$ do

13: for $c=1:C$ do

14: $u_{k,c}\xleftarrow{IFT}\hat{u}_{k,c}(\omega)$ (to compute $v_{c}$)

15: end for

16: end for

17: Update $\lambda_{c},v_{c},x_{c},$ and $\rho_{c}$:

18: for $c=1:C$ do

19: $\lambda_{c}^{(i+1)}(\omega)\leftarrow\text{Using Equation
\eqref{mlambda}}$

20: $v_{c}^{(i+1)}\leftarrow\text{Using Equation \eqref{mvjump}}$

21: $x_{c}^{(i+1)}\leftarrow\text{Using Equation \eqref{mxupdate}}$

22: $\rho_{c}^{(i+1)}\leftarrow\text{Using Equation \eqref{mrho}}$

23: $\hat{v}_{c}^{(i+1)}(\omega)\xleftarrow{FT}v_{c}^{(i+1)}$ (to update
$\hat{u}_{k,c}(\omega)$)

24: end for

25: for $c=1:C$ do

26: $r_{c}^{(i+1)}=\sum_{k}{u}_{k,c}^{(i+1)}+{v_{c}}^{(i+1)}$

27: end for

28:until Convergence:
$\sum_{c}\|r_{c}^{(i+1)}-r_{c}^{(i)}\|_{2}^{2}/\|r_{c}^{(i)}\|_{2}^{2}<\epsilon$

###  III-F Input Parameters

Input parameters play a key role in the convergence of our presented methods.
The parameters $K$ and $\alpha$ must be precisely determined in advance. A
large $K$ can lead to duplicate modes, whereas a small $K$ causes mode mixing
problems. Selecting a proper $K$ must be done by analyzing the signal in hand
a priori, and adaptive strategies have also been used to tune $K$ [25, 26,
27]. Similarly, a high value of $\alpha$ can result in incorrect noisy modes
or cause convergence issues, while a low value of $\alpha$ can lead to mode
mixing problems. Essentially, a high $\alpha$ value may be suitable when
dealing with a dense spectrum with closely spaced frequency components. We
found that setting the parameter $\alpha$ in the range of $10^{3}-10^{5}$
leads to satisfactory performance in most experiments. Additionally, different
center frequency initialization strategies can cause methods like our
presented JMD methods to diverge from the actual center frequencies in the
spectrum [2]. The parameter $\tau_{1}$ (in (15)) balances faithfully
reconstructing and denoising the input signal. Opting for a smaller $\tau_{1}$
value may be advisable when addressing significant noise levels in the input
signals. This parameter is set in the range of $0-1$. Parameter $\beta$ is a
regularization parameter that affects the extracted jump component. This
parameter should be set approximately to the inverse of the number of expected
jumps in the signal i.e., $1/\text{(number of expected jumps)}$. This ensures
that the regularization term properly balances the contribution of the jump
component without over-penalizing or under-penalizing it. Fine-tuning,
however, is then needed to achieve the best performance for the specific
applications. Parameter $\tau_{2}>1$ is typically set in the range of 1.1 to
50, with a default value of 10, for initializing $\gamma$. Parameter
$\overline{b}$ refers to the expected minimal jump height, with a default
value of 0.3. It indicates the minimal height above which every jump
$\left|(Dv)_{j}\right|$ greater than $\overline{b}$ is equally penalized. The
user should set $\overline{b}$ based on the expected minimal jump height
observed in the input data. If every non-zero derivative
$\left|(Dv)_{j}\right|>0$ satisfies $\left|(Dv)_{j}\right|>\overline{b}$, the
value of the $v$-regularization term will be equal to the number of non-zero
derivatives (jumps) in $v$.

##  IV Experiments and Results

We employ the proposed algorithm on a series of synthetic and real-world
signals to evaluate the effectiveness of our approach. Time is discretized
into $t\in[0,1]$ in all simulations. We also compare the components obtained
from the JMD and MJMD methods with the results obtained from the VMD [2], MVMD
[8], and the JOT [18] models. Note that in this section, the amplitude of
components (y axis) is scaled between $0-1$ in all figures for a better
illustration of the results. Additionally, the presented methods are evaluated
by setting the parameter $\tau_{1}=0$ as we assume that the input signal is
corrupted by noise. For the comparative methods, we set the parameter values
to ensure a fair comparison with the proposed methods. For instance, in VMD
and MVMD, we use the same values for $K$ as those used in the proposed (M)JMD
method. In other cases, we use the parameter values for the comparative
methods that produce the best results.

Figure 2: Synthetic signal (a), and decomposed components obtained from the
JMD method (rows (b) - (e)), VMD (rows (f) - (i)), and JOT (rows (j) - (l)).
In all plots, the red line refers to the extracted components and the black
dashed line indicates the original components of the signal.

###  IV-A Synthetic Signal

####  IV-A1 Example 1 (Univariate Signal)

To test the presented univariate model we start with a simple single-channel
signal $f_{1}(t)$ that is composed of three 4Hz, 80Hz, and 200 Hz oscillatory
components, one jump $v(t)$, and additive white Gaussian noise (AWGN) $\eta$
with zero mean and variances of $\sigma=0.1$ as follows:

$\displaystyle f_{1}(t)$ $\displaystyle=\cos(2\pi 4t)+\cos(2\pi 80t)+\cos(2\pi
200t)+v(t)+\eta.$ (37)

The input parameters for this experiment were set to $K=3$, $\alpha=5000$,
$\beta=0.03$, $\mkern 1.5mu\overline{\mkern-1.5mub\mkern-1.5mu}\mkern
1.5mu=0.45$, and $\tau_{2}=50$.

Figure 2 visualizes the extracted components from the proposed JMD method
including three oscillatory modes and one jump. It is evident that our method
was able to extract all components precisely. We also show the results of the
comparative methods in this figure. For instance, VMD clearly fails to extract
jumps as its formulation has not been designed to capture jumps. As the jump
component has a broad spectrum, a broad spectral overlap is expected between
the jump and the other three components. Due to this spectral overlap, VMD
struggles to capture the pure modes as it performs like a Wiener filtering of
the residual on the spectrum [2]. Although the JOT method is able to capture
jumps, it cannot obtain oscillatory modes accurately. Especially, the 80Hz and
200Hz extracted components got mixed up in one component illustrated in Fig.
2\- row (l). Because JOT decomposes the signal into three components including
jump, oscillation, and trend, it only provides three components even in cases
where there are more components than three.

####  IV-A2 Example 2 (Multivariate Signal)

To assess the performance of the proposed multivariate extension of the JMD
method, we performed a series of experiments and simulations on various
multivariate signals. The comprehensive results of these analyses are detailed
in this section. We will focus on the method’s ability to align common
frequency scales across multiple data channels, which is crucial for various
scientific and engineering applications involving multivariate data [28, 29,
6]. Here we showcase the mode-alignment property of MJMD using synthetic test
signals composed of multiple oscillatory components across different channels
and AWGN. We also compare these results with those obtained from MVMD and JOT.
In the case of JOT, the results are obtained by applying this method to each
channel separately. The synthetic multivariate data includes 3 channels (i.e.,
$C_{1},C_{2}$, and $C_{3}$) as follows:

$\displaystyle C_{1}(t)$ $\displaystyle=10\cos(2\pi 1t)+2\cos(2\pi
24t)+v(t)+\eta$ (38) $\displaystyle C_{2}(t)$ $\displaystyle=4\cos(2\pi
1t)+2\cos(2\pi 24t)+2\cos(2\pi 48t)$ $\displaystyle+2\cos(2\pi 128t)+\eta$
$\displaystyle C_{3}(t)$ $\displaystyle=10\cos(2\pi 1t)+2\cos(2\pi
48t)+2\cos(2\pi 128t)+v(t)+\eta$

The input parameters for this experiment were set to $K=3$, $\alpha=5000$,
$\beta=0.05$, $\mkern 1.5mu\overline{\mkern-1.5mub\mkern-1.5mu}\mkern
1.5mu=0.45$, and $\tau_{2}=50$.

Looking at the above formulation, it can be seen that the jump component only
exists in $C_{1}$, and $C_{3}$. The result of applying our multivariate JMD
scheme is illustrated in Fig. 3 (column A, rows (b) - (e)). This figure
visualizes the extracted jump and three multivariate oscillatory modes. Our
method precisely extracted all components. The comparative results obtained
from MVMD are also illustrated in column B. It is clear that MVMD fails to
extract jumps as it utilizes a similar basis to VMD which has not been
designed to capture jumps. As illustrated in Fig. 3 (rows (g) - (j)), the pure
modes are obtained only from channel 2 ($C_{2}$) where the jump component does
not exist. This clearly illustrates that MVMD fails to obtain meaningful
components when jump components are present in the signal. Further, channel-
wise JOT only precisely obtains the signal components of channel $C_{1}$ where
the channel consists of a jump, a trend, and an oscillation, as shown in Fig.
3 column C. It fails to purely extract all AM-FM oscillatory components for
the other channels, resulting in the third oscillatory mode (128 Hz) being
absent in channels $C_{2}$ and $C_{3}$. Moreover, tuning the input parameters
for each channel is cumbersome in channel-wise JOT.

Figure 3: Synthetic multivariate signal (rows (a) and (f)), and decomposed
components obtained from MJMD method (column A, rows (b) - (e)), MVMD (column
B, rows (g) - (j)), and channel-wise JOT (column C, rows (l) - (o)). The red
line refers to the extracted components and the black dashed line indicates
the original components of the signal in all plots.

###  IV-B Real world signal:

####  IV-B1 ECG-derived respiration

ECG-derived respiration (EDR) has demonstrated extensive utility in the
diagnosis and management of multiple conditions such as stress and sleep
disorders in recent decades [30]. Fundamentally, respiration affects the ECG
signal by adding a low-frequency component to it. Especially the presence of
jumps in ECG, which often manifest as strong harmonics in the low-frequency
range (0.3-0.7 Hz), complicates obtaining EDR. The AM-FM signal decomposition
methods fail in such scenarios. To assess the efficacy of our method under
such conditions, we employ ECG data derived from an initially clean dataset
sourced from the MIMIC Database, contaminated with simulated electrode motion
artifacts [31]. The simulated motion artifact consists of a jump component
characterized by a single rise and a single fall. The signal-to-noise ratio
(SNR) measured in dB after adding the artifact to the signal is -6 $dB$. The
decomposition results obtained from the JMD method (A), VMD (B), and JOT (C)
are depicted in Fig. 4. Notably, we decomposed the input signal into 12 modes
using both our method and VMD; however, we only illustrated the jump
components ((b) and (d)) and EDR ((c) and (e)) among them. In the case of JOT,
we similarly highlighted the components pertinent to EDR and jump. The
reference respiratory signal simultaneously recorded with ECG is plotted in
rows (c), (e), and (g) to compare the accuracy of all methods. Fig. 4 (A)
highlights the effectiveness of our proposed approach in obtaining EDR despite
the presence of jumps. Comparative results from VMD and JOT indicate that both
methods struggle to isolate pure EDR, as the weak low-frequency respiratory
signal component is obscured by other harmonics of the wide band and a strong
jump in the frequency domain. JOT cannot model AM-FM components, and VMD does
not account for the jump component. Obtaining EDR in such a scenario requires
a combination of AM-FM signal-jump modeling (i.e., the same formulation as our
method) and setting a considerably high $\alpha$.

The input parameters for this experiment were set to $K=13$, $\alpha=2\times
10^{5}$, $\beta=0.2$, $\mkern 1.5mu\overline{\mkern-1.5mub\mkern-1.5mu}\mkern
1.5mu=0.12$, and $\tau_{2}=3.6$.

Figure 4: Decomposition of the input ECG signal (record number “055m” from
MIMIC database) with additive simulated jump component (a) into a jump (b),
(d), and (f) and EDR (c), (e), and (g) using our method (Column A), VMD
(Column B), and JOT (Column C). The red line refers to the extracted
components and the black dashed line indicates the original components of the
signal in all plots. In EDR plots (c), (e), and (g), the black dashed line
refers to the reference respiratory signal simultaneously recorded with ECG
available in the MIMIC dataset. Figure 5: Decomposition of the EEG signal
(rows (a) and (g)) into a jump (rows (b) and (h)), trend (rows (c) and (i)),
alpha-band (rows (e) and (k)), and beta-band (rows (f) and (l)) using our
multivariate JMD (MJMD) scheme (column A) and MVMD (column B). Channel-wise
JOT (column C) only provides the jump, trend, and a combination of all
oscillations. Figure 6: Decomposition of the EMF signal (rows (a), (h), and
(o)) into jumps, oscillatory components, and trends using our JMD scheme
(column A), VMD (column B), and JOT (column C). The last plots of each column
(in red) show the reconstructed signals by removing the extracted jump. VMD
fails to properly extract the jump with sharp edges, resulting in a spike that
destroys the signal quality in the middle where the jump occurred. It also
mixes the trend and the jump into a single component (as shown in row (i)).

####  IV-B2 EEG sub-band time series

EEG is a multichannel recording of the brain’s electrical activity as measured
from the scalp. Typically, each channel consists of a time series of
electrical signals produced by the brain. EEG signals are known to exhibit
oscillatory behavior in specific frequency bands, such as theta (4-8 Hz),
alpha (8-12 Hz), beta (13-30 Hz), and gamma (30-150 Hz) bands. Each of these
frequency bands is associated with different underlying neural processes.
Thus, an essential task in many EEG applications is to break down these neural
oscillations into distinct frequency bands. However, the presence of jumps can
significantly impact the decomposition, potentially leading to erroneous and
unpure oscillations due to the presence of numerous harmonics across all
frequency ranges. This experiment aims to showcase the effectiveness of our
proposed approach in identifying and characterizing jumps and inherent
frequency modes in motor imagery electroencephalogram (EEG) data. The
resulting frequency bands time series hold utility for various applications,
including connectivity and power spectrum analysis.

We utilize a publicly available 16-channel EEG dataset collected during motor
imagery tasks [32], and the outcomes (i.e., the results of 4 channels as a few
examples) are illustrated in Fig. 5. Data inherently includes huge jumps that
extremely affect the results of the conventional AM-FM signal decomposition
methods. This example vividly demonstrates our method’s capability to capture
oscillatory modes despite significant jump interference. Comparative analysis
with MVMD and channel-wise JOT highlights their respective limitations in this
context: 1) MVMD struggles to isolate jumps and oscillatory modes. The modes
obtained from MVMD particularly include jump harmonics across all scales. For
instance, the alpha-band modes especially in channels $C_{2}$ and $C_{3}$
(Fig. 5, column B, row (k)) show significantly more energy than those obtained
from our MJMD method (Fig. 5, column A, row (e)). A similar observation is
nearly made for the beta band components (row (l)). This increased energy is
likely due to the jump harmonics at this scale, which can lead to incorrect
assessments, such as in alpha-band connectivity and power analysis.
Additionally, as shown in Fig. 5 (column B, row (h)), MVMD mistakenly
interprets jumps as oscillations. 2) Channel-wise JOT also fails in obtaining
AM-FM multiple frequency bands oscillations. Results obtained from channel-
wise JOT (Fig. 5, column C), aside from dealing with the mode mixing problem,
demonstrate faulty extracted jumps and trends in channels $C_{3}$ and $C_{4}$.
One may first use JOT to extract jumps from each channel and then apply the
MVMD to extract oscillatory modes associated with each band. This integration
requires considerable processing time.

The input parameters for this experiment were set to $K=15$, $\alpha=2000$,
$\beta=0.035$, $\mkern 1.5mu\overline{\mkern-1.5mub\mkern-1.5mu}\mkern
1.5mu=0.32$, and $\tau_{2}=3$.

####  IV-B3 Earth’s Magnetic and Electric Field (EMF)

Earth’s magnetic and electric fields are actively studied in geophysics for
various applications, including atmospheric electricity [33] and earthquake
prediction [34]. These signals contain both AM-FM oscillations and abrupt
changes. Geomagnetic jerks, which appear as step-like changes, are crucial for
interpreting the secular variation and secular acceleration of the magnetic
field, thereby providing insights into core dynamics and mantle conductivity
[35]. The oscillatory components are influenced by the solar cycle and its
harmonics, annual and semiannual variations, and internal secular variation
[36].

The proposed JMD method effectively extracts both AM-FM modes and step-like
abrupt changes, offering a comprehensive analysis of the Earth’s magnetic
field’s temporal behavior. AM-FM mode decomposition methods alone fail to
capture meaningful jumps with sharp edges, while methods like JOT only provide
pure jump and trend components without decomposing the AM-FM oscillations. Our
subsequent analysis will demonstrate the capability of the JMD method to
decompose all AF-FM components and jump discontinuities from the signal.
However, a detailed analysis of the extracted oscillatory components and other
transient phenomena [37] is beyond the scope of this paper.

In Fig. 6 first row, we show a sample of the Earth’s electric field measured
by the CSES-01 satellite [38]. In the second row, we present the results
obtained from our method. We also compare our results with the results
obtained from JOT and VMD. The proposed method better separates the jumps
producing cleaner oscillatory components. On the contrary, JOT cannot
decompose the multiscale oscillatory modes and VMD is not capable of
separating the jump component from the oscillations. Although the JOT
framework shows potential to be used for the pre-processing and extraction of
trends and jumps which are only two components of this kind of data, our
method can be used for the full decomposition of the data into oscillations
and jumps in this scenario. This can potentially provide better information
about the behavior of the data since the oscillatory components are also
available.

The input parameters for this experiment were set to $K=6$, $\alpha=10^{5}$,
$\beta=0.05$, $\mkern 1.5mu\overline{\mkern-1.5mub\mkern-1.5mu}\mkern
1.5mu=0.04$, and $\tau_{2}=8$.

##  V Discussion

A few comments on the performance and comparison between our JMD and other
approaches are in order. First, our (M)JMD are variational approaches capable
of decomposing the signal into both oscillatory components and jumps, whereas
no other common mode decomposition techniques, including (M)VMD, are not
designed to capture jumps. Although discontinuities in the signal can be
decomposed using jump extraction methods such as JOT, unlike our JMD, these
algorithms do not provide multiple oscillatory components. For instance, JOT
decomposes the signals into three parts: jumps, oscillations (including all
frequency components), and trends. This means one needs to utilize other
decomposition techniques alongside JOT if the extraction of AM-FM modes is
also required. In such situations, the decomposition of a signal into jumps
and oscillatory modes must be performed in multiple steps, as JOT is naturally
a two-step method. Each step in JOT includes an optimization process with its
own set of input parameters, which can be challenging to tune. This complexity
arises from using individual balancing parameters for each mathematical term
in the JOT formulation, concerning each component of the signal. The situation
becomes even more cumbersome when additional decomposition techniques such as
VMD are required along with JOT to obtain AM-FM modes of the signal. For
example, JOT and VMD have five and two key input parameters, respectively.
Tuning a total of seven input parameters can be very time-consuming, making
our JMD a better solution with only four key input parameters.

Second, the optimization formulations of our (M)JMD are non-convex and
therefore not guaranteed to converge. However, there is strong evidence
suggesting that the optimization formulations in (M)JMD will converge, as they
are based on well-established, non-convex heuristics similar to those used in
(M)VMD. Our extensive experiments have shown that our methods converge in all
the cases we have considered. Given that (M)VMD has demonstrated empirical
convergence properties for a wide range of signals, we expect the (M)JMD
algorithms to perform similarly.

Third, regarding the sensitivity of our methods to input parameters, it is
important to mention that balancing $K$ and $\alpha$ can be challenging and
can significantly affect the decomposition results in our methods. A large $K$
may lead to duplicate modes, while a small $K$ can cause mode mixing problems.
Similarly, a high value of $\alpha$ can lead to incorrect noisy modes or cause
convergence issues, whereas a low value of $\alpha$ can also result in mode
mixing. Essentially, a high value of $\alpha$ is suitable when dealing with a
dense spectrum containing frequency components that are very close to each
other. Although other input parameters must be set in all these techniques,
the methods are not highly sensitive to changes in those parameters. For
better convergence of our approaches, we suggest the following suitable ranges
for input parameters: (i) $\alpha=10^{3}-10^{5}$, (ii) $0<\beta\leq 1$, (iii)
$\gamma=0-50$, (iv) $\overline{b}=0.1-0.9$.

##  VI Conclusion

We have presented novel fully data-driven methods for decomposing multivariate
signals into jumps and their intrinsic mode functions at multiple frequency
scales. By combining tools and methods from traditional AM-FM signal
decomposition and jump extraction concepts, we have introduced JMD and MJMD
algorithms for obtaining piecewise constant functions (as jumps) and intrinsic
oscillatory components of univariate and multivariate signals, respectively.

These methods enable users to model sharp discontinuities in the signal as
jumps and separate them from the oscillatory modes. They are based on a
variational optimization formulation solved using the alternating direction
method of multipliers (ADMM). To our knowledge, no other conventional mode
decomposition methods define the signal as composed of wide-band jumps and
narrow-band AM-FM signals, incorporating a specific prior to account for jumps
in the signal.

Moreover, JMD competes well with other specific jump decomposition methods
like JOT, which cannot obtain multiple oscillatory modes of the signals.
Additionally, the JOT method is more challenging to initialize compared to our
JMD approach, as it requires more input parameters to be predefined due to its
use of different priors to account for oscillations and trends separately. In
contrast, JMD is more likely to converge since it models trends and
oscillations using only one prior.

The proposed methods have been tested on various synthetic and real data,
including real-life electroencephalogram (EEG), electrocardiogram (ECG), and
Earth’s electric field signals, demonstrating the broad applicability of these
approaches across different fields.

While the initial promise of the proposed methods has been demonstrated,
future work could focus on: i) improving the cost function of the JMD method
to better model discontinuities in either the continuous or discrete domain,
and ii) eliminating the need for presetting parameter $K$ using a successive
approach proposed in [39].

The MATLAB implementation codes for our JMD and MJMD methods are available at
the following link:
https://se.mathworks.com/matlabcentral/profile/authors/16833607.

## Acknowledgment

Funding: This work was supported by IRRAS USA Inc., San Diego, CA 92130 US.

## References

  * [1] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, S. H. H., Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis,” _Proc. R. Soc. A: Math. Phys. Eng. Sci._ , vol. 454, no. 1971, pp. 903–995, 1998.
  * [2] K. Dragomiretskiy and D. Zosso, “Variational mode decomposition,” _IEEE Trans. Signal Process._ , vol. 62, no. 3, pp. 531–544, 2014.
  * [3] I. Daubechies, J. Lu, and H.-T. Wu, “Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool,” _Appl. Comput. Harmon. Anal._ , vol. 30, no. 2, pp. 243–261, 2011.
  * [4] B. Pang, M. Nazari, and G. Tang, “Recursive variational mode extraction and its application in rolling bearing fault diagnosis,” _Mech. Syst. Signal Process._ , vol. 165, p. 108321, 2022.
  * [5] H. Akbari, M. T. Sadiq, S. Siuly, Y. Li, and P. Wen, “Identification of normal and depression eeg signals in variational mode decomposition domain,” _Health Inf. Sci. Syst._ , vol. 10, no. 1, p. 24, 2022.
  * [6] P. Gupta, K. K. Sharma, and S. D. Joshi, “Baseline wander removal of electrocardiogram signals using multivariate empirical mode decomposition,” _Healthc. Technol. Lett._ , vol. 2, no. 6, pp. 164–166, 2015.
  * [7] N. Rehman and D. Mandic, “Multivariate empirical mode decomposition,” _Proc. R. Soc. A: Math. Phys. Eng. Sci._ , vol. 466, no. 2117, pp. 1291–1302, 2010.
  * [8] N. Rehman and H. Aftab, “Multivariate variational mode decomposition,” _IEEE Trans. Signal Process._ , vol. 67, no. 23, pp. 6039–6052, 2019.
  * [9] A. Stallone, A. Cicone, and M. Materassi, “New insights and best practices for the successful use of empirical mode decomposition, iterative filtering and derived algorithms,” _Sci. Rep._ , vol. 10, no. 1, p. 15161, 2020.
  * [10] M. Storath, L. Kiefer, and A. Weinmann, “Smoothing for signals with discontinuities using higher order mumford–shah models,” _Numer. Math._ , vol. 143, pp. 423–460, 2019.
  * [11] A. Gholami and S. M. Hosseini, “A balanced combination of tikhonov and total variation regularizations for reconstruction of piecewise-smooth signals,” _Signal Process._ , vol. 93, no. 7, pp. 1945–1960, 2013.
  * [12] S.-J. Kim, K. Koh, S. Boyd, and D. Gorinevsky, “$\ell_{1}$ trend filtering,” _SIAM Rev._ , vol. 51, no. 2, pp. 339–360, 2009.
  * [13] M. Storath, A. Weinmann, and L. Demaret, “Jump-sparse and sparse recovery using potts functionals,” _IEEE Trans. Signal Process._ , vol. 62, no. 14, pp. 3654–3666, 2014.
  * [14] I. Selesnick, A. Lanza, S. Morigi, and F. Sgallari, “Non-convex total variation regularization for convex denoising of signals,” _J. Math. Imaging Vis._ , vol. 62, no. 6, pp. 825–841, 2020.
  * [15] D. Needell and R. Ward, “Stable image reconstruction using total variation minimization,” _SIAM J. Imaging Sci._ , vol. 6, no. 2, pp. 1035–1058, 2013\.
  * [16] J. A. Tropp, “Greed is good: Algorithmic results for sparse approximation,” _IEEE Trans. Inf. Theory_ , vol. 50, no. 10, pp. 2231–2242, 2004.
  * [17] T. Blumensath and M. E. Davies, “Iterative hard thresholding for compressed sensing,” _Appl. Comput. Harmon. Anal._ , vol. 27, no. 3, pp. 265–274, 2009\.
  * [18] A. Cicone, M. Huska, S.-H. Kang, and S. Morigi, “Jot: a variational signal decomposition into jump, oscillation and trend,” _IEEE Trans. Signal Process._ , vol. 70, pp. 772–784, 2022.
  * [19] Y. Meyer, _Oscillating patterns in image processing and nonlinear evolution equations: the fifteenth Dean Jacqueline B. Lewis memorial lectures_. AMS, 2001, vol. 22.
  * [20] N. ur Rehman, Y. Xia, and D. P. Mandic, “Application of multivariate empirical mode decomposition for seizure detection in eeg signals,” in _Proc. IEEE Eng. Med. Biol. Conf. 2010_. IEEE, 2010, pp. 1650–1653.
  * [21] Y. Lv, R. Yuan, and G. Song, “Multivariate empirical mode decomposition and its application to fault diagnosis of rolling bearing,” _Mech. Syst. Signal Process._ , vol. 81, pp. 219–234, 2016.
  * [22] C.-H. Zhang, “Nearly unbiased variable selection under minimax concave penalty,” 2010.
  * [23] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” _Physica D: nonlinear phenomena_ , vol. 60, no. 1-4, pp. 259–268, 1992.
  * [24] M. Huska, A. Lanza, S. Morigi, and I. Selesnick, “A convex-nonconvex variational method for the additive decomposition of functions on surfaces,” _Inverse Probl._ , vol. 35, no. 12, p. 124008, 2019.
  * [25] X. Zhang, Q. Miao, H. Zhang, and L. Wang, “A parameter-adaptive vmd method based on grasshopper optimization algorithm to analyze vibration signals from rotating machinery,” _Mech. Syst. Signal Process._ , vol. 108, pp. 58–72, 2018.
  * [26] H. Li, T. Liu, X. Wu, and Q. Chen, “An optimized vmd method and its applications in bearing fault diagnosis,” _Meas._ , vol. 166, p. 108185, 2020\.
  * [27] J. Lian, Z. Liu, H. Wang, and X. Dong, “Adaptive variational mode decomposition method for signal processing based on mode characteristic,” _Mech. Syst. Signal Process._ , vol. 107, pp. 53–77, 2018.
  * [28] N. u. Rehman, D. Looney, T. Rutkowski, and D. Mandic, “Bivariate emd-based image fusion,” in _IEEE Workshop Statistical Signal Process. 2009_. IEEE, 2009, pp. 57–60.
  * [29] H. Hao, H. Wang, and N. Rehman, “A joint framework for multivariate signal denoising using multivariate empirical mode decomposition,” _Signal Process._ , vol. 135, pp. 263–273, 2017.
  * [30] M. Nazari and S. M. Sakhaei, “Variational mode extraction: A new efficient method to derive respiratory signals from ecg,” _IEEE J. Biomed. Health Inform._ , vol. 22, pp. 1059–1067, 2017.
  * [31] A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, and H. E. Stanley, “Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals,” _Circulation_ , vol. 101, no. 23, pp. e215–e220, 2000.
  * [32] H. Cho, M. Ahn, S. Ahn, M. Kwon, and S. Jun, “Eeg datasets for motor imagery brain–computer interface,” _GigaSci._ , vol. 6, no. 7, p. gix034, 2017\.
  * [33] Y. Lu, Z. Zhou, S. Gu, D. Wu, J. Guo, and H. Tao, “Research on lightning warning method based on the characteristics of atmospheric electric field,” in _Proc. 33rd Int. Conf. Lightning Prot. (ICLP) 2016_. IEEE, 2016, pp. 1–5.
  * [34] D. Yang and K. Yang, “Multi-step prediction of strong earthquake ground motions and seismic responses of sdof systems based on emd-elm method,” _Soil Dyn. Earthq. Eng._ , vol. 85, pp. 117–129, 2016.
  * [35] M. Mandea, R. Holme, M. A. Pais, K. Pinheiro, A. Jackson, and G. Verbanac, “Geomagnetic jerks: rapid core field variations and core dynamics,” _Space Sci. Rev._ , vol. 155, pp. 147–175, 2010.
  * [36] V. Courtillot and J.-L. Le Mouel, “Time variations of the earth’s magnetic field-from daily to secular,” in _Annu. Rev. Earth Planet. Sci._ Annual Reviews Palo Alto, CA, 1988, pp. 389–476.
  * [37] P. Picozza, R. Battiston, G. Ambrosi, S. Bartocci, L. Basara, W. J. Burger, D. Campana, L. Carfora, M. Casolino, G. Castellini, _et al._ , “Scientific goals and in-orbit performance of the high-energy particle detector on board the cses,” _ApJS_ , vol. 243, no. 1, p. 16, 2019.
  * [38] P. Diego, J. Huang, M. Piersanti, D. Badoni, Z. Zeren, R. Yan, G. Rebustini, R. Ammendola, M. Candidi, Y. Guan, _et al._ , “The electric field detector on board the china seismo electromagnetic satellite—in-orbit results and validation,” _Earth Planets Space._ , vol. 72, no. 1, pp. 1–16, 2020.
  * [39] M. Nazari and S. M. Sakhaei, “Successive variational mode decomposition,” _Signal Process._ , vol. 174, p. 107610, 2020.