# Introduction to `MatchingPursuit` package

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

## 1 R version and generation time

    ## R version:    4.5.1
    ## Generated on: 23-czerwiec-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. Loading files in EDF and EDF+ formats
2. Loading files in the WFDB (WaveForm DataBase) format
3. Performing three standard EEG montages (bipolar, referential, and average)
4. Pre-filtering signals using notch, low-pass, high-pass, band-pass, and band-stop filters
5. Plotting signals
6. Executing the Matching Pursuit (MP) algorithm
7. Loading an XML-base atom dictionary
8. Executing the Orthogonal Matching Pursuit (OMP) algorithm
9. Plotting time–frequency maps based on MP and OMP results

Back to menu

## 5 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 2
    ##  $ 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

Back to menu

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

Back to menu

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

## 6 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.
    read_edf_params(file)
    ##       channel_name frequency no_of_samples length_sec
    ## 1              Fp1       256          2560         10
    ## 2              Fp2       256          2560         10
    ## 3               F3       256          2560         10
    ## 4               F4       256          2560         10
    ## 5               F7       256          2560         10
    ## 6               F8       256          2560         10
    ## 7               Fz       256          2560         10
    ## 8               C3       256          2560         10
    ## 9               C4       256          2560         10
    ## 10              Cz       256          2560         10
    ## 11              T3       256          2560         10
    ## 12              T5       256          2560         10
    ## 13              T4       256          2560         10
    ## 14              T6       256          2560         10
    ## 15              P3       256          2560         10
    ## 16              P4       256          2560         10
    ## 17              Pz       256          2560         10
    ## 18              O1       256          2560         10
    ## 19              O2       256          2560         10
    ## 20 EDF_Annotations        57           627         11

Back to menu

### Step 2 (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.

    # Original signal
    eeg <- read_edf_signals(file, resampling = FALSE, from = 0, to = 10)
    signal_eeg <- eeg$signal
    sampling_frequency <- eeg$sampling_frequency
    sampling_frequency
    ## [1] 256

    # Resampled signal
    eeg_resampled <- read_edf_signals(file, resampling = TRUE, sf_new = 64, from = 0, to = 10)
    signal_eeg_r <- eeg_resampled$signal
    sampling_frequency_r <- eeg_resampled$sampling_frequency
    sampling_frequency_r
    ## [1] 64

Back to menu

### Step 3 (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.
    eeg_m <- eeg_montage(eeg, montage_type = c("bipolar"), bipolar_pairs = pairs)
    signal_eeg_m <- eeg_m$signal

    # Original signal (first 6 rows, first 6 channels).
    signal_eeg[1:6, 1:6]
    ##         Fp1       Fp2        F3        F4         F7        F8
    ## 1 -22.71201 -25.55728 -23.71532 -24.88713 -20.368393 -24.38547
    ## 2 -12.67505 -18.02861 -14.51318 -16.68830  -7.823155 -16.02198
    ## 3 -15.18333 -19.86674 -16.68830 -18.69876 -10.833093 -18.02861
    ## 4 -13.84686 -18.86342 -15.51650 -17.52695  -9.328124 -17.02529
    ## 5 -14.51318 -19.53358 -16.18665 -18.19710 -10.162942 -17.69162
    ## 6 -14.01153 -19.20041 -15.85349 -17.86011  -9.661285 -17.18996

    # Signal after banana montage (first 6 rows, first 6 channels).
    signal_eeg_m[1:6, 1:6]
    ##       Fp2_F4      F4_C4     C4_P4     P4_O2   Fp1_F3     F3_C3
    ## 1 -0.6701516 -0.1684953 -1.167979 0.3331611 1.003313 -1.003313
    ## 2 -1.3403032 -0.3331611 -2.006625 0.6663222 1.838130 -2.006625
    ## 3 -1.1679785 -0.3369905 -1.838130 0.5016563 1.504969 -1.673464
    ## 4 -1.3364738 -0.3369905 -1.838130 0.6701516 1.669635 -1.838130
    ## 5 -1.3364738 -0.1684953 -2.006625 0.5016563 1.673464 -1.673464
    ## 6 -1.3403032 -0.3331611 -2.010455 0.5054858 1.841960 -1.841960

Back to menu

### Step 4 (usually required)

EEG signals are rarely analyzed without prior filtering. Using the `filters_coeff()` 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 <- filters_coeff(
       sf = sampling_frequency,
       notch = c(49, 51),
       lowpass = 40,
       highpass = 1,
    )

    # Filtering input signals.
    signal_eeg_m_f <- signal_eeg_m

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

For the selected channel (after the bipolar montage, it is named Fp2_F4), both the original and filtered signals are displayed. The effect of filtering is clearly beneficial, as noise—typically without diagnostic significance—has been removed from the signal.

Back to menu

### 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 <- list(signal_eeg_m_f, sampling_frequency)
    # names(sig) <- c("signal", "sampling_frequency")

    # empi_class <- empi_execute (
    #   signal = sig,
    #   empi_options = NULL,
    #   write_to_file = TRUE,
    #   path = NULL,
    #   file_name = "EEG_bipolar_filtered.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_bipolar_filtered.db", package = "MatchingPursuit")

    # Generate time-frequency map based on MP 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: 256 Hz
    ## Epoch size (in points): 2560
    ## Signal length (in seconds): 10
    ##
    ## Energy of the original signal:      231102.49
    ## Energy of the reconstructed signal: 221487.18
    ## reconstruction / original %:        95.84

### Step 7 (optional)

The package also includes a function for displaying EEG signals. To do this, you can use `plot_eeg()`, the `S3` method for the generic `plot()` function. This function requires an object of class `edf`, created with `read_edf_signals()`.

`panel_height` parameter controls the vertical distance between individual signals. If `NULL`, its value is set automatically so that all signals are clearly visible and do not overlap. The calculated value of parameter `panel_height` is displayed below the figure.

    temp <- eeg_m
    temp$signal <- signal_eeg_m_f

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

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

In the next example, we show a plot of the same EEG signal as above. However, this time, we display the unfiltered version. Undesirable artifacts (noise and drift) are clearly visible.

    temp <- eeg_m
    temp$signal <- signal_eeg_m

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

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

Back to menu

## 7 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_ecg_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 <- list(
      signal = data.frame(out_ecg$signal),
      sampling_frequency = out_ecg$sampling_frequency
    )

    str(signal)
    ## List of 2
    ##  $ 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

### 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 MP 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_ecg()`, the `S3` method for the generic `plot()` function. This function requires an object of class `ecg`, created with `read_ecg_signals()`.

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

## 8 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 MP algorithm, which prevents a continuous representation of a signal with smoothly varying frequency. This limitation (and, in some respects, a drawback) of the MP 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

## 9 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\) MP finds a sparse signal representation

\[ f \approx \sum_{n = 0}^{N} 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 size(D)$. Also, \(g_{\gamma}\) is the dictionary atom indexed by \(\gamma\) and \(\Gamma\) is the set of all indices in the directory.

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 MP 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 MP 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\}\) we find \(g_{\gamma_n} \in D\) such that the following inner product \(\langle \cdot, \cdot \rangle\) is maximized

\[ g_{\gamma_n} = \operatorname*{arg\,min}_{\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} \langle R^n, g_{\gamma_n} \rangle g_{\gamma_n} = \sum_{n = 0}^{N} a_n g_{\gamma_n} \tag{9} \]

or equivalently

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

The procedure stops when \(\|R^{n+1}\|\) is below a threshold or after a fixed number of iterations.

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.

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

Back to menu

## 10 The Orthogonal Matching Pursuit (OMP) in a nutshell

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

In the classical MP algorithm, after selecting an atom, only the residual signal is updated (8). Previously selected atoms and their coefficients \(a_n\) remain unchanged. and are not recalculated in the next steps.v

OMP addresses this limitation by recomputing the coefficients \(a_n\) of all selected atoms after each iteration. This makes the residue orthogonal to all selected atoms.

At each iteration, the next atom \(g_{\gamma_n}\) is selected according to (5). At this point, a fundamental difference with the MP algorithm appears. The coefficients are obtained from a least-squares problem:

\[ a_n = \operatorname*{arg\,min}_{\mathbf{c}} \left\| f - D_n \mathbf{c} \right\|_2^2 \]

and, assuming \(D_n\) has full column rank \[ a_n = (D_n^T D_n)^{-1}D_n^Tf \]

where \(D_k\) is the matrix containing the selected atoms. This is the core of OMP. Instead of just adding the contribution of the new atom (see (7)), OMP projects the original signal \(f\) orthogonally onto the space spanned by all currently selected atoms in \(D_n\). Because the residual is orthogonal to the span of the selected atoms, a previously selected atom has zero correlation with the residual and therefore cannot be selected again.

Now, the new residual is calculated by subtracting the new approximation from the original signal. The new residual is now completely orthogonal to all columns that have been selected so far

\[ R^{n+1} = f - D_n a_n \]

Compared with MP, OMP usually provides:

1. more accurate signal reconstruction,
2. sparser representations,
3. fewer selected atoms for a comparable approximation quality.

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 ours), the least-squares problem is not solved from scratch at each iteration. Instead, an incremental Cholesky factorization is used to update the solution efficiently.

Back to menu

## 11 OMP workflow in the package

The OMP implementation uses a three-stage procedure.

### Step 1

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

Note: Strictly speaking, the `read_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 13 for details.

    sig_file <- system.file("extdata", "sample3.csv", package = "MatchingPursuit")
    sample3 <- read_csv_signals(sig_file, col_names_in_csv = TRUE)

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

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

    atoms_dict <- read_dict(
      xml_file,
      sf,
      duration,
      verbose = TRUE)
    ## Number of blocks: 7
    ## ===================================
    ## Block: 1
    ## windowLen = 17
    ## windowShift = 1
    ## fftSize = 32
    ## Atoms in block: 3840
    ## ===================================
    ## Block: 2
    ## windowLen = 27
    ## windowShift = 1
    ## fftSize = 64
    ## Atoms in block: 7360
    ## ===================================
    ## Block: 3
    ## windowLen = 41
    ## windowShift = 2
    ## fftSize = 128
    ## Atoms in block: 6912
    ## ===================================
    ## Block: 4
    ## windowLen = 65
    ## windowShift = 3
    ## fftSize = 128
    ## Atoms in block: 4096
    ## ===================================
    ## Block: 5
    ## windowLen = 103
    ## windowShift = 6
    ## fftSize = 256
    ## Atoms in block: 3328
    ## ===================================
    ## Block: 6
    ## windowLen = 161
    ## windowShift = 9
    ## fftSize = 512
    ## Atoms in block: 2816
    ## ===================================
    ## Block: 7
    ## windowLen = 255
    ## windowShift = 15
    ## fftSize = 512
    ## Atoms in block: 256
    ## ===================================
    ## Total atoms: 28608

### 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,
      sf = sf,
      topk = 5000,
      verbose = TRUE
    )
    ## topk_atoms(), step 1, calculating 28608 inner products...
    ## topk_atoms(), step 1 finished.
    ## topk_atoms(), step 2, signal 1 finished.
    ## topk_atoms(), step 2, signal 2 finished.
    ## topk_atoms(), step 2, signal 3 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 `omp_execute()` function.

    fit <- omp_execute(
      dictionary = dict_topk,
      signal = signal,
      sf = sf,
      n_nonzero_coefs = 50,
      verbose = TRUE
    )
    ## omp_execute(): channel 1 processed.
    ## omp_execute(): channel 2 processed.
    ## omp_execute(): channel 3 processed.

### Step 4

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

    plot(fit, channel = 3)
    ## Channel number: 3
    ## Total channels: 3
    ## Number of atoms: 50
    ## Sampling frequency: 128 Hz
    ## Epoch size (in points): 256
    ## Signal length (in seconds): 2
    ##
    ## Energy of the original signal:      21.7
    ## Energy of the reconstructed signal: 21.35
    ## reconstruction / original %:        98.39

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

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

    out <- run_omp_pipeline(
       sig_file = sig_file,
       col_names_in_csv = TRUE,
       xml_file = xml_file,
       topk = 5000,
       n_nonzero_coefs = 50,
       verbose = FALSE
     )

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

## 12 Relationship between OMP and EMPI

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

For meaningful comparisons between MP results generated by EMPI and OMP 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).

Summarising, the typical workflow is as follows:

`read_dict()` \(\rightarrow\) `topk_atoms()` \(\rightarrow\) `omp_execute()` \(\rightarrow\) `plot()`/`tf_map()`

Back to menu

## 13 XML-base 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

## Bibliogaphy

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.