# Introduction to `MatchingPursuit` package

#### authors: Artur Gramacki and Jarosław Gramacki

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

### Step 1

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

### Step 2

The input data (`signal`) is passed as an argument to the `empi_execute()` function, which generates the final output file in `SQLite` format (`sample1.db`) containing all atom parameters.

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 can not be executed here. That is why the `sample1.db` file has been generated in advance and included in the package. In the next chunk, the `empi2ft()` function uses this file as input (the `x` parameter).

Notice the `empi_options` parameter in the `empi_execute()` function. You can specify `NULL` for 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 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.

    # 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"
    # )

### Step 3

It is now time to generate the final time-frequency (T-F) map for the selected channel. 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.

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.

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 atoms) is calculated. The results show that 94.05% 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%.

    # 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

### Step 3.bis

To display the T–F map, you can also use `plot_empi()`, the `S3` method for the generic `plot()` function. This function requires an object of class `empi`, created with `empi_execute()`. Due to CRAN regulations, the code below has been commented out. See the explanation here. Uncomment it to view the result.

    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)

Back to menu

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

### Step 1

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 (usually required)

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 (sometimes required)

Sometimes it is necessary to change the sampling frequency (increase — upsampling or decrease — downsampling). The `read_edf_signals()` function provides this functionality. In the example below, the original sampling frequency is reduced from 256 Hz to 64 Hz.

    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

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 (usually required)

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

### Step 5

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"

    # 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"
    # )

### Step 6

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 (97.75%) is “explained” by the generated atoms.

    # 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

### Step 7 (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 `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_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

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

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

### Step 1

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

### Step 2

The input data (`signal`) is passed as an argument to the `empi_execute()` function, which generates the final output file in `SQLite` format (`00001_lr.db`) containing all atom parameters.

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

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

### Step 3

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

    # 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

### Step 4 (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

## 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()`.

### Step 1

Dictionary generation. A dictionary of candidate Gabor atoms is generated with `read_gabor_dict()`. The dictionary may contain tens or hundreds of thousands of atoms.

Note: Strictly speaking, the `read_gabor_dict()` function does not create an atom dictionary itself. Instead, it reads an XML-based atom dictionary specification that defines the parameters and structure used to construct the dictionary. The format is compatible with the MPTK program program and can also be generated by the EMPI program; see Section 15 for details.

    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

### Step 2

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

    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.

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.

### Step 3

OMP decomposition. The reduced dictionary is passed to `mp_omp_execute()` function.

    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.

### Step 4

The result is an object of class `mp`, making it fully compatible with the same visualization functions used for Matching Pursuit decompositions:

    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

### Combine steps 1–4 into a single pipeline

For simplicity, all four steps described above have been encapsulated in a single function, `mp_omp_pipeline()`, which returns an object of class `mp`.

    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"

Back to menu

## 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\)

### 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\), the number of positive frequency bins is

\[ N_{\text{freq}} = \frac{32}{2} = 16 \] Only positive frequencies are considered, excluding the Nyquist component. \(512\) is omitted, so we have \(16\) different frequencies. Therefore, the frequencies are:

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

### The total number of atoms

\[ N_{\text{atoms}} = N_{\text{windows}} \times N_{\text{freq}} = 1\;008 \times 16 = 16\;128 \]

In other words, for each of the \(1\;008\) possible positions of the \(17\)-sample window, \(16\) 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 table below. In total, the XML file defines parameters for \(564\;416\) 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>
    <block>
    <param name="windowLen" value="2559"/>
    <param name="windowShift" value="154"/>
    <param name="fftSize" value="8192"/>
    </block>
    </dict>

windowLen | windowShift | fftSize | number of windows | number of frequencies | number of atoms
--- | --- | --- | --- | --- | ---
17 | 1 | 32 | 2544 | 16 | 40 704
27 | 1 | 64 | 2534 | 32 | 81 088
41 | 2 | 128 | 1260 | 64 | 80 640
65 | 3 | 128 | 832 | 64 | 53 248
103 | 6 | 256 | 410 | 128 | 52 480
163 | 9 | 512 | 267 | 256 | 68 352
259 | 15 | 512 | 154 | 256 | 39 424
409 | 24 | 1024 | 90 | 512 | 46 080
647 | 39 | 2048 | 50 | 1024 | 51 200
1023 | 61 | 2048 | 26 | 1024 | 26 624
1619 | 97 | 4096 | 10 | 2048 | 20 480
2559 | 154 | 8192 | 1 | 4096 | 4 096

An interesting special case is the last block:

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

For a signal containing \(2\;560\) samples, a window of length \(2\;559\) can be placed only once. Consequently, this block generates exactly \(1 \times 4\;096 = 4\;096\) atoms. The total number of atoms is

\[ 40\;704 + 81\;088 + \cdots + 4\;096 =564\;416 \]

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

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

## Bibliography

Durka, P. 2007. Matching Pursuit and Unification in EEG Analysis. Engineering in Medicine and Biology. Boston: Artech House.
Elad, M. 2010. Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing. Springer. https://doi.org/10.1007/978-1-4419-7011-4.
Kunik, M., and A. Gramacki. 2025. “Deep Learning Epileptic Seizure Detection based on the Matching Pursuit Algorithm and its Time–Frequency Graphical Representation.” International Journal of Applied Mathematics and Computer Science 35 (4): 617–30. https://doi.org/10.61822/amcs-2025-0044.
Mallat, S., and Z. Zhang. 1993. “Matching Pursuits With Time-Frequency Dictionaries.” IEEE Transactions on Signal Processing 41 (12): 3397–3415. https://doi.org/10.1109/78.258082.
Pati, Y. C., R. Rezaiifar, and P. S. Krishnaprasad. 1993. “Orthogonal Matching Pursuit: Recursive Function Approximation with Applications to Wavelet Decomposition.” In Proceedings of the 27th Asilomar Conference on Signals, Systems and Computers, 40–44. https://doi.org/10.1109/ACSSC.1993.342465.
Różański, P. T. 2024. “empi: GPU-Accelerated Matching Pursuit with Continuous Dictionaries.” ACM Transactions on Mathematical Software 50 (3): 1–17. https://doi.org/10.1145/3674832.