Introduction to MatchingPursuit package

Step 1 - Input signal preparation

Let us construct an example signal by combining seven non-stationary components. The resulting signal (highlighted in blue below) will be used to demonstrate the basic functionality of the package.

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)

Data must be stored in a data frame: rows represent samples for all channels, and columns represent channels. Our first demo dataset consists of only one channel (one column). The read_csv_signals() function checks whether the data has the correct structure. The first line of the file must contain two numbers: the sampling rate in Hz (freq) and the signal length in seconds (sec). This allows verification that the file actually contains freq*sec samples.

# 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)
signal
## Signal object (class 'sig')
## --------------------------------------------------
## Samples:            1024
## Channels:           1
## Sampling frequency: 1024 Hz
## Duration:           1 s

Step 2 - Native R MP decomposition

The input data (signal) is passed as an argument to the mp_omp_execute() function, which generates the object of class mp containing all atom parameters.

For Gabor-based decomposition, mp_omp_execute() provides the high-level native R interface. It prepares the Gabor dictionary, selects candidate atoms, performs MP or OMP decomposition, and returns an object of class "mp". In this example, classical Matching Pursuit is used; setting mode = "omp" runs Orthogonal Matching Pursuit instead.

fit_mp <- mp_omp_execute(
  mode = "mp",
  signal = signal,
  n_nonzero_coefs = 25
)
## mp_omp_execute(): method = "mp", channel = 1 Successfully processed.

summary(fit_mp)
## Summary of Matching Pursuit object (class 'mp')
## --------------------------------------------------------
## Samples:            1024
## Channels:           1
## Sampling frequency: 1024 Hz
## Duration:           1 s
## 
## Per-channel decomposition:
## 
## Explained energy = 1 - residual energy / signal energy
## For MP, this measure is preferred because the reconstruction 
## and residualare not generally orthogonal, so reconstruction 
## energy / signal energy is not equivalent to explained energy.
## 
## Channel 1 (Channel 1): 25 atoms
## Signal energy:         2746.16
## Reconstruction energy: 2594.91
## Residual energy:       231.225
## Explained energy:      91.58%

Step 3 - Time-frequency map

It is now time to generate the final time-frequency (T-F) map for the selected channel. To display the T–F map, you can also use plot.mp(), the S3 method for the generic plot() function. This function requires an object of class mp, such as an object returned by mp_omp_execute() or empi_execute().

plot(fit_mp)

By comparing the two signal waveforms below the T-F map, it can be seen that the original signal and the reconstructed signal differ only minimally. Also, basic signal parameters are displayed, along with information about the number of atoms into which the input signal was decomposed. Additionally, the energy of the input signal and the reconstructed signal (from the atoms) is calculated. The results show that 91.58% of the original signal’s energy is “explained” by the generated atoms. Naturally, increasing the number of generated atoms will likely bring the energy of the reconstructed signal closer to 100%.

Alternatively, you can use the tf_map() function to display a time-frequency map. The function provides greater control over the visualization process and offers additional functionality, such as saving time-frequency maps to a file. See the function documentation for details.

In the below plot, the centers of the atoms (in terms of time and frequency coordinates) are marked with the numbers of successive atoms, sorted from highest to lowest energy.

out <- tf_map(
  x = fit_mp,
  channel = 1,
  freq_divide = 4,
  atom_centers = "numbers"
)

Step 4 - Optional EMPI-based decomposition

The same signal can also be decomposed using the optional external Enhanced Matching Pursuit Implementation (EMPI) backend. Note that EMPI must be installed separately before empi_execute() can be used.

Important note: The code in the chunk below has been commented out because CRAN’s verification rules prohibit automatic binary downloads. empi_execute() function requires that EMPI is installed, otherwise it terminates with an error message. Therefore, this function cannot be executed here. This is why the output object from the empi_execute() has been generated in advance and included in the package as rds file. In the next chunk, the plot() and tf_map() functions use this file as input.

Notice the empi_options parameter in the empi_execute() function. You can omit this parameter, and the EMPI program will run with the default values set in the function ("-o local --gabor -i 50"). It is also worth noting that the EMPI program offers a wide range of configuration options. Details can be found in the README.md file located in the directory where the EMPI program was installed. In our example, parameters were set to instruct the program to find 25 atoms.

# sample1_empi_out <- empi_execute(
#   signal = signal,
#   empi_options = "-o local --gabor -i 25",
# )

# saveRDS(sample1_empi_out, file = "sample1.rds")

Now, let’s load a previously generated result (stored in the sample1.rds file).

file <- system.file("extdata", "sample1.rds", package = "MatchingPursuit")
sample1_empi_out <- readRDS(file = file)

plot(sample1_empi_out)

The T-F map is similar to the previous one, although some differences in atom selection can be observed. The reconstructed signal explains a slightly larger proportion of the original signal energy, reaching 92.69%.

Back to menu

Step 1 - Reading EEG data

Reading an example EEG signal (EDF file). The signal is 10 seconds long and consists of 20 channels (19 plus a special channel called the annotation channel that does not contain EEG signal data). The sampling rate is 256 Hz for each channel with EEG data. The channels have standard names.

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

Step 2 - Filtering

EEG signals are rarely analyzed without prior filtering. Using the design_filters() function, you can define the filter parameters and then apply the filter to the signal. The filter parameters listed below use typical values recommended in the literature for EEG signal analysis.

# 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
}

Step 3 - Resampling

Sometimes it is necessary to change the sampling frequency (increase — upsampling or decrease — downsampling). In the example below, the original sampling frequency is reduced from 256 Hz to 128 Hz.

signal_eeg_f_r <- resample_signal(signal = signal_eeg_f, p = 1, q = 2)
time_128 <- seq(0, nrow(signal_eeg_f_r) - 1) / (sampling_frequency / 2)
sampling_frequency_r <- 128

The effect of the preprocessing (filtering and resampling) is illustrated below. Although filtering removes baseline drift and high-frequency noise and downsampling reduces the sampling frequency, the overall morphology of the EEG waveform is well preserved.

Step 4 - Montage construction

A bipolar montage is created (the classical double banana montage), where each channel compares two adjacent electrodes. In the first step, you define the pairs of electrodes to be connected using the pairs list. In the second step, the eeg_montage() function generates the required montage.

# 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.009169519 -0.00948931 -0.007902186 -0.009404775 -0.004150085 -0.006489879
## 2 -0.007879380 -0.02545971 -0.017611303 -0.019963222 -0.007967950 -0.027958014
## 3  0.078689947  0.17633851  0.135141209  0.153953511  0.073793375  0.189075962
## 4 -4.405347474 -5.25777306 -4.292400132 -5.062491209 -2.274244020 -4.099115184
## 5 -7.123860061 -9.75573535 -7.216335597 -8.494470758 -3.112619764 -7.582372940
## 6 -4.904350386 -7.37213153 -4.761021228 -5.880097268 -1.128596436 -5.172483701

# 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,] -8.453465e-05 -0.009038717 -0.003565737  0.0001613474 -0.001267334 -0.0003177297
## [2,] -5.496491e-03  0.030559501  0.005492579  0.0012699229  0.009731924 -0.0059268445
## [3,]  2.238500e-02 -0.150521483 -0.025526483 -0.0073991820 -0.056451262  0.0364474510
## [4,] -1.952819e-01 -2.676997475 -1.303608531  0.1215588493 -0.112947342 -0.4015900412
## [5,] -1.261265e+00 -2.573360624 -2.040162295  0.2925304240  0.092475536 -0.8818594888
## [6,] -1.492034e+00 -2.397153899 -1.985208885  0.2483389685 -0.143329158 -0.6516086820

Step 5 - MP decomposition

Important note: The code below has been commented out. See the explanation given here.

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

# # To make the RDS file smaller (CRAN requirement), we select only one channel.
# sig <- as_sig(signal_eeg_f_r_m[, 2], sampling_frequency_r)

# eeg_empi_out <- empi_execute (
#   signal = sig,
#   empi_options = NULL
# )
# 
# saveRDS(eeg_empi_out, file = "eeg_empi_out.rds")

Step 6 - Time-frequency map

Generating the final time-frequency map for the selected channel (Fp2_F4). The centers of atoms (in the sense of time and frequency coordinates) are now marked with white crosses.

Comparing the two signal waveforms under the T-F map, it can be seen that the original and reconstructed signals differ only minimally.

Below the plot, basic signal parameters are displayed, along with information about the number of atoms into which the input signal was decomposed.

Additionally, the energy of the input signal and the reconstructed signal (from the generated atoms) is calculated. The results show that nearly all of the energy of the original signal (96.18%) is “explained” by the generated atoms.

file <- system.file("extdata", "eeg_empi_out.rds", package = "MatchingPursuit")
eeg_empi_out <- readRDS(file = file)

plot(eeg_empi_out)

Step 7 - Displaying EEG signals using plot() (optional)

The package also provides the plot.edf() function, an S3 method for plot() function, for visualization of multichannel EEG recordings. The function accepts an object of class edf returned by read_edf_signals() and displays the selected time interval using a conventional stacked EEG layout.

Since the preprocessing functions are designed to operate on generic multichannel signal matrices, they are independent of the edf class and can be applied to signals imported from different file formats. Therefore, before using plot.edf(), the processed signals are combined with the corresponding metadata into an object of class edf.

The panel_height argument controls the vertical spacing between adjacent channels. If it is NULL, an appropriate value is determined automatically to prevent overlap between signals. The selected value is reported in the R console.

edf_processed <- structure(
  list(
    signal = as.data.frame(signal_eeg_f_r_m),
    sampling_frequency = sampling_frequency_r,
    time = time_128,
    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: 67.5

The same approach can be used to visualize the original EEG recording before filtering, resampling, and montage construction.

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

Back to menu

Step 1 - Reading ECG data

Reading an example ECG signal (.dat and .hea files). The signal is 10 seconds long and consists of 12 channels. The sampling rate is 100 Hz. The channels have standard names. The data comes from the repository available at PhysioNet.

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

out_ecg <- read_wfdb_signals(file)

head(out_ecg$signal)
##           I     II   III   AVR    AVL   AVF     V1     V2     V3     V4     V5     V6
## [1,] -0.119 -0.055 0.064 0.086 -0.091 0.004 -0.069 -0.031  0.000 -0.026 -0.039 -0.079
## [2,] -0.116 -0.051 0.065 0.083 -0.090 0.006 -0.064 -0.036 -0.003 -0.031 -0.034 -0.074
## [3,] -0.120 -0.044 0.076 0.082 -0.098 0.016 -0.058 -0.034 -0.010 -0.028 -0.029 -0.069
## [4,] -0.117 -0.038 0.080 0.077 -0.098 0.021 -0.050 -0.030 -0.015 -0.023 -0.022 -0.064
## [5,] -0.103 -0.031 0.072 0.066 -0.087 0.021 -0.045 -0.027 -0.020 -0.019 -0.018 -0.058
## [6,] -0.097 -0.025 0.071 0.061 -0.084 0.023 -0.036 -0.025 -0.009 -0.014 -0.012 -0.052
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"

Step 2 - MP decomposition

The input data (signal) is passed as an argument to the empi_execute() function, which generates the object of class mp containing all atom parameters.

Important note: The code below has been commented out. See the explanation given here.

# ecg_empi_out <- empi_execute(
#   signal = out_ecg
# )
# 
# saveRDS(ecg_empi_out, file = "00001_lr.rds")

Step 3 - Time-frequency map

It is now time to generate the final time-frequency (T-F) map for the selected channel.

# Read the previously generated decomposition result.
file <- system.file("extdata", "00001_lr.rds", package = "MatchingPursuit")
ecg_empi_out <- readRDS(file = file)

# Create time-frequency map based on atoms.
out <- tf_map(
  x = ecg_empi_out,
  channel = 1,
  verbose = TRUE
)
## Channel number:             1
## Total channels:             12
## Number of atoms:            50
## Sampling frequency:         100 Hz
## Epoch size (in points):     1000
## Signal length (in seconds): 10
## 
## Signal energy:              11.89
## Reconstruction energy:      11.95
## Residual energy:            0.55
## Explained_energy:           95.41%

Step 4 - ECG printout (optional)

The package also includes a function for displaying ECG signals in a layout corresponding to standard paper ECG printouts. A typical ECG paper layout was used, with a small grid of 0.04 s × 0.1 mV and a large grid of 0.20 s × 0.5 mV. To do this, you can use plot.wfdb(), the S3 method for the generic plot() function. This function requires an object of class ecg, created with read_wfdb_signals().

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

Back to menu

Core MP and OMP functions

The native R implementation is built around two core functions, mp_core() and omp_core(), which implement the Matching Pursuit and Orthogonal Matching Pursuit algorithms, respectively. These functions form the general-purpose sparse decomposition layer of the package and are independent of any particular atom family or dictionary construction method.

Both functions operate on user-defined dictionaries represented as numeric matrices, with candidate atoms stored in columns. This design allows the core algorithms to be used not only within the Gabor-based workflow provided by the package, but also with arbitrary dictionaries constructed independently by the user.

To demonstrate that the core solvers are independent of the Gabor-specific workflow, the following examples use a simple user-defined matrix dictionary.

dictionary <- matrix(
  c(
     0.098, -0.308, -0.342, -0.894,
     0.928,  0.674, -0.270,  0.283,
     0.145, -0.326,  0.810, -0.114,
    -0.170, -0.466,  0.377, -0.036,
     0.281, -0.357,  0.105, -0.327
  ),
  nrow = 5,
  byrow = TRUE
)

# Although mp_core() and omp_core() normalize dictionary atoms internally, 
# the atoms are normalized explicitly here so that the coefficients used to 
# construct the signal have a direct interpretation.
dictionary <- sweep(dictionary, 2, sqrt(colSums(dictionary^2)), "/")

# Signal constructed from two correlated dictionary atoms
signal <- 1.5 * dictionary[, 1] - 0.8 * dictionary[, 2]

After defining the dictionary and signal, the MP and OMP decompositions can be performed.

fit_mp <- mp_core(
  dictionary = dictionary,
  signal = signal,
  n_nonzero_coefs = 2
)

fit_omp <- omp_core(
  dictionary = dictionary,
  signal = signal,
  n_nonzero_coefs = 2
)

The decomposition results can then be compared in terms of the selected atoms, estimated coefficients, and residual error.

fit_mp$support
## [1] 1 2
fit_omp$support
## [1] 1 2

fit_mp$coefs
## [1]  1.0783989 -0.5778156
fit_omp$coefs
## [1]  1.5 -0.8

sum(fit_mp$residual^2)
## [1] 0.1283816
sum(fit_omp$residual^2)
## [1] 4.622232e-31

The two examples illustrate the different behavior of MP and OMP when dictionary atoms are correlated. Both algorithms identify the same relevant atoms, but OMP re-estimates all active coefficients after each new atom is added. As a result, OMP recovers the original coefficients more accurately and produces a smaller residual error. In contrast, MP updates the approximation sequentially and does not revise previously estimated coefficients, which can lead to a less accurate reconstruction for the same number of selected atoms.

High-level decomposition interface

The mp_omp_execute() function provides a convenient high-level interface for Gabor-based decomposition. After loading the signal, the user only needs to select the decomposition mode and specify the main analysis parameters.

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

The decomposition can then be performed using either OMP or classical MP. The returned object is ready for direct visualization with plot().

fit <- mp_omp_execute(
   mode = 'omp',          # or "mp" for classical Matching Pursuit
   signal = signal,
   n_nonzero_coefs = 50,
   topk = 10000,
   verbose = FALSE
 )
## mp_omp_execute(): method = "omp", channel = 1 Successfully processed.

plot(fit)

Gabor dictionary preparation

For educational and experimental purposes, the individual steps performed internally by the mp_omp_execute() function can also be carried out explicitly. These steps include generating an XML file containing the dictionary parameters using generate_xml_dict(), reading this file using read_gabor_dict(), and preselecting candidate atoms using topk_gabor_atoms().

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

xml_file <- tempfile(fileext = ".xml")

dict <- generate_xml_dict(
   N = signal_length,
   file = xml_file
)

Reading the Gabor dictionary

The read_gabor_dict() function reads the XML-based dictionary specification and returns the corresponding candidate Gabor atom definitions. Depending on the specification, the resulting set may contain tens or hundreds of thousands of candidate atoms. For format details, see the section XML-based atom dictionary specification.

atoms_dict <- read_gabor_dict(
  xml_file = xml_file, 
  sampling_frequency = sampling_frequency, 
  duration = duration, 
  full_atoms_in_signal = FALSE,
  verbose = FALSE
)

dim(atoms_dict)
## [1] 322154      7
head(atoms_dict)
##      block time_sample   time_sec freq_bin freq_hz window_len fft_size
## [1,]     1          -8 -0.0078125        0       0         17       64
## [2,]     1          -8 -0.0078125        1      16         17       64
## [3,]     1          -8 -0.0078125        2      32         17       64
## [4,]     1          -8 -0.0078125        3      48         17       64
## [5,]     1          -8 -0.0078125        4      64         17       64
## [6,]     1          -8 -0.0078125        5      80         17       64

Signal-dependent Gabor atom selection

The topk_gabor_atoms() function evaluates all dictionary atoms using phase-invariant complex projections and selects the atoms with the highest similarities to the analysed signal.

This step substantially reduces the size of the optimization problem. Only the selected atoms are retained and converted into real-valued Gabor atoms with optimal phase estimates.

dict_topk <- topk_gabor_atoms(
  atoms_dict = atoms_dict,
  signal = signal,
  topk = 10000,
  verbose = TRUE
)
## topk_gabor_atoms(), step 1, calculating 322154 inner products...
## topk_gabor_atoms(), step 1 finished.
## topk_gabor_atoms(), step 2 finished.
## 10000 out of 322154 atoms selected successfully.

Validation using a synthetic sparse signal

To validate the behavior of the general sparse decomposition solvers, we use a synthetic signal with a known sparse representation. The procedure applies to both MP and OMP; here, OMP is used as a representative example.

First we construct the dictionary. The dictionary contains sinusoidal atoms at several frequencies, while only three atoms are used to generate the signal. This allows the selected atoms, estimated coefficients, and reconstruction error to be compared directly with the known ground truth.

Each column of the matrix represents one dictionary atom. The atoms are normalized to unit L2 norm so that the coefficients used to generate the signal can be compared directly with the coefficients estimated by OMP.

set.seed(1)

N <- 256
t <- (0:(N - 1)) / N

dictionary <- cbind(
  sin(2 * pi * 3  * t),  cos(2 * pi * 3  * t),  sin(2 * pi * 7  * t),
  cos(2 * pi * 7  * t),  sin(2 * pi * 12 * t),  cos(2 * pi * 12 * t),
  sin(2 * pi * 20 * t),  cos(2 * pi * 20 * t)
)

# Unit-L2 normalization of dictionary atoms.
dictionary <- sweep(dictionary, 2, sqrt(colSums(dictionary^2)), "/")

colnames(dictionary) <- c(
  "sin_3", "cos_3",  "sin_7", "cos_7",  
  "sin_12", "cos_12",  "sin_20", "cos_20"
)

Second, we generate a signal with known sparse representation. The signal is generated using only three non-zero coefficients. Therefore, the true support consists of the atoms cos_3, sin_12, and cos_20.

true_coef <- c(0, 1.5, 0, 0, -0.8, 0, 0, 0.5)
true_atoms <- c("cos_3", "sin_12", "cos_20")

signal <- as.vector(dictionary %*% true_coef)

Now, we can do the OMP decomposition. Because the signal was generated exactly from three dictionary atoms, OMP should recover the same support and the corresponding coefficients.

out <- omp_core(
  dictionary = dictionary,
  signal = signal,
  n_nonzero_coefs = 3
)

colnames(out$selected_atoms)
## [1] "cos_3"  "sin_12" "cos_20"

out$coefs
## [1]  1.5 -0.8  0.5

For the noise-free signal, the selected support should match the known support and the relative reconstruction error should be close to zero. The reconstruction error is at the level of numerical precision, indicating exact recovery of the noise-free signal.

setequal(colnames(out$selected_atoms), true_atoms)
## [1] TRUE

relative_error <- sqrt(sum((signal - out$reconstruction)^2)) / sqrt(sum(signal^2))
relative_error
## [1] 7.733811e-16

Now, we can examine the robustness to additive noise. The noisy example demonstrates whether the same sparse support can still be identified when the observed signal is perturbed. Since the clean signal is known in this synthetic example, the reconstructed signal can also be compared directly with the noise-free ground truth.

signal_noisy <- signal + rnorm(N, sd = 0.1)

out_noisy <- omp_core(
  dictionary = dictionary,
  signal = signal_noisy,
  n_nonzero_coefs = 3
)

colnames(out_noisy$selected_atoms)
## [1] "cos_3"  "sin_12" "cos_20"

out_noisy$coefs
## [1]  1.6412326 -0.7505156  0.3159247

setequal(colnames(out_noisy$selected_atoms), true_atoms)
## [1] TRUE

Despite additive Gaussian noise, OMP recovered the same three-atom support, while the reconstruction remained close to the known clean signal.

relative_residual_error <- sqrt(
  sum((signal_noisy - out_noisy$reconstruction)^2)
  ) / sqrt(sum(signal_noisy^2))

relative_residual_error
## [1] 0.6353605

relative_clean_error <- sqrt(
  sum((signal - out_noisy$reconstruction)^2)
  ) / sqrt(sum(signal^2))

relative_clean_error
## [1] 0.1338779

Back to menu

Relationship between MP-R, OMP-R, and EMPI backends

MP-R and OMP-R implement Matching Pursuit and Orthogonal Matching Pursuit, respectively. EMPI is a high-performance C++ backend specialized in Gabor-based Matching Pursuit and is designed for substantially faster decomposition of large dictionaries.

The diagram below summarizes the three available decomposition workflows and their relationship within the package. Blue boxes denote R functions, whereas peach boxes represent workflow components, inputs, or returned objects. Panel (b) shows the main steps performed internally by mp_omp_execute().

The general matrix-based MP/OMP workflow operates directly on a user-defined numeric dictionary and input signal through mp_core() or omp_core(). It is independent of the Gabor-specific dictionary construction utilities and can therefore be used with arbitrary matrix dictionaries. These functions return low-level decomposition results, giving the user direct control over the dictionary and the decomposition procedure.

The native R Gabor-based MP/OMP workflow, implemented by mp_omp_execute(), provides a higher-level interface for Gabor-based MP and OMP decomposition. It integrates dictionary specification or generation, dictionary reading, channel-specific atom preselection, decomposition with mp_core() or omp_core(), and aggregation of the results into an object of class "mp". If no XML dictionary specification is supplied, generate_xml_dict() creates one internally; otherwise, the user-provided XML file is used. The specification is processed by read_gabor_dict(), after which topk_gabor_atoms() selects candidate Gabor atoms for each signal channel.

The external EMPI-based MP workflow provides an alternative MP implementation through empi_execute(). In this case, decomposition is performed by the external high-performance C++ EMPI backend rather than by the native R MP/OMP core functions. The wrapper integrates the external decomposition results into the same package-level "mp" representation, allowing them to be handled using the same downstream visualization functions, including plot() and tf_map().

Thus, the three workflows differ mainly in the dictionary representation and decomposition backend: direct matrix-based MP/OMP for arbitrary user-defined dictionaries, the integrated native Gabor MP/OMP workflow, and the external EMPI-based Gabor MP workflow. The latter two return a common "mp" object, facilitating consistent visualization and interpretation within the package.

Back to menu

Number of window positions

With a shift of one sample (\(S = 1\)), the number of possible window positions is

\[ N_{\text{windows}} = N - L + 1 = 1\;024 - 17 + 1 = 1\;008 \]

Number of frequency bins

For an FFT size of \(32\), and because the analyzed signal is real-valued, it is sufficient to consider the non-negative frequencies from 0 to the Nyquist frequency. Thus, the number of frequency bins is

\[ N_{\text{freq}} = \frac{32}{2} + 1 = 17 \] For a sampling frequency of \(1024\) Hz, the frequency resolution is

\[ \Delta f = \frac{1024}{32} = 32\ \text{Hz}. \]

Therefore, the frequencies are:

\[ 0, 32, 64, 96, 128, 160, 192, 224, 256, 288, 320, 352, 384, 416, 448, 480, 512\ \text{Hz} \]

The total number of atoms

\[ N_{\text{atoms}} = N_{\text{windows}} \times N_{\text{freq}} = 1\;008 \times 17 = 17\;136 \]

In other words, for each of the \(1\;008\) possible positions of the \(17\)-sample window, \(17\) atoms corresponding to different FFT frequencies are generated.

Multiple-block dictionaries

A practical dictionary usually consists of multiple blocks with different window lengths, shifts, and FFT sizes. The calculations for each block are analogous to those shown above. The results for the example dictionary are summarized in the table below. As above, the sampling frequency is (\(f_s = 1\;024\ \text{Hz}\)) and the signal duration is (\(T = 1\ \text{sec.}\)), giving the number of samples (\(N = 1\;024\)). In total, the XML file defines parameters for \(185\;365\) atoms.

<?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>
</dict>

An interesting special case are the two last blocks:

<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>

The last two blocks illustrate a boundary case resulting from the finite signal length. For windowLen = \(1\;023\), the window can be placed only once within the 1024-sample signal, resulting in \(1\;025\) atoms corresponding to the available frequency bins. For the final block (windowLen = \(1\;619\)), the window length exceeds the signal length, so no fully contained atom can be generated. This behavior follows directly from the full_atoms_in_signal = TRUE convention, under which only atoms whose complete support lies within the signal boundaries are included.

A regular pattern can be observed in the dictionary construction: the ratio of windowShift to windowLen is nearly constant across all blocks (approximately 0.06). This means that consecutive windows overlap by about 94%.

Furthermore, as the window length increases, the frequency resolution also increases. This is a direct consequence of the time–frequency uncertainty principle - the longer the window, the better the ability to distinguish closely spaced frequencies, but at the expense of time localization.

Therefore, it is beneficial to analyze a larger number of frequency components for long windows, as they provide meaningful frequency information. In contrast, for short windows, a very dense frequency grid would not contribute much additional information because the frequency resolution is fundamentally limited by the window length itself.

In summary, the dictionary provides an approximately uniform coverage of the time–frequency plane - short windows are associated with many temporal positions and relatively few frequency bins, whereas long windows have fewer time positions but a much denser frequency sampling.

Back to menu