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

[Back to menu](#top)

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

[Back to menu](#top)

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

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

[Back to menu](#top)

### Step 2 (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}
# Original signal
eeg <- read_edf_signals(file, resampling = FALSE, from = 0, to = 10)
signal_eeg <- eeg$signal
sampling_frequency <- eeg$sampling_frequency
sampling_frequency

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

```{r sampling.results, include = TRUE, echo = FALSE, warning = FALSE, fig.width = 7, fig.height = 4}
par(mfrow = c(2, 1), pty = "m", mai = c(0.8, 0.5, 0.5, 0.5))

# Not-filtered signal (raw signal).
plot(
  signal_eeg[, 1],
  type = "l",
  panel_first = grid(),
  main = "Original signal (256 Hz), channel #1",
  xlab = "sample points",
  ylab = "",
  col = "blue"
)

# Signal after filtering.
plot(
  signal_eeg_r[, 1],
  type = "l",
  panel_first = grid(),
  main = "Signal after downsampling (64 Hz), channel #1",
  xlab = "sample points",
  ylab = "",
  col = "blue"
)

par(old_par)
```

[Back to menu](#top)

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

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

```

[Back to menu](#top)

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

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.

```{r eeg_filtering, include = TRUE, echo = TRUE, warning = FALSE}
# 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
`r ch <- 1; colnames(signal_eeg_m_f)[ch]`), 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.

```{r no_filtering_and_filtering_plot, include = TRUE, echo = FALSE, warning = FALSE, fig.width = 7, fig.height = 4}
ch <- 1
par(mfrow = c(2, 1), pty = "m", mai = c(0.8, 0.5, 0.5, 0.5))

# Not-filtered signal (raw signal).
plot(
  signal_eeg_m[, ch],
  type = "l",
  panel_first = grid(),
  main = paste(colnames(signal_eeg_m)[ch], " (raw signal, channel #1)", sep = ""),
  xlab = "sample points",
  ylab = "",
  col = "blue"
)

# Signal after filtering.
plot(
  signal_eeg_m_f[, ch],
  type = "l",
  panel.first = grid(),
  main = paste(colnames(signal_eeg_m)[ch], " (filtered signal, channel #1)", sep = ""),
  xlab = "sample points",
  ylab = "",
  col = "blue"
)

par(old_par)
```

[Back to menu](#top)

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

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

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

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.

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

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.

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

[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_ecg_signals, include = TRUE, echo = TRUE, warning = FALSE}
file <- system.file("extdata", "00001_lr.hea", package = "MatchingPursuit")

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

```

## 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_fun(n, f, mean = 2, phase = 0, sigma = 1, frequency = 50, normalization = F)
g2 <- gabor_fun(n, f, mean = 8, phase = 0, sigma = 1.5, frequency = 25, normalization = F)
g3 <- gabor_fun(n, f, mean = 2, phase = 0, sigma = 2, frequency = 30, normalization = F)
g4 <- gabor_fun(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 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.

```{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$ 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
\ll 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-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_fun, 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_fun(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$cosinus, type = "l", col = "grey")
  lines(gb$t, gb$gabor, col="blue", lwd = 2)
  
  s <- s + gb$gabor
}

par(old_par)
```

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

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

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

## OMP workflow in the package 

The OMP implementation uses a three-stage procedure.

### Step 1 {.unnumbered}

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

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

### 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,
  sf = sf,
  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 `omp_execute()` function.

```{r omp_execution, include = TRUE, echo = TRUE, warning = FALSE}
fit <- omp_execute(
  dictionary = dict_topk,
  signal = signal,
  sf = sf,
  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 MP decompositions:

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

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

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

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

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

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

## Bibliogaphy