Introduction to MatchingPursuit package

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:

more accurate signal reconstruction,
sparser representations,
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

Introduction to `MatchingPursuit` package

authors: Artur Gramacki and Jarosław Gramacki

1 R version and generation time

2 Loading the package

3 Installing the required external software

4 The purpose and functionality of the package

5 Quick start

Step 1

Step 2

Step 3

Step 3.bis

6 Electroencephalogram (EEG) analysis

Step 1

Step 2 (sometimes required)

Step 3 (usually required)

Step 4 (usually required)

Step 5

Step 6

Step 7 (optional)

7 Electrocardiography (ECG) analysis

Step 1

Step 2

Step 3

Step 4 (optional)

8 One specific example

9 The Matching Pursuit algorithm in a nutshell

10 The Orthogonal Matching Pursuit (OMP) in a nutshell

11 OMP workflow in the package

Step 1

Step 2

Step 3

Step 4

Combine steps 1–4 into a single pipeline

12 Relationship between OMP and EMPI

13 XML-base atom dictionary specification

Number of window positions

Number of frequency bins

The total number of atoms

Multiple-block dictionaries

Bibliogaphy