# Variational Mode Decomposition (VMD) using R

11 Jun 2022This blog post is a slight modification of the R package Vignette. 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/).

## Example Use Cases

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.

### Noise Removal (1-Dimensional VMD)

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)
alpha = 2000 # moderate bandwidth constraint
tau = 0 # noise-tolerance (no strict fidelity enforcement)
K = 9 # 9 modes
DC = FALSE # no DC part imposed
init = 1 # initialize omegas uniformly
tol = 1e-6
vec_arrhythmia = arrhythmia[['MLII']]
set.seed(1)
arr_vmd = vmd(data = vec_arrhythmia,
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 ...
op <- par(mfrow = c(3, K/3))
for (item in K:1) {
item_mode = glue::glue("IMF {K-item+1}")
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,

```
imfs_noise_free = rowSums(arr_vmd$u[, -c(1,K)])
op <- par(mfrow = c(1,1))
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)
```

### Synthetic image segmentation (2-Dimensional VMD)

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)
pth_texture = system.file('Matlab', 'VMD_2D', 'texture.mat', package = 'VMDecomp')
data = R.matlab::readMat(pth_texture)
data = data$f
dim(data)
# [1] 256 256
alpha = 1000 # bandwidth constraint
tau = 0.25 # Lagrangian multipliers dual ascent time step
K = 5 # number of modes
DC = TRUE # includes DC part (first mode at DC)
init = 1 # initialize omegas randomly, may need multiple runs!
tol = 1e-7 # tolerance (for convergence)
set.seed(2)
res_2d = vmd(data = data,
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.
end_dims = dim(res_2d$u)[3]
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,

```
res_2d_lst = lapply(1:end_dims, function(x) res_2d$u[,, x])
res_2d_lst = append(res_2d_lst, list(data), after = 0)
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 ...
init_plt = OpenImageR::GaborFeatureExtract$new()
mlt_plt = init_plt$plot_multi_images(list_images = res_2d_lst,
par_ROWS = 3,
par_COLS = 2,
axes = TRUE,
titles = c("Input Image", glue::glue("Imf_{1:end_dims}")))
```

### Package Installation & Citation:

To install the package from CRAN use,

```
install.packages("VMDecomp")
```

and to download the latest version of the package from Github,

```
remotes::install_github('mlampros/VMDecomp')
```

If you use the **VMDecomp** R package in your paper or research please
cite both **VMDecomp** and the **original articles / software**
`https://cran.r-project.org/web/packages/VMDecomp/citation.html`

:

```
@Manual{,
title = {VMDecomp: Variational Mode Decomposition using R},
author = {Lampros Mouselimis},
year = {2022},
note = {R package version 1.0.0},
url = {https://CRAN.R-project.org/package=VMDecomp},
}
```