The **VMDecomp** R package is an *RcppArmadillo* implementation of the Matlab Code for *Variational Mode Decomposition* (1- and 2-dimensional) based on the papers,

- “Variational Mode Decomposition” by K. Dragomiretskiy and D. Zosso in IEEE Transactions on Signal Processing, vol. 62, no. 3, pp. 531-544, Feb.1, 2014, https://doi.org/10.1109/TSP.2013.2288675
- “Two-Dimensional Variational Mode Decomposition” by Dragomiretskiy, K., Zosso, D. (2015), In: Tai, XC., Bae, E., Chan, T.F., Lysaker, M. (eds) Energy Minimization Methods in Computer Vision and Pattern Recognition. EMMCVPR 2015. Lecture Notes in Computer Science, vol 8932. Springer, https://dx.doi.org/10.1007/978-3-319-14612-6_15

The Matlab code is available to download in the Author’s website (https://math.montana.edu/dzosso/code/).

VMD has been used in many scientific areas with true or synthetic data. I’ll explain the functionality of the *VMDecomp* R package using two simple use cases.

In the Authors’ paper “Variational Mode Decomposition” (page 13 / 16) is mentioned, *“A ‘clean’ electrocardiogram (ECG) signal can be reconstructed by summing all but the first and last VMD modes, thus discarding the low-frequency baseline oscillation and most of the high-frequency noise.”* I’ll use the *‘arrhythmia’* data (especially the first column which corresponds to the *modified limb lead II (MLII)*) for this purpose. This same example appears also online (Matlab Mathworks).

```
require(VMDecomp)
require(glue)
data(arrhythmia)
= 2000 # moderate bandwidth constraint
alpha = 0 # noise-tolerance (no strict fidelity enforcement)
tau = 9 # 9 modes
K = FALSE # no DC part imposed
DC = 1 # initialize omegas uniformly
init = 1e-6
tol
= arrhythmia[['MLII']]
vec_arrhythmia
set.seed(1)
= vmd(data = vec_arrhythmia,
arr_vmd alpha = alpha,
tau = tau,
K = K,
DC = DC,
init = init,
tol = tol,
verbose = TRUE)
# 1-dimensional VMD starts ...
# --------------------------------
# The 1-dimensional VMD starts ...
# --------------------------------
# Iteration: 10 uDiff: 2.28024
# Iteration: 20 uDiff: 1.05777
# Iteration: 30 uDiff: 0.515047
# Iteration: 40 uDiff: 0.346114
# Iteration: 50 uDiff: 0.277273
# Iteration: 60 uDiff: 0.240478
# Iteration: 70 uDiff: 0.212126
# Iteration: 80 uDiff: 0.20821
# Iteration: 90 uDiff: 0.196179
# ...............
# Iteration: 470 uDiff: 0.00728391
# Iteration: 480 uDiff: 0.00698528
# Iteration: 490 uDiff: 0.00663912
# Iteration: 500 uDiff: 0.00629761
# -----------------------------------------
# The algorithm converged in iteration: 500
# -----------------------------------------
# Elapsed time: 0 hours and 0 minutes and 19 seconds.
str(arr_vmd)
# List of 3
# $ u : num [1:10000, 1:9] -0.0695 -0.0697 -0.0702 -0.0707 -0.0711 ...
# $ u_hat: cplx [1:10000, 1:9] 6.54e-04+0.0i -3.83e-03+5.1e-04i -1.99e-03+8.5e-04i ...
# $ omega: num [1:500, 1:9] 0 0.00561 0.00595 0.00596 0.00586 ...
<- par(mfrow = c(3, K/3))
op
for (item in K:1) {
= glue::glue("IMF {K-item+1}")
item_mode plot(x = arr_vmd$u[, item], type = 'l', main = item_mode, xlab = 'Time', ylab = '')
}
```

The previous plot shows the intrinsic mode functions (IMF’s). Based on the Authors’ statement if we exclude the first and last Modes (IMF’s) we’ll receive a noise-free ECG,

```
= rowSums(arr_vmd$u[, -c(1,K)])
imfs_noise_free
<- par(mfrow = c(1,1))
op
plot(x = vec_arrhythmia, type = 'l', col = "blue", xlab = 'Time', ylab = 'Signal')
lines(x = imfs_noise_free, col = "orange")
legend("topright",
legend = c("Original ECG", "Clean ECG"),
col = c("blue", "orange"),
lty = 1:2,
cex = 0.8)
```

To explain the functionality of the *“vmd()”* R function for 2-dimensional data I’ll load the synthetic image used in section 3.1 of the Authors’ second paper where it’s mentioned,

*“The synthetic image is a composition of spatially overlapping basic shapes, more precisely six ellipses and a rectangle, with frequency patterns varying in both periodicity and direction… The spectrum is ideal for segmentation due to modes being deliberately both well isolated and narrow-banded. The resolution of the synthetic image is 256 x 256 and the experiment was run with parameters Alpha = 1000 and K = 5. This experiment converged in 520 iterations which took 45 seconds on a standard PC. The algorithm has no problems in accuracy nor timeliness in segmenting the image into its five constituent sub-images. The first is the DC component of the image - a solid ellipse and rectangle, while the four remaining decompositions show clear separation of the patterned ellipses. Due to the solid pieces having sharp edges, their spectra are not band limited, and only smoothed versions are recovered.”*

```
require(R.matlab)
require(OpenImageR)
= system.file('Matlab', 'VMD_2D', 'texture.mat', package = 'VMDecomp')
pth_texture = R.matlab::readMat(pth_texture)
data = data$f
data dim(data)
# [1] 256 256
= 1000 # bandwidth constraint
alpha = 0.25 # Lagrangian multipliers dual ascent time step
tau = 5 # number of modes
K = TRUE # includes DC part (first mode at DC)
DC = 1 # initialize omegas randomly, may need multiple runs!
init = 1e-7 # tolerance (for convergence)
tol
set.seed(2)
= vmd(data = data,
res_2d alpha = alpha,
tau = tau,
K = K,
DC = DC,
init = init,
tol = tol,
verbose = TRUE)
# --------------------------------
# The 2-dimensional VMD starts ...
# --------------------------------
# Iteration: 10 uDiff: 42.7451 omegaDiff: 0.00012395
# Iteration: 20 uDiff: 1.05264 omegaDiff: 1.15039e-09
# Iteration: 30 uDiff: 0.31892 omegaDiff: 4.14469e-10
# Iteration: 40 uDiff: 0.117648 omegaDiff: 1.45234e-10
# Iteration: 50 uDiff: 0.0491776 omegaDiff: 6.399e-11
# Iteration: 60 uDiff: 0.0222704 omegaDiff: 3.33701e-11
# Iteration: 70 uDiff: 0.0106591 omegaDiff: 1.87647e-11
# Iteration: 80 uDiff: 0.0053165 omegaDiff: 1.08343e-11
# Iteration: 90 uDiff: 0.00273909 omegaDiff: 6.30121e-12
# ............
# Iteration: 250 uDiff: 6.31694e-07 omegaDiff: 1.22605e-15
# Iteration: 260 uDiff: 4.33964e-07 omegaDiff: 8.05652e-16
# Iteration: 270 uDiff: 3.08055e-07 omegaDiff: 5.66624e-16
# Iteration: 280 uDiff: 2.26559e-07 omegaDiff: 4.31241e-16
# Iteration: 290 uDiff: 1.728e-07 omegaDiff: 3.54727e-16
# Iteration: 300 uDiff: 1.36554e-07 omegaDiff: 3.11406e-16
# Iteration: 310 uDiff: 1.11508e-07 omegaDiff: 2.86634e-16
# -----------------------------------------
# The algorithm converged in iteration: 317
# -----------------------------------------
# Elapsed time: 0 hours and 0 minutes and 41 seconds.
= dim(res_2d$u)[3]
end_dims str(res_2d)
# List of 3
# $ u : num [1:256, 1:256, 1:5] 3.16e-05 1.80e-05 -3.32e-06 -3.17e-05 -6.32e-05 ...
# $ u_hat: cplx [1:256, 1:256, 1:5] 7.65e-04+0.00i -6.63e-04-3.28e-04i 7.74e-04-3.15e-04i ...
# $ omega: num [1:317, 1:2, 1:5] 0 0 0 0 0 0 0 0 0 0 ...
```

To the output IMF’s I’ll also include the input image and then I’ll visualize all images in a multi-plot to verify that each IMF represents a different shape of the input image,

```
= lapply(1:end_dims, function(x) res_2d$u[,, x])
res_2d_lst = append(res_2d_lst, list(data), after = 0)
res_2d_lst str(res_2d_lst)
# List of 6
# $ : num [1:256, 1:256] 0 0 0 0 0 0 0 0 0 0 ...
# $ : num [1:256, 1:256] 3.16e-05 1.80e-05 -3.32e-06 -3.17e-05 -6.32e-05 ...
# $ : num [1:256, 1:256] -4.36e-05 -3.08e-05 -1.19e-05 2.28e-05 6.01e-05 ...
# $ : num [1:256, 1:256] -2.67e-07 2.53e-06 1.03e-06 1.94e-06 2.44e-07 ...
# $ : num [1:256, 1:256] 3.31e-06 1.21e-06 2.15e-06 1.66e-07 4.97e-07 ...
# $ : num [1:256, 1:256] 6.72e-06 9.27e-06 1.14e-05 8.79e-06 1.45e-06 ...
= OpenImageR::GaborFeatureExtract$new()
init_plt = init_plt$plot_multi_images(list_images = res_2d_lst,
mlt_plt par_ROWS = 3,
par_COLS = 2,
axes = TRUE,
titles = c("Input Image", glue::glue("Imf_{1:end_dims}")))
```