Introduction to MatchingPursuit package

1 R version and generation time

## R version:    4.5.1
## Generated on: 19-sierpień-2026

2 Loading the package

library(MatchingPursuit)

3 Installing the required external software

The below empi_install() function downloads Enhanced Matching Pursuit Implementation external program (or EMPI for short), see Różański (2024), and stores it in the cache directory. The function downloads the EMPI program in a version compatible with the operating system used (Windows, Linux, MacOS-arm64).

First, the user can see where EMPI can be downloaded from.

empi_locate()
## $url
## [1] "https://github.com/develancer/empi/releases/download/1.0.4/empi-1.0.4-windows-x64.zip"
## 
## $fname
## [1] "empi-1.0.4-windows-x64.zip"

The code in the below chunk has been commented out because CRAN’s verification rules prohibit automatic binary downloads. Therefore, this function cannot be executed while generating this vignette.

# empi_install()

User can check whether the EMPI program is installed; if not, an error message is displayed indicating that installation is required.

empi_check()
## [1] "C:\\Users\\Artur\\AppData\\Local/R/cache/R/MatchingPursuit/empi/empi.exe"

4 The purpose and functionality of the package

The presented package enables the analysis of time-series signals using the Matching Pursuit (MP) and Orthogonal Matching Pursuit (OMP) algorithms (see Mallat and Zhang (1993), Pati, Rezaiifar, and Krishnaprasad (1993), Durka (2007), Elad (2010)). Additionally, it supports working with EEG (electroencephalogram) and ECG (electrocardiograph) signals. The functionality of the package includes:

1. Execution of the Matching Pursuit algorithm
   • optimized C++ backend, Enhanced Matching Pursuit Implementation (EMPI)
     
   • reference R implementation (MP)
     
2. Execution of the Orthogonal Matching Pursuit (OMP) algorithm
   
3. Support for Gabor dictionaries via XML-based atom definitions
   
4. Support for biomedical signal formats
   • EDF and EDF+ file support
     
   • WFDB record support
     
   • Standard EEG montages(bipolar, referential, average-reference)
     
5. Visualization of time–frequency maps based on MP, OMP, and EMPI results
   
6. Signal pre-processing using notch, low-pass, high-pass, band-pass and band-stop filters

5 Naming convention

In this vignette, we use the following naming convention:

• The terms MP-R and OMP-R with a suffix R refer to the corresponding backends rather than to algorithms. The terms MP and OMP without a suffix R refer to the corresponding algorithms rather than to a specific implementation.
• EMPI refers to the external backend implementing the Matching Pursuit algorithm, developed by Piotr T. Różański and available on GitHub.

Back to menu

6 Quick start

fs <- 1024
T <- 1
t <- seq(0, T - 1 / fs, 1 / fs)
N <- length(t)

# 7 non-stationary signals.
x1 <- sin(2 * pi * (10 + 40 * t) * t)                            # linear chirp
x2 <- sin(2 * pi * (20 * t^2) * t)                               # nonlinear chirp
x3 <- (1 + 0.5 * sin(2 * pi * 2 * t))  *  sin(2 * pi * 30 * t)   # AM
x4 <- sin(2 * pi * 50 * t + 5 * sin(2 * pi * 3 * t))             # FM
x5 <- exp(-2 * t)  *  sin(2 * pi * 60 * t)                       # decreasing amplitude
x6 <- sin(2 * pi * (5 + 20 * sin(2 * pi * t)) * t)               # frequency modulated sine wave
x7 <- t * sin(2 * pi * 40 * t)                                   # increasing amplitude

signal <- data.frame(x = x1 + x2 + x3 + x4 + x5 + x6 + x7)

# The sample1.csv file contains exactly the same data as shown in Step 1.
file <- system.file("extdata", "sample1.csv", package = "MatchingPursuit")

# The first line of the file contains two values:
# the sampling rate in Hz (1024 Hz here) and the signal duration
# in seconds (1 s here).
out <- read.csv(file, header = FALSE)
head(out)
##           V1
## 1     1024 1
## 2 0.00000000
## 3 1.02492083
## 4 1.93756420
## 5 2.63875099
## 6 3.05949332

signal <- read_csv_signals(file)
str(signal)
## List of 3
##  $ signal            :'data.frame':  1024 obs. of  1 variable:
##   ..$ v1: num [1:1024] 0 1.02 1.94 2.64 3.06 ...
##  $ sampling_frequency: num 1024
##  $ time              : num [1:1024] 0 0.000977 0.001953 0.00293 0.003906 ...
##  - attr(*, "class")= chr "sig"

# file <- system.file("extdata", "sample1.csv", package = "MatchingPursuit")
# signal <- read_csv_signals(file)

# empi_class <- empi_execute (
#   signal = signal,
#   empi_options = "-o local --gabor -i 25",
#   write_to_file = TRUE,
#   path = NULL,
#   file_name = "sample1.db"
# )

# Reading a SQLite file in which all generated atom parameters are stored.
file <- system.file("extdata", "sample1.db", package = "MatchingPursuit")

# Create time-frequency map based on atoms.
out <- tf_map(
  x = file,
  channel = 1,
  mode = "sqrt",
  freq_divide = 4,
  increase_factor = 4,
  display_crosses = FALSE,
  display_atom_numbers = TRUE,
  out_mode = "plot"
)
## Channel number: 1
## Total channels: 1
## Number of atoms: 25
## Sampling frequency: 1024 Hz
## Epoch size (in points): 1024
## Signal length (in seconds): 1
## 
## Energy of the original signal:      2746.16
## Energy of the reconstructed signal: 2582.88
## reconstruction / original %:        94.05

file <- system.file("extdata", "sample1.csv", package = "MatchingPursuit")
signal <- read_csv_signals(file, col_names = "ch1")

# empi_class <- empi_execute(signal = signal)

# plot(empi_class)

7 Electroencephalogram (EEG) analysis

In this section, we demonstrate how to analyze electroencephalography (EEG) signals using the Matching Pursuit algorithm. The package provides a dedicated function for reading files in EDF and EDF+ (European Data Format). It also supports three types of EEG montages and allows for signal filtering.

file <- system.file("extdata", "EEG.edf", package = "MatchingPursuit")

# Read signal parameters and display them in a tabular form.
eeg <- read_edf_signals(file)
signal_eeg <- eeg$signal
sampling_frequency <- eeg$sampling_frequency
time <- seq(0, nrow(signal_eeg) - 1) / sampling_frequency

# Filter parameters that will be used (quite typical in filtering EEG signals).
fc <- design_filters(
   sampling_frequency = sampling_frequency,
   notch = c(49, 51),
   lowpass = 40,
   highpass = 1,
)

# Filtering input signals.
signal_eeg_f <- signal_eeg

for (m in 1:ncol(signal_eeg_f)) {
  signal_eeg_f[, m] = signal::filtfilt(fc$notch, signal_eeg[, m])      # 50Hz notch filter
  signal_eeg_f[, m] = signal::filtfilt(fc$lowpass, signal_eeg_f[, m])  # Low pass IIR Butterworth
  signal_eeg_f[, m] = signal::filtfilt(fc$highpass, signal_eeg_f[, m]) # High pass IIR Butterwoth
}

signal_eeg_f_r <- resample_signal(signal = signal_eeg_f, p = 1, q = 4)
time_64 <- seq(0, nrow(signal_eeg_f_r) - 1) / (sampling_frequency / 4)
sampling_frequency_r <- 64

# Pairs of signals for bipolar montage (so called "double banana").
pairs <- list(
  c("Fp2", "F4"), c("F4", "C4"), c("C4", "P4"), c("P4", "O2"), c("Fp1", "F3"), c("F3", "C3"),  
  c("C3", "P3"), c("P3", "O1"), c("Fp2", "F8"), c("F8", "T4"), c("T4", "T6"), c("T6", "O2"),
  c("Fp1", "F7"), c("F7", "T3"), c("T3", "T5"), c("T5", "O1"), c("Fz", "Cz"), c("Cz", "Pz")
)

# Make the bipolar montage.
signal_eeg_f_r_m <- eeg_montage(
  signal_eeg_f_r,
  montage_type = c("bipolar"),
  bipolar_pairs = pairs
)

# Original signal (first 6 rows, first 6 channels).
signal_eeg_f_r[1:6, 1:6]
##           Fp1         Fp2          F3          F4           F7           F8
## 1  0.01377308  0.01425343  0.01186948  0.01412645  0.006233639  0.009748129
## 2 -1.09734194 -1.25577649 -1.03276959 -1.22017066 -0.548599014 -0.948150122
## 3 -6.73256051 -9.18488110 -6.72161000 -7.98515438 -2.771428391 -6.997555308
## 4 -5.24066880 -7.80164906 -5.09483104 -6.53407891 -1.282696156 -5.476552630
## 5 -4.99119418 -7.31485652 -4.47433883 -6.18484001 -0.891544395 -4.937985774
## 6 -4.50858820 -6.63299757 -3.60798202 -5.48054191 -0.498481654 -4.119310078

# Signal after banana montage (first 6 rows, first 6 channels).
signal_eeg_f_r_m[1:6, 1:6]
##             Fp2_F4       F4_C4        C4_P4         P4_O2       Fp1_F3         F3_C3
## [1,]  0.0001269753  0.01357661  0.005355918 -0.0002423519  0.001903599  0.0004772461
## [2,] -0.0356058328 -0.78857509 -0.353528714  0.0268757282 -0.064572356 -0.0811772988
## [3,] -1.1997267231 -2.88036501 -2.089476968  0.2729456302 -0.010950512 -0.8066486927
## [4,] -1.2675701509 -3.18196004 -2.197927437  0.2499240153 -0.145837757 -0.7927429590
## [5,] -1.1300165091 -3.44430422 -2.306539993  0.1891060065 -0.516855348 -0.7023765071
## [6,] -1.1524556568 -3.69005644 -2.352597831  0.1263696300 -0.900606184 -0.6442897802

# The empi_options parameter is NULL, so the EMPI program is 
# run with the parameters "-o local --gabor -i 50"

# sig <- as_sig(signal_eeg_f_r_m, sampling_frequency_r)

# empi_class <- empi_execute (
#   signal = sig,
#   empi_options = NULL,
#   write_to_file = TRUE,
#   path = NULL,
#   file_name = "EEG_filter_resample_montage.db"
# )

# Reading a SQLite file where all the generated atom's parameters are stored.
file <- system.file("extdata", "EEG_filter_resample_montage.db", package = "MatchingPursuit")

# Generate time-frequency map based on atoms.
out <- tf_map(
  x = file,
  channel = 2,
  mode = "sqrt",
  increase_factor= 8,
  display_crosses = TRUE,
  display_atom_numbers = FALSE,
  out_mode = "plot"
)
## Channel number: 2
## Total channels: 18
## Number of atoms: 50
## Sampling frequency: 64 Hz
## Epoch size (in points): 640
## Signal length (in seconds): 10
## 
## Energy of the original signal:      74035.04
## Energy of the reconstructed signal: 71586.87
## reconstruction / original %:        96.69

edf_processed <- structure(
  list(
    signal = as.data.frame(signal_eeg_f_r_m),
    sampling_frequency = sampling_frequency_r,
    time = time_64,
    signal_names = colnames(signal_eeg_f_r_m),
    record_name = basename(file)
  ),
  class = "edf"
)

plot(
  x = edf_processed,
  begin = 0,
  end = 10,
  panel_height = NULL,
  rainbow = FALSE,
  bg_colour = "white",
  txt_col = "blue",
  zero_line = TRUE,
  main = "EEG after filtering, resampling, and double banana montage"
)

## Actual value of 'panel_height' parameter is: 74

plot(
  x = eeg,
  begin = 0,
  end = 10,
  panel_height = NULL,
  rainbow = FALSE,
  bg_colour = "white",
  txt_col = "blue",
  zero_line = TRUE,
  main = "Original EEG before preprocessing"
)

## Actual value of 'panel_height' parameter is: 75.5

8 Electrocardiography (ECG) analysis

In this section, we demonstrate how to analyze electrocardiography (ECG) signals using the Matching Pursuit algorithm. The package provides a dedicated function for reading files in WFDB (WaveForm DataBase) format. Once the ECG data has been loaded, further analysis is essentially the same as demonstrated in previous chapters.

file <- system.file("extdata", "00001_lr.hea", package = "MatchingPursuit")

out_ecg <- read_wfdb_signals(file)

head(out_ecg$signals)
## NULL
out_ecg$sampling_frequency
## [1] 100
out_ecg$lead_names
##  [1] "I"   "II"  "III" "AVR" "AVL" "AVF" "V1"  "V2"  "V3"  "V4"  "V5"  "V6"
out_ecg$record_name
## [1] "00001_lr"

# Create a list compatible with the empi_execute() function.
signal <- as_sig(
  signal = data.frame(out_ecg$signal),
  sampling_frequency = out_ecg$sampling_frequency
)

str(signal)
## List of 3
##  $ signal            :'data.frame':  1000 obs. of  12 variables:
##   ..$ I  : num [1:1000] -0.119 -0.116 -0.12 -0.117 -0.103 -0.097 -0.119 -0.096 -0.048 -0.037 ...
##   ..$ II : num [1:1000] -0.055 -0.051 -0.044 -0.038 -0.031 -0.025 -0.014 0.008 0.044 0.045 ...
##   ..$ III: num [1:1000] 0.064 0.065 0.076 0.08 0.072 0.071 0.106 0.104 0.092 0.081 ...
##   ..$ AVR: num [1:1000] 0.086 0.083 0.082 0.077 0.066 0.061 0.066 0.044 0.002 -0.004 ...
##   ..$ AVL: num [1:1000] -0.091 -0.09 -0.098 -0.098 -0.087 -0.084 -0.112 -0.1 -0.07 -0.059 ...
##   ..$ AVF: num [1:1000] 0.004 0.006 0.016 0.021 0.021 0.023 0.046 0.056 0.068 0.063 ...
##   ..$ V1 : num [1:1000] -0.069 -0.064 -0.058 -0.05 -0.045 -0.036 -0.029 -0.023 -0.015 -0.05 ...
##   ..$ V2 : num [1:1000] -0.031 -0.036 -0.034 -0.03 -0.027 -0.025 -0.012 0.003 0.018 0.009 ...
##   ..$ V3 : num [1:1000] 0 -0.003 -0.01 -0.015 -0.02 -0.009 0.005 0.018 0.021 0.018 ...
##   ..$ V4 : num [1:1000] -0.026 -0.031 -0.028 -0.023 -0.019 -0.014 -0.008 0.002 0.009 0.022 ...
##   ..$ V5 : num [1:1000] -0.039 -0.034 -0.029 -0.022 -0.018 -0.012 -0.007 -0.001 0.005 0.009 ...
##   ..$ V6 : num [1:1000] -0.079 -0.074 -0.069 -0.064 -0.058 -0.052 -0.048 -0.041 -0.038 -0.033 ...
##  $ sampling_frequency: num 100
##  $ time              : num [1:1000] 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 ...
##  - attr(*, "class")= chr "sig"

# empi_class <- empi_execute (
#   signal = signal,
#   write_to_file = TRUE,
#   path = NULL,
#   file_name = "00001_lr.db"
# )

# Reading a SQLite file in which all generated atom parameters are stored.
file <- system.file("extdata", "00001_lr.db", package = "MatchingPursuit")

# Create time-frequency map based on atoms.
out <- tf_map(
  x = file,
  channel = 2,
  increase_factor = 8
)
## Channel number: 2
## Total channels: 12
## Number of atoms: 50
## Sampling frequency: 100 Hz
## Epoch size (in points): 1000
## Signal length (in seconds): 10
## 
## Energy of the original signal:      6.93
## Energy of the reconstructed signal: 7.06
## reconstruction / original %:        101.88

plot(
  x = out_ecg,
  begin = 0,
  end = 10,
  panel_height = 1,
  zero_line = FALSE,
  small_squares = TRUE
)

9 One specific example

In this chapter, we present a specific data example adapted from the work of Durka (2007). The signal consists of a mixture of seven components: (a) four Gabor functions with different parameters, (b) a unit impulse, (c) a sinusoidal waveform, and (d) a chirp signal.

In the T–F map, all signal components—except for the chirp—are represented clearly and accurately (i.e., blobs for Gabor functions, a horizontal line for the sine wave, and a vertical line for the unit impulse). However, the chirp signal is decomposed into several separate blobs. This behavior arises from the discrete nature of the atom dictionary used in the Matching Pursuit algorithm, which prevents a continuous representation of a signal with smoothly varying frequency. This limitation (and, in some respects, a drawback) of the Matching Pursuit algorithm should be taken into account.

file <- system.file("extdata", "sample2.db", package = "MatchingPursuit")

out <- tf_map(
  x = file,
  channel = 1,
  mode = "sqrt",
  freq_divide = 1,
  increase_factor= 4,
  display_crosses = TRUE,
  display_atom_numbers = FALSE,
  out_mode = "plot",
  plot_signals = FALSE
)
## Channel number: 1
## Total channels: 1
## Number of atoms: 32
## Sampling frequency: 128 Hz
## Epoch size (in points): 1280
## Signal length (in seconds): 10
## 
## Energy of the original signal:      903.48
## Energy of the reconstructed signal: 910.88
## reconstruction / original %:        100.82

Back to menu

10 The Matching Pursuit algorithm in a nutshell

The Matching Pursuit algorithm is well-known and described in the literature. Its purpose is to approximate the analyzed signal using so-called atoms. (the text below is adapted from Kunik and Gramacki (2025)).

Given a signal \(f \in \mathbb{R}^n\), and a (possibly overcomplete) large redundant dictionary \(D =\{g_{\gamma}\}_{\gamma \in \Gamma}\) of normalized atoms \(\|g_{\gamma}\|=1\) Matching Pursuit finds a sparse signal representation

\[ f \approx \sum_{n = 0}^{N-1} a_n g_{\gamma_n}, \tag{1} \] where \(a_n \in \mathbb{R}\) are coefficients, \(g_{\gamma_n} \in D\) are atoms selected from the dictionary and \(N\) is the desired number of iterations (or stopping threshold). In most practical cases \(N \ll size(D)\). Also, \(g_{\gamma}\) is the dictionary atom indexed by \(\gamma\) and \(\Gamma\) is the set of all indices in the dictionary.

In the ideal case, the linear expansion (1) should include all atoms \(g_{\gamma_n}\) that represent the relevant structures of the signal \(f\). For real signals, such an ideal scenario is rarely possible, and some form of approximation is required. This task can be accomplished elegantly using the Matching Pursuit algorithm, which was first proposed by Mallat and Zhang (1993) in the context of signal processing.

Each atom \(g_{\gamma}\) is typically a time-frequency shifted, scaled version of a prototype function, such as the Gabor function (often called a Gaussian-windowed sinusoid). The dictionary is constructed to cover a wide range of time and frequency characteristics. The real-valued Gabor function has the following form:

\[ g_{\gamma}(t) = K(\gamma) e^{- \pi \left( \frac{t-\mu}{\sigma} \right) ^2} \cos(\omega (t - \mu) + \phi), \tag{2} \]

where \(\gamma = (\mu, \omega, \sigma, \phi)\) constitute a four-dimmensional space and \(K(\gamma)\) is such that \(||g_{\gamma}|| = 1\). It is easy to see that Gabor functions are constructed by multiplying Gaussian envelopes with cosine oscillations of different frequencies \(\omega\) and phases offset \(\phi\). By multiplying these two functions, we can obtain a wide variety of shapes depending on their parameters. A few examples of Gabor function are presented in figure below (in blue). The sinusoidal plane wave (in gray) is modulated by a Gaussian envelope (in red).

In the Matching Pursuit algorithm, the decomposition process is iterative. At each step, the algorithm selects an atom \(g_{\gamma_n}\) from the dictionary \(D\) that best matches the current residual signal \(R\). Formally, starting with the signal \(f\) at iteration \(n=0\) the initial residual and initial function approximation are

\[ R^0 = f \tag{3} \] and

\[ f^0 = 0. \tag{4} \]

For each iteration \(n = \{0,1,\ldots, N-1\}\) we find \(g_{\gamma_n} \in D\) such that the following inner product \(\langle \cdot, \cdot \rangle\) is maximized

\[ g_{\gamma_n} = \operatorname*{arg\,max}_{\gamma \in \Gamma} | \langle R^{n}, g_{\gamma} \rangle |. \tag{5} \] The coefficients \(a_n\) in (1) are

\[ a_n = \langle R^{n}, g_{\gamma_n} \rangle \tag{6} \]

and updated function approximation is defined as

\[ f^{n+1} = f^{n} + a_n g_{\gamma_n}. \tag{7} \]

Similarly, updated residual is defined as

\[ R^{n+1} = R^{n} - a_n g_{\gamma_n}. \tag{8} \]

After \(N\) iterations, the signal \(f\) is approximated as

\[ f \approx \sum_{n = 0}^{N-1} \langle R^n, g_{\gamma_n} \rangle g_{\gamma_n} = \sum_{n = 0}^{N-1} a_n g_{\gamma_n} \tag{9} \]

or equivalently

\[ f = \sum_{n = 0}^{N-1} \langle R^n, g_{\gamma_n} \rangle g_{\gamma_n} + R^{N}. \tag{10} \]

The procedure stops when \(\|R^{n+1}\|_2\) falls below a predefined threshold or when a fixed number of iterations has been reached

It should be noted that finding an optimal approximation (1) is an NP-hard problem. A suboptimal solution can be obtained using an iterative procedure, such as the MP algorithm.

Another key property of MP is energy conservation: the total energy of the signal is preserved in the MP decomposition.

Because the selected atoms are normalized, each MP iteration removes an amount of squared \(\ell_2\) energy equal to \(|\langle R^n,g_{\gamma_n}\rangle|^2\) from the current residual. Consequently, the signal energy can be decomposed as

\[ ||f||^2_2 = \sum_{n = 0}^{N-1} |\langle R^n, g_{\gamma_n} \rangle |^2 + ||R^{N}||^2_2. \tag{11} \]

This identity is sometimes referred to as the energy conservation property of Matching Pursuit.

Back to menu

11 The Orthogonal Matching Pursuit in a nutshell

The package also provides an implementation of the Orthogonal Matching Pursuit (OMP) algorithm. OMP algorithm is closely related to MP, but differs in how approximation coefficients are estimated.

In the classical MP algorithm, after an atom has been selected, only the residual signal is updated (8). The previously selected atoms and their coefficients remain unchanged and are not re-estimated in subsequent iterations.

OMP addresses this limitation by recomputing the coefficients of all selected atoms after each atom selection. As a result, the new residual is orthogonal to the subspace spanned by all selected atoms.

At each iteration, OMP selects the next atom \(g_{\gamma_n}\) by maximizing its correlation with the current residual, as in (7). The crucial difference from MP lies in the subsequent estimation of the coefficients. After the new atom has been added to the selected set, the coefficients are obtained by solving a least-squares problem

\[ \mathbf{a}^{(n)} = \operatorname*{arg\,min}_{\mathbf{c}} \left\| f - D_n \mathbf{c} \right\|_2^2, \tag{12} \]

where \(D_n\) denotes the matrix whose columns are the atoms selected up to and including iteration \(n\):

\[ D_n = [g_{\gamma_0},g_{\gamma_1},\ldots,g_{\gamma_n}]. \tag{13} \]

Assuming that \(D_n\) has full column rank, the least-squares solution is

\[ \mathbf{a}^{(n)} = (D_n^T D_n)^{-1}D_n^T f. \tag{14} \]

This coefficient re-estimation is the core difference between OMP and MP. Instead of simply adding the contribution of the newly selected atom, as in (7), OMP projects the original signal \(f\) orthogonally onto the subspace spanned by all currently selected atoms.

The new residual is then calculated as

\[ R^{n+1} = f - D_n \mathbf{a}^{(n)}. \tag{15} \]

Because the residual is orthogonal to the span of the selected atoms, its correlation with every selected atom is zero. Consequently, a previously selected atom cannot be selected again.

Compared with MP, OMP generally provides: 1. More accurate signal reconstruction for a given number of selected atoms than MP, 2. Potentially fewer selected atoms for a prescribed reconstruction accuracy, 3. More stable coefficient estimates in the presence of correlated atoms.

The increased accuracy comes at the cost of higher computational complexity because a least-squares problem must be solved after each iteration. In practical OMP implementations, including the implementation provided by our package, the least-squares problem need not be solved from scratch at each iteration. Instead, an incremental Cholesky factorization can be used to update the solution efficiently as new atoms are added to the selected set.

Back to menu

12 Pure R implementation of the MP and OMP algorithms

For educational and experimental purposes, the package also provides pure R implementations of the MP and OMP algorithms through the mp_core() and omp_core() functions, respectively.

The MP-R and OMP-R backends use a three-stage procedure consisting of dictionary generation or import with read_gabor_dict(), signal-dependent candidate atom selection with topk_atoms(), and sparse decomposition with mp_omp_execute().

sig_file <- system.file("extdata", "sample1.csv", package = "MatchingPursuit")
signal <- read_csv_signals(sig_file, col_names_in_csv = FALSE)

sampling_frequency <- signal$sampling_frequency
duration <- nrow(signal$signal) / sampling_frequency

xml_file <- system.file("extdata", "sample1.xml", package = "MatchingPursuit")

atoms_dict <- read_gabor_dict(
  xml_file, 
  sampling_frequency, 
  duration, 
  verbose = TRUE)
## Number of blocks: 10
## ===================================
## Block: 1
## windowLen = 17
## windowShift = 1
## fftSize = 32
## Atoms in block: 16128
## ===================================
## Block: 2
## windowLen = 27
## windowShift = 1
## fftSize = 64
## Atoms in block: 31936
## ===================================
## Block: 3
## windowLen = 41
## windowShift = 2
## fftSize = 128
## Atoms in block: 31488
## ===================================
## Block: 4
## windowLen = 65
## windowShift = 3
## fftSize = 128
## Atoms in block: 20480
## ===================================
## Block: 5
## windowLen = 103
## windowShift = 6
## fftSize = 256
## Atoms in block: 19712
## ===================================
## Block: 6
## windowLen = 163
## windowShift = 9
## fftSize = 512
## Atoms in block: 24576
## ===================================
## Block: 7
## windowLen = 259
## windowShift = 15
## fftSize = 512
## Atoms in block: 13312
## ===================================
## Block: 8
## windowLen = 409
## windowShift = 24
## fftSize = 1024
## Atoms in block: 13312
## ===================================
## Block: 9
## windowLen = 647
## windowShift = 38
## fftSize = 2048
## Atoms in block: 10240
## ===================================
## Block: 10
## windowLen = 1023
## windowShift = 61
## fftSize = 2048
## Atoms in block: 1024
## ===================================
## Total atoms: 182208

dict_topk <- topk_atoms(
  atoms_dict = atoms_dict,
  signal = signal,
  topk = 5000,
  verbose = TRUE
)
## topk_atoms(), step 1, calculating 182208 inner products...
## topk_atoms(), step 1 finished.
## topk_atoms(), step 2, signal 1 finished.

fit <- mp_omp_execute(
  mode = 'omp',          # or "mp" for classical Matching Pursuit
  dictionary = dict_topk,
  signal = signal,
  n_nonzero_coefs = 50,
  verbose = TRUE
)
## iteration: 1, selected atom: 1
## iteration: 2, selected atom: 75
## iteration: 3, selected atom: 200
## iteration: 4, selected atom: 248
## iteration: 5, selected atom: 719
## iteration: 6, selected atom: 543
## iteration: 7, selected atom: 1156
## iteration: 8, selected atom: 1559
## iteration: 9, selected atom: 2446
## iteration: 10, selected atom: 3084
## iteration: 11, selected atom: 4145
## iteration: 12, selected atom: 809
## iteration: 13, selected atom: 737
## iteration: 14, selected atom: 4487
## iteration: 15, selected atom: 1918
## iteration: 16, selected atom: 4984
## iteration: 17, selected atom: 4971
## iteration: 18, selected atom: 2873
## iteration: 19, selected atom: 2179
## iteration: 20, selected atom: 614
## iteration: 21, selected atom: 838
## iteration: 22, selected atom: 4108
## iteration: 23, selected atom: 4810
## iteration: 24, selected atom: 3922
## iteration: 25, selected atom: 4920
## iteration: 26, selected atom: 4417
## iteration: 27, selected atom: 4745
## iteration: 28, selected atom: 4740
## iteration: 29, selected atom: 4855
## iteration: 30, selected atom: 4948
## iteration: 31, selected atom: 3405
## iteration: 32, selected atom: 4323
## iteration: 33, selected atom: 4866
## iteration: 34, selected atom: 4017
## iteration: 35, selected atom: 4751
## iteration: 36, selected atom: 4735
## iteration: 37, selected atom: 3960
## iteration: 38, selected atom: 4124
## iteration: 39, selected atom: 4564
## iteration: 40, selected atom: 4925
## iteration: 41, selected atom: 4896
## iteration: 42, selected atom: 4378
## iteration: 43, selected atom: 4686
## iteration: 44, selected atom: 3160
## iteration: 45, selected atom: 4945
## iteration: 46, selected atom: 3699
## iteration: 47, selected atom: 4987
## iteration: 48, selected atom: 3678
## iteration: 49, selected atom: 4388
## iteration: 50, selected atom: 3866
## mp_omp_execute(): channel 1 processed.

plot(
  fit, 
  channel = 1, 
  freq_divide = 4
)
## Channel number: 1
## Total channels: 1
## Number of atoms: 50
## Sampling frequency: 1024 Hz
## Epoch size (in points): 1024
## Signal length (in seconds): 1
## 
## Energy of the original signal:      2746.16
## Energy of the reconstructed signal: 2491.18
## reconstruction / original %:        90.72

out <- mp_omp_pipeline(
   mode = 'omp',          # or "mp" for classical Matching Pursuit
   sig_file = sig_file,
   col_names_in_csv = FALSE,
   xml_file = xml_file,
   topk = 5000,
   n_nonzero_coefs = 50,
   verbose = FALSE
 )

 class(out)
## [1] "mp"

13 Relationship between MP-R / OMP-R and EMPI backends

The topk_atoms() function is designed to work with dictionaries exported from the EMPI program. Dictionaries generated by EMPI can be read using read_gabor_dict() and subsequently used as input for topk_atoms() and mp_omp_execute().

For meaningful comparisons between results generated by EMPI and MP-R/OMP-R results generated in R, it is recommended to run EMPI with the options

-o none --full-atoms-in-signal

which disable additional EMPI-specific optimization procedures and force generation of complete atoms within the analyzed signal. Details of these options can be found in the EMPI documentation (README.md file).

The external EMPI program is a multi-threaded C++ implementation of the Matching Pursuit algorithm with GPU support and CPU parallelization. The implementation supports optimal dictionaries, including the simulation of continuous (quasi-infinite) dictionaries (switch -o global). By default, atoms are allowed to extend beyond the signal boundaries, with the signal assumed to be zero outside its observed range. This behavior can be changed using the --full-atoms-in-signal switch, which restricts atoms to lie entirely within the signal.

The diagram below illustrates the three main workflows for running the MP-R, OMP-R, and EMPI backends.

                                       START
                                         |
                      -----------------------------------------
                      |                                       |
           MP-R / OMP-R workflow                         EMPI workflow
                      |                                       | 
         ---------------------------                          |            
         |                         |                          | 
 mp_omp_pipeline()         read_*_signals()            read_*_signals()
         |                         │                          │
         |                 read_gabor_dict()                  |             
         |                         │                          │
         |                   topk_atoms()                     │
         |                         │                          │
         |                 mp_omp_execute()             empi_execute()
         |                         |                          |
         ---------------------------                          |
                     |                                        |
                     ------------------------------------------
                                          |
                                   plot() / tf_map()                                   
                                   
read_*_signals() - select the appropriate function depending on the file format:
read_csv_signals(), read_edf_signals(), read_wfdb_signals()

Back to menu

14 XML-based atom dictionary specification

The XML file encodes the structure of a dictionary of basis functions (Gabor atoms), including window lengths, window shifts, and frequency grids used during signal analysis. A simple example illustrates how atom parameters are encoded in the XML file. Consider the following block:

<block>
<param name="windowLen" value="17"/>
<param name="windowShift" value="1"/>
<param name="fftSize" value="32"/>
</block>

Assume that the analyzed signal is sampled at \(1\;024\) Hz and has a duration of \(1\) second. We define:

• Sampling frequency: (\(f_s = 1\;024\ \text{Hz}\))
• Signal duration: (\(T = 1\ \text{sec.}\))
• Number of samples: (\(N = 1\;024\))
• Window length: (\(L = 17\))
• Window shift: (\(S = 1\))
• FFT size: \(32\)

<?xml version="1.0" encoding="ISO-8859-1"?>
<dict>
<block>
<param name="windowLen" value="17"/>
<param name="windowShift" value="1"/>
<param name="fftSize" value="32"/>
</block>
<block>
<param name="windowLen" value="27"/>
<param name="windowShift" value="1"/>
<param name="fftSize" value="64"/>
</block>
<block>
<param name="windowLen" value="41"/>
<param name="windowShift" value="2"/>
<param name="fftSize" value="128"/>
</block>
<block>
<param name="windowLen" value="65"/>
<param name="windowShift" value="3"/>
<param name="fftSize" value="128"/>
</block>
<block>
<param name="windowLen" value="103"/>
<param name="windowShift" value="6"/>
<param name="fftSize" value="256"/>
</block>
<block>
<param name="windowLen" value="163"/>
<param name="windowShift" value="9"/>
<param name="fftSize" value="512"/>
</block>
<block>
<param name="windowLen" value="259"/>
<param name="windowShift" value="15"/>
<param name="fftSize" value="512"/>
</block>
<block>
<param name="windowLen" value="409"/>
<param name="windowShift" value="24"/>
<param name="fftSize" value="1024"/>
</block>
<block>
<param name="windowLen" value="647"/>
<param name="windowShift" value="39"/>
<param name="fftSize" value="2048"/>
</block>
<block>
<param name="windowLen" value="1023"/>
<param name="windowShift" value="61"/>
<param name="fftSize" value="2048"/>
</block>
<block>
<param name="windowLen" value="1619"/>
<param name="windowShift" value="97"/>
<param name="fftSize" value="4096"/>
</block>
<block>
<param name="windowLen" value="2559"/>
<param name="windowShift" value="154"/>
<param name="fftSize" value="8192"/>
</block>
</dict>

<block>
  <param name="windowLen" value="2559"/>
  <param name="windowShift" value="154"/>
  <param name="fftSize" value="8192"/>
</block>

15 XML-based atom dictionary generator

The package includes a utility for generating XML-based atom dictionaries for MPTK-like sparse decomposition algorithms. The generator creates multiscale Gabor dictionaries with logarithmically distributed window lengths,automatically selecting window shifts and FFT sizes for each atom scale.

The generated dictionaries support multiresolution signal analysis by combining short atoms for transient components and long atoms for slowly varying structures. Long atoms are allowed through zero-padding, following the strategy commonly used in Matching Pursuit implementations.

# Generate a dictionary for a 4096-sample signal
dict <- generate_xml_dict (
   N = 4096,
   file = tempfile(fileext = ".xml")
)

dict
##    windowLen windowShift fftSize
## 1         17           1      64
## 2         25           2      64
## 3         37           2     128
## 4         55           3     128
## 5         83           5     256
## 6        121           7     256
## 7        179          11     512
## 8        265          16    1024
## 9        393          24    1024
## 10       577          35    2048
## 11       857          51    2048
## 12      1265          76    4096
## 13      1873         112    4096
## 14      2769         166    8192
## 15      4095         246    8192

Back to menu