---
title: "Introduction to `MatchingPursuit` package"
author: "authors: Artur Gramacki and Jarosław Gramacki"
bibliography: bibliography.bib
output: 
  rmarkdown::html_vignette:
    toc: true
    number_sections: true  
    fig_caption: true
vignette: >
  %\VignetteIndexEntry{Introduction to `MatchingPursuit` package}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 72
---

<style>

body {
  background-color: #ffffff;
}

h1, h2, h3, h4 {
  border-bottom: none;
}

pre {
  border: 1px solid #3b82f6;
  border-radius: 6px;  
}

</style>

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE
)

old_par <- par("mfrow", "mai", "pty")
options(width = 999)
```

## R version and generation time

```{r R_version, include = TRUE, echo = FALSE}
out <- version
cat(
  "R version:    ", 
  out$major, ".", 
  out$minor, 
  "\n", 
  "Generated on: ", 
  format(Sys.time(), "%d-%B-%Y"), 
  sep = ""
)
```

## Loading the package

```{r library}
library(MatchingPursuit)
```

## Installing the required external software

The below `empi_install()` function downloads **Enhanced Matching
Pursuit Implementation** external program (or **EMPI** for short), see
@Rozanski-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.

```{r empi_locate, include = TRUE, echo = TRUE}
empi_locate()
```

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.

```{r empi_install, include = TRUE, echo = TRUE}
# empi_install()
```

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

```{r empi_check, include = TRUE, echo = TRUE}
empi_check()
```

## 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-1993, @Pati-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  

## 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](https://github.com/develancer/empi).

[Back to menu](#top)

## Quick start

### Step 1 {.unnumbered}

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.

```{r 7_non_stationary, include = TRUE, echo = TRUE, warning = FALSE}
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)
```

```{r 7_non_stationary_plot, include = TRUE, echo = FALSE, warning = FALSE, fig.width = 7, fig.height = 7}
range <- range(signal)
par(mfrow = c(8, 1), pty = "m", mai = c(0.2, 0.4, 0.2, 0.1))

plot(t, signal$x, type = "l", col = "blue", main = "", xlab = "", ylab = "")
plot(t, x1, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "Linear chirp")
plot(t, x2, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "Nonlinear chirp")
plot(t, x3, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "AM")
plot(t, x4, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "FM")
plot(t, x5, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "Decreasing amplitude")
plot(t, x6, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "Frequency modulated sine wave")
plot(t, x7, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "Increasing amplitude")

par(old_par)
```

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.

```{r sample1_csv, include = TRUE, echo = TRUE, warning = FALSE}
# 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)

signal <- read_csv_signals(file)
str(signal)
```

### Step 2 {.unnumbered}

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

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

```{r empi_execute, include = TRUE, echo = TRUE, warning = FALSE}
# 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 {.unnumbered}

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

```{r sample1_db, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7}
# 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"
)
```

### Step 3.bis {.unnumbered}

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](#important_note). Uncomment it to view the result.

```{r plot_empi, include = TRUE, echo = TRUE, warning = FALSE}
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](#top)

## 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 {.unnumbered}

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.

```{r read_edf_params, include = TRUE, echo = TRUE, warning = FALSE}
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) {.unnumbered}

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.

```{r eeg_filtering, include = TRUE, echo = TRUE, warning = FALSE}
# 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) {.unnumbered}

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.

```{r read_edf_signal_resampling, include = TRUE, echo = TRUE, warning = FALSE}
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.

```{r filtering_and__resampling_plot, include = TRUE, echo = FALSE, warning = FALSE, fig.width = 7, fig.height = 5}
old_par <- par("mfrow", "mai", "mgp")
par(mfrow = c(3, 1), pty = "m", mai = c(0.35, 0.5, 0.4, 0.1), mgp = c(1.5, 0.5, 0))

rg <- range(c(signal_eeg[, 1], signal_eeg_f[, 1], signal_eeg_f_r[, 1]))

plot(
  time,
  signal_eeg[, 1],
  type = "l",
  xlab = "Time (s)",
  ylab = "Amplitude",
  main = "Original EEG signal, channel Fp1",
  col = "blue",
  ylim = rg,
  panel.first = grid()
)
abline(h = 0, col = "gray")

plot(
  time,
  signal_eeg_f[, 1],
  type = "l",
  xlab = "Time (s)",
  ylab = "Amplitude",
  main = "Filtered EEG signal, channel Fp1",
  col = "blue",
  ylim = rg,
  panel.first = grid()
)
abline(h = 0, col = "gray")

plot(
  time_64,
  signal_eeg_f_r[, 1],
  type = "l",
  xlab = "Time (s)",
  ylab = "Amplitude",
  main = "Filtered and downsampled EEG signal (64 Hz), channel Fp1",
  col = "blue",
  ylim = rg,
  panel.first = grid()
)
abline(h = 0, col = "gray")
par(old_par)
```

### Step 4 (usually required) {.unnumbered}

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.

```{r eeg_montage, include = TRUE, echo = TRUE, warning = FALSE}
# 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]

# Signal after banana montage (first 6 rows, first 6 channels).
signal_eeg_f_r_m[1:6, 1:6]
```

### Step 5 {.unnumbered}

**Important note:** The code below has been commented out. See the
explanation given [here](#important_note).

```{r eeg_empi_execute, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7}
# 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 {.unnumbered}

Generating the final time-frequency map for the selected channel
(`r ch <- 1; colnames(signal_eeg_f_r_m)[ch]`). 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.

```{r eeg_tf_map, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7}
# 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"
)
```

### Step 7 (optional) {.unnumbered}

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.

```{r plot_eeg_filt_resamp_mont, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7}
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"
)
```

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

```{r plot_orig_eeg, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7}

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

[Back to menu](#top)

## 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 {.unnumbered}

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](https://physionet.org/content/ptb-xl/1.0.3/).

```{r read_wfdb_signals, include = TRUE, echo = TRUE, warning = FALSE}
file <- system.file("extdata", "00001_lr.hea", package = "MatchingPursuit")

out_ecg <- read_wfdb_signals(file)

head(out_ecg$signals)
out_ecg$sampling_frequency
out_ecg$lead_names
out_ecg$record_name
```

```{r ecg_signal, include = TRUE, echo = TRUE, warning = FALSE}
# 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)
```

### Step 2 {.unnumbered}

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

```{r ecg_empi_execute, include = TRUE, echo = TRUE, warning = FALSE}
# empi_class <- empi_execute (
#   signal = signal,
#   write_to_file = TRUE,
#   path = NULL,
#   file_name = "00001_lr.db"
# )
```

### Step 3 {.unnumbered}

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

```{r 00001_lr, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7}
# 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
)
```

### Step 4 (optional) {.unnumbered}

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

```{r plot_ecg, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7}
plot(
  x = out_ecg,
  begin = 0,
  end = 10,
  panel_height = 1,
  zero_line = FALSE,
  small_squares = TRUE
)

```

[Back to menu](#top)

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

```{r chirp_def, include = TRUE, echo = FALSE, warning = FALSE}
n <- 1280
f <- 128
t1 = n / f
t <- seq(from = 0, to = n - 1, by = 1) / f

g1 <- gabor_atom(n, f, mean = 2, phase = 0, sigma = 1, frequency = 50, normalization = F)
g2 <- gabor_atom(n, f, mean = 8, phase = 0, sigma = 1.5, frequency = 25, normalization = F)
g3 <- gabor_atom(n, f, mean = 2, phase = 0, sigma = 2, frequency = 30, normalization = F)
g4 <- gabor_atom(n, f, mean = 7, phase = 0, sigma = 0.5, frequency = 15, normalization = F)
imp <- rep(0, n)
imp[n / 2] <- 4
sine <- 0.2 * sin(2 * pi * 3 * t)
chirp <- signal::chirp(t = t, f0 = 0, t1 = t1, f1 = 50, form = c("linear"), phase = 0)
signal <- g1$gabor + g2$gabor + g3$gabor + g4$gabor + imp  + chirp + sine
```

```{r chirp_plot, include = TRUE, echo = FALSE, warning = FALSE, fig.width = 7, fig.height = 7}
par(mfcol = c(8, 1), pty = "m", mai = c(0.2, 0.1, 0.2, 0.1)) 

ylim <- c(-2, 2)
plot(t, signal,   ylim = ylim, xaxt = "s", yaxt = "n", bty = "o", col = "blue",  type = "l", xlab = "", ylab = "", main = "The sum of the below functions")
plot(t, g1$gabor, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "",   main = "Gabor: f = 50Hz, mean = 2, sigma = 1, phase = 0")
plot(t, g2$gabor, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "",   main = "Gabor: f = 25Hz, mean = 8, sigma = 1.5, phase = 0")
plot(t, g3$gabor, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "",   main = "Gabor: f = 30Hz, mean = 2, sigma = 2, phase = 0")
plot(t, g4$gabor, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "",   main = "Gabor: f = 15Hz, mean = 7, sigma = 0.5, phase = 0")
plot(t, imp,      ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "",   main = "Unit impulse")
plot(t, sine,     ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "",   main = "Sine wave: 3 Hz")
plot(t, chirp,    ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "",   main = "Linear-frequency chirp (0Hz - 50Hz)")

par(old_par)
```

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.

```{r chirp_plot_tf, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 6}
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
)

```

[Back to menu](#top)

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

```{r gabor_atom, include = TRUE, echo = FALSE, fig.width = 7, fig.height = 5, fig.align = 'center'}
N <- 512
fs <- 256
# normalization = T --> signal's norm = 1

par(mfrow = c(2,2), pty = "m", mai = c(0.4, 0.4, 0.3, 0.2))

n <- 4

sigmas <- c(0.5, 0.2, 0.8, 0.5)
frequencies <- c(14, 8, 4, 1)
phases <- c(0, 1, 1.5, -2)
means = c(0.5, 0.8, 1, 1.5)

s <- rep(0, N)

for (i in c(1:n)) {
  sigma <- sigmas[i]
  frequency <- frequencies[i]
  phase <- phases[i]
  mean <- means[i]
  main <- latex2exp::TeX(paste(
    "$\\mu=$", means[i], ", ", 
    "$\\sigma=$", sigmas[i], ", ",
    "$\\f=$", frequencies[i], ", ",
    "$\\phi=$", phases[i], 
    sep = ""
    )
  )
  
  gb <- gabor_atom(N, fs, mean, phase, sigma, frequency, normalization = F)

  plot(
    gb$t, 
    gb$gauss, type="l", ylim = c(-1, 1), col = "red", 
    xlab= "", ylab= "",
    xaxt = "t", yaxt = "t", bty = "or",
    cex.axis = 1, lwd = 2,
    main = main)
  lines(gb$t, gb$cosine, type = "l", col = "grey")
  lines(gb$t, gb$gabor, col="blue", lwd = 2)
  
  s <- s + gb$gabor
}

par(old_par)
```

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.


```{r restore_par, include = FALSE}
par(old_par)
```

[Back to menu](#top)

## 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](#top)

## 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 {.unnumbered}

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](https://mptk.irisa.fr/)  program 
and can also be generated by the EMPI program; see Section 15 for details.

```{r read_gabor_dict, include = TRUE, echo = TRUE, warning = FALSE}
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)
```

### Step 2 {.unnumbered}

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.

```{r topk_atoms, include = TRUE, echo = TRUE, warning = FALSE}
dict_topk <- topk_atoms(
  atoms_dict = atoms_dict,
  signal = signal,
  topk = 5000,
  verbose = TRUE
)
```

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 {.unnumbered}

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

```{r omp_execution, include = TRUE, echo = TRUE, warning = FALSE}
fit <- mp_omp_execute(
  mode = 'omp',          # or "mp" for classical Matching Pursuit
  dictionary = dict_topk,
  signal = signal,
  n_nonzero_coefs = 50,
  verbose = TRUE
)
```

### Step 4 {.unnumbered}

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

```{r omp_plot, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7}
plot(
  fit, 
  channel = 1, 
  freq_divide = 4
)

```

### Combine steps 1–4 into a single pipeline {.unnumbered}

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

```{r mp_omp_pipeline, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7}
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)
```

[Back to menu](#top)

## 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](#top)

## 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 {.unnumbered}

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 {.unnumbered}


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 {.unnumbered}


$$
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 {.unnumbered}


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](#top)

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

```{r generate_dict, include = TRUE, echo = TRUE, warning = FALSE}
# Generate a dictionary for a 4096-sample signal
dict <- generate_xml_dict (
   N = 4096,
   file = tempfile(fileext = ".xml")
)

dict
```

[Back to menu](#top)

## Bibliography
