Package {MatchingPursuit}

Type:	Package
Title:	Processing Time Series Data Using the Matching Pursuit Algorithm
Version:	1.3.0
Maintainer:	Artur Gramacki <a.gramacki@gmail.com>
Description:	Provides tools for analysing and decomposing time series data using the Matching Pursuit (MP) algorithm, a greedy signal decomposition technique that represents complex signals as a linear combination of simpler functions (called atoms) selected from a redundant dictionary. Support for the Orthogonal Matching Pursuit (OMP) variant of the classical MP algorithm is also provided. For more details see Mallat and Zhang (1993) <doi:10.1109/78.258082>, Pati et al. (1993) <doi:10.1109/ACSSC.1993.342465>, Elad (2010) <doi:10.1007/978-1-4419-7011-4> and Różański (2024) <doi:10.1145/3674832>.
SystemRequirements:	external tool (installed via empi_install() function). The package uses the implementation of the Matching Pursuit algorithm (Enhanced Matching Pursuit Implementation; EMPI) by Piotr T. Różański, available at https://github.com/develancer/empi.
Imports:	edf, signal, RSQLite, imager, raster, graphics, grDevices, utils, digest, EGM, xml2
Suggests:	knitr, rmarkdown, latex2exp, remotes
VignetteBuilder:	knitr
Depends:	R (≥ 4.0.0)
License:	GPL-3
BugReports:	https://github.com/artur-gramacki/MatchingPursuit/issues
Encoding:	UTF-8
RoxygenNote:	7.3.3
NeedsCompilation:	no
Packaged:	2026-09-17 10:21:08 UTC; Artur
Author:	Artur Gramacki [aut, cre], Jarosław Gramacki [ctb], Piotr T. Różański [ctb] (provided technical guidance on the operation of the EMPI software)
Repository:	CRAN
Date/Publication:	2026-09-17 10:50:07 UTC

Sparse Time-Series Decomposition Using Matching Pursuit and Orthogonal Matching Pursuit

Description

Tools for sparse decomposition and analysis of time-series signals using Matching Pursuit (MP) and Orthogonal Matching Pursuit (OMP). Both algorithms represent a signal as a sparse linear combination of atoms selected from a dictionary.

Details

The package provides native R implementations of both MP and OMP through mp_core() and omp_core(). These functions operate on arbitrary user-defined dictionaries represented as numeric matrices, with candidate atoms stored in columns.

For Gabor-based time-frequency analysis, the package provides a dedicated workflow based on generate_xml_dict(), read_gabor_dict(), topk_gabor_atoms(), and the high-level mp_omp_execute() interface. Decomposition results can be visualized using plot() and tf_map().

In addition to the native R implementations, the package can use Enhanced Matching Pursuit Implementation (EMPI) as an optional external high-performance backend for Gabor-based Matching Pursuit. EMPI is implemented in C++ and supports optimized CPU and GPU execution.

The package also supports direct loading and processing of physiological signals stored in EDF/EDF(+) and WFDB (WaveForm DataBase) formats, facilitating analysis of EEG and ECG recordings.

Author(s)

Maintainer: Artur Gramacki a.gramacki@gmail.com (ORCID)

Other contributors:

• Jarosław Gramacki j.gramacki@gmail.com (ORCID) [contributor]

• Piotr T. Różański piotr@develancer.pl (ORCID) (provided technical guidance on the operation of the EMPI software) [contributor]

References

Durka, P. J. (2007). Matching Pursuit and Unification in EEG Analysis. Artech House, Engineering in Medicine and Biology. Boston. ISBN: 978-1596932497

Elad, M. (2010). Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing. Springer. ISBN 978-1-4419-7010-7, doi:10.1007/978-1-4419-7011-4

Gramacki, A. & Kunik, M. (2025). Deep learning epileptic seizure detection based on matching pursuit algorithm and its time-frequency graphical representation. International Journal of Applied Mathematics & Computer Science, vol. 35, no. 4, pp. 617-630, doi:10.61822/amcs-2025-0044

Mallat, S. & Zhang, Z. (1993). Matching Pursuits with Time-Frequency Dictionaries. IEEE Transactions on Signal Processing, vol. 41, no. 12, pp. 3397-3415, doi:10.1109/78.258082

Pati, Y.C. & Rezaiifar, R. & Krishnaprasad, P.S. (1993). Orthogonal Matching Pursuit: Recursive Function Approximation with Applications to Wavelet Decomposition. Proceedings of the 27th Asilomar Conference on Signals, Systems and Computers, vol. 1, pp. 40-44 doi:10.1109/ACSSC.1993.342465

Różański, P.T. (2024). empi: GPU-Accelerated Match ing Pursuit with Continuous Dictionaries. ACM Transactions on Mathematical Software, vol.50, no. 3, pp. 1-17, doi:10.1145/3674832

See Also

Useful links:

• Report bugs at https://github.com/artur-gramacki/MatchingPursuit/issues

Convert a signal to a sig object

Description

Creates an object of class sig from signal data already available in R. The function provides a convenient way to prepare signals for subsequent processing and decomposition without importing them from a file.

Usage

as_sig(signal, sampling_frequency)

Arguments

Details

The time vector is generated automatically from the number of signal samples and the specified sampling frequency. The resulting object has the same basic structure as objects returned by read_csv_signals().

Value

An object of class sig, which is a list containing:

• signal: A data frame containing the signal values.

• sampling_frequency: The sampling frequency in Hz.

• time: A numeric vector containing the time coordinates of the signal samples in seconds, starting at 0.

See Also

read_csv_signals, read_edf_signals, read_wfdb_signals

Examples

# Single-channel signal
x <- rnorm(1000)
sig <- as_sig(x, sampling_frequency = 100)
str(sig)

# Multi-channel signal
x <- cbind(
  channel1 = rnorm(1000),
  channel2 = rnorm(1000)
)
sig <- as_sig(x, sampling_frequency = 100)
str(sig)

Clear MatchingPursuit cache

Description

Deletes all files in the MatchingPursuit cache directory.

Usage

clear_cache()

Value

Logical scalar. Returns TRUE if all files were successfully removed, and FALSE otherwise. The return value is invisible.

Examples

if (interactive()) {
  clear_cache()
}

Design Butterworth filters

Description

Designs notch, low-pass, high-pass, band-pass, and band-stop Butterworth filters for a specified sampling frequency.

Usage

design_filters(
  sampling_frequency = 256,
  notch = c(49, 51),
  notch_order = 2,
  lowpass = 30,
  lowpass_order = 4,
  highpass = 1,
  highpass_order = 4,
  bandpass = c(0.5, 40),
  bandpass_order = 4,
  bandstop = c(0.5, 40),
  bandstop_order = 4
)

Arguments

Value

A list containing the designed Butterworth filter objects:

notch

Notch filter used to remove a specific narrow frequency band.

lowpass

Low-pass filter that attenuates high-frequency components.

highpass

High-pass filter that attenuates low-frequency components.

bandpass

Band-pass filter that retains frequencies within a selected range.

bandstop

Band-stop filter that removes frequencies within a selected range.

Examples

file <- system.file("extdata", "EEG.edf", package = "MatchingPursuit")
out <- read_edf_signals(file, resampling = FALSE)
signal <- out$signal
sampling_frequency <- out$sampling_frequency

fc <- design_filters(
  sampling_frequency = sampling_frequency,
  notch = c(49, 51),
  lowpass = 40,
  highpass = 1,
  bandpass = c(0.5, 40),
  bandstop = c(10, 50)
)

print(fc)

signal::freqz(fc$notch, Fs =  sampling_frequency)
signal::freqz(fc$lowpass, Fs =  sampling_frequency)
signal::freqz(fc$highpass, Fs =  sampling_frequency)
signal::freqz(fc$bandpass, Fs =  sampling_frequency)
signal::freqz(fc$bandstop, Fs =  sampling_frequency)

plot(signal[, 1], type = "l", panel.first = grid())

signal_filt <- signal

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

plot(signal_filt[, 1], type = "l", panel.first = grid())

Methods for EDF Objects

Description

Methods for printing, summarizing, and plotting objects of class "edf", together with a print method for objects returned by summary().

Usage

## S3 method for class 'edf'
plot(
  x,
  begin = NULL,
  end = NULL,
  panel_height = NULL,
  rainbow = FALSE,
  bg_colour = "white",
  txt_col = "black",
  zero_line = TRUE,
  main = NULL,
  ...
)

## S3 method for class 'edf'
print(x, ...)

## S3 method for class 'edf'
summary(object, ...)

## S3 method for class 'summary.edf'
print(x, ...)

Arguments

Details

Objects of class "edf" represent multi-channel signals imported from EDF files, together with sampling information, channel names, a time vector, and the source record name.

Value

print.edf() and print.summary.edf() return their input object invisibly.

summary.edf() returns an object of class "summary.edf" containing record information, signal dimensions, sampling information, channel names, and basic descriptive statistics for each channel.

plot.edf() is called for its side effect and returns no value.

See Also

read_edf_signals

Examples

file <- system.file("extdata", "EEG.edf", package = "MatchingPursuit")
x <- read_edf_signals(file, resampling = FALSE)

print(x)
summary(x)

plot(
  x,
  begin = 0,
  end = 10,
  panel_height = NULL,
  rainbow = TRUE,
  bg_colour = "black",
  txt_col = "white",
  zero_line = TRUE,
  main = "EEG signals stored in the EEG.edf file"
)

plot(
  x,
  begin = 0,
  end = 10,
  panel_height = NULL,
  rainbow = FALSE,
  bg_colour = "white",
  txt_col = "black",
  zero_line = TRUE,
  main = "EEG signals stored in the EEG.edf file"
)

Performs bipolar, reference or average EEG montage

Description

An EEG montage refers to the arrangement of EEG electrodes and the way their signals are displayed relative to one another during electroencephalogram interpretation. The same EEG recording may appear very different depending on the montage used. This function implements the three montage methods most commonly used in practice: 1) Bipolar Montage, 2) Referential (Monopolar) Montage, and 3) Average Reference Montage.

Usage

eeg_montage(
  x,
  montage_type = c("average", "reference", "bipolar"),
  ref_channel = NULL,
  bipolar_pairs = NULL
)

Arguments

Details

To check the channel names in the analysed EEG recording, use the read_edf_params() function.

Value

An object of class edf, which is a list with fields:

Examples

file <- system.file("extdata", "EEG.edf", package = "MatchingPursuit")
out <- read_edf_signals(file, resampling = FALSE, from = 0, to = 10)

read_edf_params(file)

# The classical double banana montage.
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")
)

signal_bip_mont <- eeg_montage(out, montage_type = "bipolar", bipolar_pairs = pairs)
signal_ref_mont <- eeg_montage(out, montage_type = "reference", ref_channel = "O1")
signal_avg_mont <- eeg_montage(out, montage_type = "average")

head(signal_bip_mont$signal)
head(signal_ref_mont$signal)
head(signal_avg_mont$signal)

Check whether EMPI is installed

Description

The EMPI program is installed using the empi_install() function and stored in the cache directory. This function checks whether the EMPI program is still available there (users have full access to the cache directory and may remove its contents at any time).

Usage

empi_check()

Value

A character string containing the full path to the EMPI executable if found. If EMPI is not available, invisibly returns NULL and displays a message suggesting installation with empi_install().

See Also

empi_install, empi_locate, empi_execute, plot.mp

Examples

if (interactive()) {
  empi_check()
}

Launches the empi program

Description

Runs the EMPI program for the given data (signal).

Usage

empi_execute(
  signal,
  empi_options = NULL,
  write_to_file = FALSE,
  path = NULL,
  file_name = NULL,
  ...
)

Arguments

Details

The EMPI program (source code and binary files for multiple operating systems) can be downloaded from https://github.com/develancer/empi. Details are presented in the journal paper: Różański, P. T. (2024). empi: GPU-Accelerated Matching Pursuit with Continuous Dictionaries. ACM Transactions on Mathematical Software, Volume 50, Issue 3, Article No. 17, pp. 1-17, doi:10.1145/3674832.

Value

Results of signal decomposition using the MP algorithm. An object of class mp is returned. If write_to_file = TRUE, the results are also written to a SQLite file in the path directory.

See Also

empi_check, empi_install, empi_locate, plot.mp

Examples

if (interactive()) {
file <- system.file("extdata", "sample1.csv", package = "MatchingPursuit")
signal <- read_csv_signals(file)

out_empi <- empi_execute(
  signal = signal
)

# Suppress standard output and standard error
out_empi <- empi_execute(
  signal = signal,
  ignore.stdout = TRUE,
  ignore.stderr = TRUE
)

# The default EMPI options have been modified; see the EMPI README.md for details.
# The '--full-atoms-in-signal' option restricts the decomposition to atoms
# fully contained within the analyzed signal.
# The decomposition results are saved to a SQLite database file.
out_empi <- empi_execute(
  signal = signal,
  empi_options = "-o local --full-atoms-in-signal -i 50 --gabor",
  write_to_file = TRUE,
  path = NULL,
  file_name = "my_decomposition.db"
)

plot(out_empi, freq_divide = 4)
}

Installs the EMPI external program

Description

Downloads the Enhanced Matching Pursuit Implementation (EMPI) external program compatible with the current operating system and stores it in the package cache directory.

Usage

empi_install()

Details

The function detects the operating system (Windows, Linux, macOS arm64), downloads the appropriate archive from the official repository, verifies its integrity using a checksum, and extracts it.

Value

The function downloads the EMPI program in a version compatible with the operating system used (Windows, Linux, MacOS-x64, MacOS-arm64) and stores it in the package cache directory.

See Also

empi_check, empi_locate, empi_execute, plot.mp

Examples

if (interactive()) {
  empi_install()
}

Get required external software localization

Description

Returns Enhanced Matching Pursuit Implementation binary locations for the following operating systems: Windows, Linux, macOS-arm64.

Usage

empi_locate()

Value

A list containing:

• url: URL of the EMPI binary archive,

• fname: archive file name.

See Also

empi_check, empi_install, empi_execute, plot.mp

Examples

sys <- Sys.info()[["sysname"]]
mach <- Sys.info()[["machine"]]

if (sys %in% c("Windows", "Linux") ||
    (sys == "Darwin" && mach == "arm64")) {
  empi_locate()
}

Generate a Gabor atom

Description

Generates a real-valued Gabor atom consisting of a sinusoidal component localized by a Gaussian envelope. Gabor atoms provide simultaneous localization in time and frequency and are commonly used in time-frequency dictionaries for Matching Pursuit decomposition.

Usage

gabor_atom(
  number_of_samples,
  sampling_frequency,
  mean,
  phase,
  sigma,
  frequency,
  normalization = TRUE
)

Arguments

Value

A list containing four numeric vectors of length number_of_samples:

Examples

number_of_samples <- 512
sampling_frequency <- 256.0
mean <- 1
phase <- pi
sigma <- 0.5
frequency <- 5.0
normalization = TRUE

out <- gabor_atom(
  number_of_samples,
  sampling_frequency,
  mean,
  phase,
  sigma,
  frequency,
  normalization
)

# Verify unit-norm normalization
sqrt(sum(out$gabor^2))

plot(out$time, out$gabor, type = "l", xlab = "t", ylab = "gabor", panel.first = grid())

FFT-based computation of projections onto Gabor atoms

Description

Computes complex projection coefficients between one or more signals and a set of Gabor atoms using FFT-based frequency-domain operations.

Usage

gabor_projection_fft(block, signal, sigma_divisor = NULL)

Arguments

Details

For each time position defined in block, the corresponding signal segment is multiplied by a normalized Gaussian envelope and transformed using the Fast Fourier Transform. Only Fourier coefficients corresponding to the frequencies specified in the dictionary are retained.

Atom supports may extend beyond the observed signal boundaries. In this case, samples outside the signal are treated as zero. The complete Gaussian envelope is normalized before boundary truncation; the part overlapping the observed signal is not renormalized.

The Gaussian envelope is constructed over the complete atom support and normalized to unit L2 norm before it is applied to the signal. If the atom support extends before the first signal sample or beyond the last signal sample, only the overlapping signal samples contribute to the projection; values outside the observed signal are implicitly treated as zero.

Consequently, the visible part of a boundary-crossing atom is not renormalized. This preserves the normalization of the complete atom independently of its position relative to the signal boundaries.

Value

A list with two matrices:

Note

This function is primarily intended for internal use by topk_gabor_atoms(), but it is exported to support advanced experiments and methodological testing.

Examples

signal <- as.matrix(rnorm(256))
sampling_frequency <- 256
duration <- 1

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

block <- read_gabor_dict(
  xml_file = xml_file,
  full_atoms_in_signal = FALSE,
  sampling_frequency = sampling_frequency,
  duration = duration,
  verbose = TRUE
)

out <- gabor_projection_fft(block, signal)

pmm <- out$proj_mod_mtx
scm <- out$fft_bin_mtx

head(scm)
head(pmm)

# Projection magnitudes are the moduli of the complex coefficients
head(Mod(scm))

Generate an EMPI-compatible Gabor dictionary

Description

Generates an XML dictionary file containing a set of Gabor atoms.

Usage

generate_xml_dict(N, file, max_window_length = c("3N", "N-1"))

Arguments

Details

The generated dictionary contains multiple Gabor blocks with logarithmically distributed window lengths. The smallest window length is fixed to 17 samples and the maximum window length is determined by max_window_length. It can be set to the largest odd integer not exceeding either 3*N or N-1. The 3*N upper limit is chosen heuristically.

The window lengths are generated on a logarithmic scale and then quantized to obtain a set of practical window sizes. Window lengths are forced to be odd, which provides an exact temporal centre for symmetric Gaussian/Gabor windows. The window shift is estimated as approximately 6 percent of the window length:

windowShift = round(0.06 * windowLen)

The FFT size is selected as the smallest power of two satisfying:

fftSize >= 2 * windowLen

This corresponds to the zero-padding strategy commonly used in MPTK-like Gabor dictionaries.

The dictionary generation procedure is based on the following rules. The rules were derived from an analysis of XML files generated by the EMPI program using the --dictionary-output option (see the read_gabor_dict() function and the corresponding vignette available on CRAN).

1. Minimum window length: 17 samples.

2. Maximum window length: determined by max_window_length.

3. Number of scales:

   K = ceil(log2(N)) + 3

4. Logarithmic spacing of scales:

   L_k = 17 r^k

   where

   r=(L_{max}/17)^{1/(K-1)}

5. Quantization of window lengths depending on their size.

6. Enforcement of odd window lengths.

7. Window shift proportional to window length.

8. FFT size selected as a power of two.

Value

The function writes an XML dictionary file to file. Invisibly returns a data frame containing the generated dictionary parameters:

• windowLen - Gabor window length in samples,

• windowShift - temporal shift between consecutive atoms,

• fftSize - FFT size used for the block.

See Also

read_gabor_dict

Examples

# Generate a dictionary for a 256-sample signal
xml_file <- tempfile(fileext = ".xml")

# Default setting: maximum window length based on 3*N
dict <- generate_xml_dict(
  N = 256,
  file = xml_file
)

dict

# Read the generated Gabor dictionary for a 2-second,
# 128-Hz signal (256 samples)
atoms_dict <- read_gabor_dict(
  xml_file,
  sampling_frequency = 128,
  duration = 2,
  full_atoms_in_signal = FALSE,
  verbose = TRUE
)

head(atoms_dict)
tail(atoms_dict)

# Alternative setting: limit the maximum window length to N-1
dict <- generate_xml_dict(
  N = 256,
  file = xml_file,
  max_window_length = "N-1"
)

dict

Methods for Matching Pursuit Objects

Description

Methods for printing, summarizing, and plotting objects of class "mp", together with a print method for objects returned by summary().

Usage

## S3 method for class 'mp'
plot(
  x,
  channel = 1,
  mode = "sqrt",
  freq_divide = NULL,
  increase_factor = 8,
  shortening_factor_x = 2,
  shortening_factor_y = 2,
  atom_centers = "crosses",
  display_grid = FALSE,
  color = "white",
  palette = "my custom palette",
  plot_signals = TRUE,
  verbose = FALSE,
  ...
)

## S3 method for class 'mp'
print(x, ...)

## S3 method for class 'mp'
summary(object, ...)

## S3 method for class 'summary.mp'
print(x, ...)

Arguments

Details

Objects of class "mp" represent the result of a Matching Pursuit or Orthogonal Matching Pursuit decomposition. They contain the original signal, its reconstruction, parameters of the selected atoms, individual selected atom waveforms, the corresponding time vector, and sampling frequency.

The plotting method displays a time-frequency (T-F) map to visualize the decomposition. It is a wrapper around tf_map() with out_mode = "plot".

Value

print.mp() and print.summary.mp() return their input object invisibly.

summary.mp() returns an object of class "summary.mp" containing signal dimensions, sampling information, the total number of selected atoms, the number of selected atoms per channel, #' signal, residual, and explained energy information, and ranges of selected atom parameters.

plot.mp() is called for its side effect and returns no value.

See Also

tf_map, empi_execute, mp_omp_execute

Examples

## Not run: 
file <- system.file("extdata", "sample1.csv", package = "MatchingPursuit")
signal <- read_csv_signals(file, col_names = "ch1")

# Execute the MP algorithm.
out_empi <- empi_execute(signal = signal)

# Print and summarize the decomposition.
print(out_empi)
summary(out_empi)

# Plot a time-frequency map based on MP atoms.
plot(out_empi)

## End(Not run)

Implements the classical Matching Pursuit (MP) algorithm

Description

Computes a sparse representation of a signal using the classical Matching Pursuit (MP) algorithm and a dictionary of atoms.

Usage

mp_core(
  dictionary,
  signal,
  channel = NULL,
  n_nonzero_coefs = NULL,
  tol = NULL,
  verbose = FALSE
)

Arguments

Details

This is a native R implementation of the classical MP algorithm supporting arbitrary matrix-based dictionaries. It provides direct access to the decomposition procedure and can be used independently of the Gabor-specific workflow. For high-performance Matching Pursuit decomposition using the external EMPI backend, see empi_locate(), empi_install(), empi_check(), and empi_execute().

Dictionary atoms are normalized internally to unit L2 norm. For classical MP, the reduction in squared residual norm at iteration k is therefore theoretically equal to the squared MP coefficient. The former definition is used explicitly to maintain a consistent interpretation of energy across MP and OMP decompositions.

Value

A list containing the result of the Matching Pursuit decomposition with the following elements:

See Also

read_gabor_dict, topk_gabor_atoms, mp_omp_execute

Examples

dictionary <- matrix(
c(
  1.0,  0.9,  0.1,  1.0, -0.2,  0.3,  0.7, -0.5,  1.2,  0.4,
  0.2,  1.0,  0.8, -0.3,  1.0, -0.6,  0.5,  0.9, -0.1,  0.8,
  0.0,  0.1,  1.0,  0.5,  0.7,  1.1, -0.4,  0.2,  0.6, -0.7,
  0.9, -0.2,  0.4,  1.3,  0.1,  0.0,  0.8, -0.9,  0.5,  1.0,
 -0.3,  0.6,  1.1, -0.4,  0.2,  0.7, -0.8,  1.0,  0.3,  0.9),
nrow = 5, byrow = TRUE
)

signal <- matrix(
c(
  4, 3, 5, 2,
  2, 1, 2, 3,
  3, 2, 4, 1,
  5, 4, 3, 2,
  1, 3, 2, 4),
nrow = 5, byrow = TRUE
)

fit <- mp_core(
  dictionary = dictionary,
  signal = signal,
  channel = 1,
  n_nonzero_coefs = 3,
  verbose = TRUE
)

fit$coefs
# [1] 6.274348 2.535423 1.870568

fit$support
# [1] 9 5 7

fit$relative_residual_energy
# [1] 1.0000000 0.2842283 0.1673489 0.1037303

# More realistic example, see mp_omp_execute() examples.

Matching Pursuit (MP) or Orthogonal Matching Pursuit (OMP) decomposition for multi-channel signals

Description

Performs sparse signal decomposition using either the Matching Pursuit (MP) or Orthogonal Matching Pursuit (OMP) algorithm, as specified by the mode parameter. The decomposition is performed independently for each signal channel using a dictionary of candidate atoms selected by topk_gabor_atoms().

Usage

mp_omp_execute(
  signal,
  mode = NULL,
  dictionary = NULL,
  topk = NULL,
  full_atoms_in_signal = FALSE,
  n_nonzero_coefs = NULL,
  tol = NULL,
  verbose = FALSE
)

Arguments

Details

The returned object is of class "mp" and can be visualized using plot() and tf_map().

The XML dictionary specification is first processed by read_gabor_dict(), after which topk_gabor_atoms() selects a channel-specific subset of candidate Gabor atoms. The selected dictionary is then passed to mp_core() or omp_core() independently for each signal channel.

The results from all channels are combined into an object of class "mp", which can be visualized using plot() and tf_map().

Value

An object of class "mp" containing:

The atoms data frame contains:

• channel_id — signal channel identifier,

• atom_number — atom index within the channel,

• energy — atom energy contribution,

• envelope — envelope type,

• frequency — atom frequency (Hz),

• phase — atom phase (radians),

• scale — atom scale (seconds),

• position — atom centre position (seconds).

See Also

read_gabor_dict, topk_gabor_atoms, generate_xml_dict, omp_core, mp_core

Examples

# +-------------------------------------------------------------+
# | Read signal                                                 |
# +-------------------------------------------------------------+
file <- system.file(
  "extdata",
  "sample1.csv",
  package = "MatchingPursuit"
)

signal <- read_csv_signals(
  file,
  col_names_in_csv = FALSE
)

# +-------------------------------------------------------------+
# | Run Matching Pursuit (MP-R backend)                         |
# +-------------------------------------------------------------+
# set  "mode = omp" to run Orthogonal Matching Pursuit (OMP-R backend)
#
# topk is set to a relatively small value to reduce computation time.
# If this parameter is omitted, the function sets topk by default to 10%
# of the total number of atoms in the dictionary.
#
# In practical applications, n_nonzero_coefs is typically set to a larger
# value, e.g. 50.
#
fit_mp <- mp_omp_execute(
  mode = "mp",
  signal = signal,
  n_nonzero_coefs = 25,
  topk = 5000,
  verbose = TRUE
)

plot(fit_mp, freq_divide = 4)

# ### NOTE ###
# Additional examples are provided below for illustration.
# They are commented out because they may take longer to run.

# The '--full-atoms-in-signal' option restricts the
# decomposition to atoms fully contained within the analyzed
# signal. Compare the two time-frequency maps obtained with
# and without this option.

# fit_mp <- mp_omp_execute(
#   mode = "mp",
#   signal = signal,
#   full_atoms_in_signal = TRUE,
#   n_nonzero_coefs = 50,
#   verbose = TRUE
# )
# plot(fit_mp, freq_divide = 4)

# +-------------------------------------------------------------+
# | Run Orthogonal Matching Pursuit (OMP-R backend)             |
# +-------------------------------------------------------------+
# fit_omp <- mp_omp_execute(
#   mode = "omp",
#   signal = signal,
#   n_nonzero_coefs = 50,
#   verbose = TRUE
# )

# plot(fit_omp, freq_divide = 4)

# +-------------------------------------------------------------+
# | Use an external XML dictionary specification                |
# +-------------------------------------------------------------+
# xml_file <- system.file(
#   "extdata",
#   "sample1.xml",
#    package = "MatchingPursuit"
# )

# fit_mp_xml <- mp_omp_execute(
#   mode = "mp",
#   signal = signal,
#   dictionary = xml_file,
#   n_nonzero_coefs = 50,
#   verbose = TRUE
# )

# plot(fit_mp_xml, freq_divide = 4)

Implements Orthogonal Matching Pursuit (OMP) algorithm

Description

Performs Orthogonal Matching Pursuit (OMP) to obtain a sparse representation of a signal using a dictionary of candidate atoms.

Usage

omp_core(
  dictionary,
  signal,
  channel = NULL,
  n_nonzero_coefs = NULL,
  tol = NULL,
  verbose = FALSE
)

Arguments

Details

Unlike classical Matching Pursuit, OMP recomputes the coefficients of all previously selected atoms at each iteration by solving a least-squares problem. This makes the residual orthogonal to the subspace spanned by the selected atoms and generally provides a more accurate approximation for a given number of atoms.

The least-squares problem is solved efficiently using incremental Cholesky factorization.

At each iteration, the atom with the largest absolute correlation with the current residual is selected. All coefficients associated with the selected atoms are then recomputed simultaneously by least squares, and the residual is updated.

For OMP, the quantity stored in energy is not calculated from the final coefficient of an individual atom. Selected atoms can be strongly correlated, and therefore quantities such as coefs^2 * colSums(selected_atoms^2) are not additive contributions to the energy of the reconstruction.

Instead, the energy associated with the atom selected at iteration k is defined as the reduction in residual energy:

E_k = \|r_{k-1}\|_2^2 - \|r_k\|_2^2.

Thus, energy[k] quantifies the reduction in residual energy associated with the k-th OMP iteration, after adding a new atom and re-estimating all active coefficients.

\sum_k E_k = \|x\|_2^2 - \|r_K\|_2^2.

Value

A list containing the result of the Orthogonal Matching Pursuit decomposition with the following elements:

See Also

topk_gabor_atoms, mp_core, mp_omp_execute, tf_map

Examples

dictionary <- matrix(
c(
  1.0,  0.9,  0.1,  1.0, -0.2,  0.3,  0.7, -0.5,  1.2,  0.4,
  0.2,  1.0,  0.8, -0.3,  1.0, -0.6,  0.5,  0.9, -0.1,  0.8,
  0.0,  0.1,  1.0,  0.5,  0.7,  1.1, -0.4,  0.2,  0.6, -0.7,
  0.9, -0.2,  0.4,  1.3,  0.1,  0.0,  0.8, -0.9,  0.5,  1.0,
 -0.3,  0.6,  1.1, -0.4,  0.2,  0.7, -0.8,  1.0,  0.3,  0.9),
nrow = 5, byrow = TRUE
)

signal <- matrix(
c(
  4, 3, 5, 2,
  2, 1, 2, 3,
  3, 2, 4, 1,
  5, 4, 3, 2,
  1, 3, 2, 4),
nrow = 5, byrow = TRUE
)

fit <- omp_core(
  dictionary = dictionary,
  signal = signal,
  channel = 1,
  n_nonzero_coefs = 3,
  verbose = TRUE
)

fit$coefs
# [1] 5.282278 2.637693 2.195920

fit$support
# [1] 9 5 7

fit$relative_residual_energy
# [1] 1.00000000 0.28422833 0.16609350 0.08795047

# For a complete Gabor decomposition workflow, see mp_omp_execute().

Reference implementation of Orthogonal Matching Pursuit (OMP)

Description

A straightforward reference implementation of the Orthogonal Matching Pursuit (OMP) algorithm that closely follows its mathematical formulation. The function is intended for reference and illustrative purposes and does not use computational optimizations. The least-squares problem is solved explicitly using the normal-equation formula.

Usage

omp_reference(
  dictionary,
  signal,
  n_nonzero_coefs = NULL,
  tol = NULL,
  verbose = FALSE
)

Arguments

Details

This implementation is intended for small illustrative examples rather than large-scale computations. In particular, when verbose = TRUE, detailed information is printed at every iteration, including residual and reconstruction vectors. For large signals or a large number of iterations, this may produce a substantial amount of console output.

This implementation operates on a single signal represented by a numeric vector. The signal argument must therefore be a numeric vector, not a matrix containing multiple signals or channels.

Because the least-squares problem is solved explicitly through the normal equations, the implementation may fail for singular or nearly singular selected subdictionaries and is not intended for numerically demanding applications.

Dictionary atoms are internally normalized to unit Euclidean norm before decomposition. The returned coefficients are provided both for the normalized working dictionary and for the original dictionary scaling.

Value

A list containing:

See Also

omp_core, mp_omp_execute

Examples

dictionary <- matrix(
c(
  1.0,  0.9,  0.1,  1.0, -0.2,  0.3,  0.7, -0.5,  1.2,  0.4,
  0.2,  1.0,  0.8, -0.3,  1.0, -0.6,  0.5,  0.9, -0.1,  0.8,
  0.0,  0.1,  1.0,  0.5,  0.7,  1.1, -0.4,  0.2,  0.6, -0.7,
  0.9, -0.2,  0.4,  1.3,  0.1,  0.0,  0.8, -0.9,  0.5,  1.0,
  -0.3,  0.6,  1.1, -0.4,  0.2,  0.7, -0.8,  1.0,  0.3,  0.9),
nrow = 5, byrow = TRUE
)

signal <- c(4, 2, 3, 5, 1)

out <- omp_reference(
  dictionary = dictionary,
  signal = signal,
  n_nonzero_coefs = 3,
  verbose = TRUE
)

out

Read atom parameters from a SQLite database

Description

Reads the atom parameters from a SQLite database produced by empi_execute().

Usage

read_atom_params(db_file)

Arguments

Value

A data frame containing the atom parameters stored in the database:

Examples

# Example database containing data from 18 channels
file <- system.file("extdata", "EEG_filter_resample_montage.db", package = "MatchingPursuit")
out <- read_atom_params(file)
out[which(out$channel_id == 1), ]
out[which(out$channel_id == 18), ]

Reads and validates a CSV file structure

Description

Reads and validates a CSV file structure

Usage

read_csv_signals(file, col_names = NULL, col_names_in_csv = FALSE)

Arguments

Value

A list containing:

signal

Data frame containing all signals (rows = samples, columns = channels).

sampling_frequency

Sampling frequency.

time

Time vector corresponding to signal samples.

Examples

file <- system.file("extdata", "sample1.csv", package = "MatchingPursuit")

# The first line of the file must contain two numbers:
# a) the sampling frequency in Hz
# b) the signal length in seconds
out <- read.csv(file, header = FALSE)
head(out)

signal <- read_csv_signals(file, col_names = "signal_1")
head(signal$signal)
signal$sampling_frequency
head(signal$time)
tail(signal$time)

file <- system.file("extdata", "sample2.csv", package = "MatchingPursuit")
signal <- read_csv_signals(file, col_names = c("signal_1"))
head(signal$signal)
signal$sampling_frequency

# Now, the csv file contains signal names in the second line
file <- system.file("extdata", "sample3.csv", package = "MatchingPursuit")
signal <- read_csv_signals(file, col_names_in_csv = TRUE)
head(signal$signal)
signal$sampling_frequency

Reads a selected EDF or EDF+ file and returns signal parameters

Description

Reads a selected EDF or EDF+ file and returns basic signal parameters (channel names, sampling frequency of each channel, number of samples per channel, and signal duration in seconds). Additional information stored in EDF+ files (such as interrupted recordings or time-stamped annotations) is not used by the package and is therefore not read.

Usage

read_edf_params(file)

Arguments

Value

A data frame containing the basic parameters of the EDF / EDF+ file:

Examples

file <- system.file("extdata", "EEG.edf", package = "MatchingPursuit")
read_edf_params(file)

Reads a selected EDF or EDF+ file and returns signal data

Description

The function reads a selected EDF or EDF+ file. Optionally, resampling can be performed (upsampling or downsampling).

Usage

read_edf_signals(
  file,
  resampling = FALSE,
  sf_new = NULL,
  from = NULL,
  to = NULL,
  verbose = FALSE
)

Arguments

Details

If resampling = TRUE, signals are resampled according to the target frequency specified by f.new. Since the EDF standard allows different sampling rates per channel, some channels may be upsampled while others are downsampled. The function does not support independent resampling of individual channels.

Value

An object of class edf, which is a list with fields:

Examples

# Read EDF signals without resampling
file <- system.file("extdata", "EEG.edf", package = "MatchingPursuit")
out1  <- read_edf_signals(file, resampling = FALSE)

lapply(out1, class)
out1$sampling_frequency

# Read EDF signals and resample them to 128 Hz
out2 <- read_edf_signals(file, resampling = TRUE, sf_new = 128, verbose = TRUE)

lapply(out2, class)
out2$sampling_frequency

Read EMPI decomposition results from a SQLite database

Description

Reads data from a SQLite file (.db) created by the Matching Pursuit algorithm. The reconstructed signal(s) and Gabor function(s) are also returned.

Usage

read_empi_db(db_file)

Arguments

Value

An object of class "mp" containing:

atoms

A data frame describing the selected atoms.

signal

Matrix containing the original signal(s).

reconstruction

Matrix containing the reconstructed signal(s).

selected_atoms

List of matrices containing selected atoms for each channel.

time

Time vector corresponding to signal samples.

sampling_frequency

Sampling frequency.

Examples

file <- system.file("extdata", "EEG_filter_resample_montage.db", package = "MatchingPursuit")
out <- read_empi_db(file)

n_channels <- ncol(out$signal)
signal <- out$signal
reconstruction <- out$reconstruction
t <- out$time
sampling_frequency <- out$sampling_frequency

old.par <- par("mfrow", "pty", "mai")

par(mfrow = c(2, 1))
par(pty = "m")
par(mai = c(0.9, 0.5, 0.3, 0.4))

plot(
  signal[,1], type = "l", col = "blue",
  main = paste("channel: ", 1, " / " , n_channels, " (original signal)",  sep = ""),
  xaxt = "n", ylab = "", xlab = "time [sec]"
)

len <- length(signal[, 1])
lab <- seq(t[1], t[len] + 1 / sampling_frequency, length.out = 11)
axis(side = 1, las = 1, cex.axis = 0.9, at = seq(0, len, length.out = 11), labels = lab)

plot(
  reconstruction[,1], type = "l", col = "blue",
  main = paste("channel: ", 1, " / " , n_channels, " (reconstructed signal)",  sep = ""),
  xaxt = "n", ylab = "", xlab = "time [sec]"
)

axis(side = 1, las = 1, cex.axis = 0.9, at = seq(0, len, length.out = 11), labels = lab)

par(old.par)

Read a Gabor dictionary from an XML file

Description

The function parses an XML file describing a multiscale Gabor dictionary.

Usage

read_gabor_dict(
  xml_file,
  sampling_frequency,
  duration,
  verbose = FALSE,
  full_atoms_in_signal = FALSE
)

Arguments

Details

Each <block> in the XML file defines a time-frequency scale of atoms using three parameters:

• windowLen — length of the analysis window (in samples),

• windowShift — step size between consecutive windows,

• fftSize — FFT size defining frequency resolution.

The function assumes an XML structure containing param nodes with name and value attributes. An example XML file is shown below. For simplicity, the example contains only one block; in practice, dictionary files usually contain multiple blocks.

<?xml version="1.0" encoding="ISO-8859-1"?>
<dict>
 <block>
   <param name="windowLen" value="31"/>
   <param name="windowShift" value="2"/>
   <param name="fftSize" value="64"/>
 </block>
</dict>

Each block generates a grid of atoms over time and frequency bins, forming a multiresolution Gabor dictionary. Smaller windows provide better time resolution, while larger windows improve frequency resolution.

The treatment of atoms near signal boundaries is controlled by full_atoms_in_signal.

If TRUE, only atoms fully contained within the signal support are generated. In this case, atom start positions satisfy

0 \leq t \leq N - L,

where t is the atom start position, N is the signal length, and L is the window length. If the signal is shorter than the window length, no time positions are generated for that block.

If FALSE, atom centres are allowed at positions throughout the signal, and the support of an atom may extend beyond the signal boundaries.

Value

A matrix where each row describes a Gabor atom with the following columns:

Usage in sparse decomposition workflow

The output of read_gabor_dict() is a low-level description of the Gabor time-frequency grid. It serves as input to topk_gabor_atoms(), which:

• evaluates complex Gabor atoms,

• computes phase-invariant projections onto the signal,

• selects the best topk atoms for each channel,

• constructs real-valued atom representations using optimal phases.

The resulting "topk" object contains channel-specific atom matrices and associated metadata. Individual atom matrices can subsequently be passed to mp_core() or omp_core() for sparse decomposition. The higher-level mp_omp_execute() function performs these preparation steps internally.

EMPI compatibility

XML dictionary definitions exported by EMPI can be read directly by this function. Additional XML elements not used by read_gabor_dict() are ignored. The argument full_atoms_in_signal controls the boundary convention corresponding to the EMPI --full-atoms-in-signal option.

For full details on the EMPI options and their behavior, see the EMPI documentation in README.md.

See Also

topk_gabor_atoms, mp_omp_execute, omp_core, mp_core, generate_xml_dict

Examples

# +-------------------------------------------------------------+
# | Step 1: Read signal                                         |
# +-------------------------------------------------------------+
sig_file <- system.file(
  "extdata",
  "sample3.csv",
  package = "MatchingPursuit"
)

sample3 <- read_csv_signals(
  sig_file,
  col_names_in_csv = TRUE
)

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

# +-------------------------------------------------------------+
# | Step 2: Read dictionary                                     |
# +-------------------------------------------------------------+
xml_file <- system.file(
  "extdata",
  "sample3.xml",
  package = "MatchingPursuit"
)

# +-------------------------------------------------------------+
# | Step 3: Compare boundary conventions                        |
# +-------------------------------------------------------------+
# Generate only atoms whose complete support lies within
# the signal boundaries.
atoms_full <- read_gabor_dict(
  xml_file = xml_file,
  sampling_frequency = sampling_frequency,
  duration = duration,
  full_atoms_in_signal = TRUE,
  verbose = TRUE
)

# Allow atom support to extend beyond the signal boundaries.
# Atom centres still remain within the signal.
atoms_overstep <- read_gabor_dict(
  xml_file = xml_file,
  sampling_frequency = sampling_frequency,
  duration = duration,
  full_atoms_in_signal = FALSE,
  verbose = TRUE
)

# Allowing boundary overstep increases the number of dictionary atoms.
nrow(atoms_full)
nrow(atoms_overstep)

# With full_atoms_in_signal = TRUE, atom start positions
# are always non-negative.
range(atoms_full[, "time_sample"])

# With full_atoms_in_signal = FALSE, atoms centred near the beginning
# of the signal may have negative start positions.
range(atoms_overstep[, "time_sample"])

Reads WFDB-compatible signal and header files

Description

WFDB (WaveForm DataBase) is a standard file format for storing, reading, and analyzing physiological time-series signals. It is widely used for signals such as ECG, EEG, blood pressure, respiration, and other biomedical waveforms. It is the file format used by the PhysioNet project and is commonly used in research datasets.

Usage

read_wfdb_signals(file)

Arguments

Details

A WFDB record typically consists of two main files: .dat - binary signal samples (waveform values), and .hea - a header file describing how to interpret the data. In some cases, additional annotation files such as .atr may be present, containing beat labels or rhythm annotations.

Value

An object of class wfdb. The returned value is a list containing:

signal

Matrix of signals stored in the WFDB file.

sampling_frequency

Sampling frequency.

time

Time vector corresponding to signal samples.

lead_names

Names of the WFDB leads (channels).

record_name

Name of the file.

Note

The function EGM::read_wfdb() from version 0.2.0 of the EGM package does not support multi-frequency signals. Consequently, records containing different numbers of samples per frame, as indicated by the 16x2, 16x4, and 16x1 specifications below, cannot be read correctly.

⁠ multi_freq_test 3 100 1000 multi_freq_test.dat 16x2 200.0(0)/mV 16 0 0 258 0 ECG multi_freq_test.dat 16x4 400.0(0)/mmHg 16 0 400 57824 0 ABP multi_freq_test.dat 16x1 100.0(0)/pm 16 0 0 18204 0 RESP ⁠

Examples

# ECG data comes from https://physionet.org/content/ptb-xl/1.0.3/
file <- system.file("extdata", "00001_lr.hea", package = "MatchingPursuit")

out <- read_wfdb_signals(file)
head(out$signal)
out$sampling_frequency
out$lead_names

plot(out, begin = 0, end = 10, panel_height = 1.5)

Resample a signal (upsampling or downsampling)

Description

Resamples one or more dimensional numeric signals using signal::resample().

Usage

resample_signal(signal, p, q, d = 5)

Arguments

Details

The new sampling frequency is determined by the ratio p/q:

f_{\mathrm{new}} = f_{\mathrm{old}} \frac{p}{q}.

For matrices and data frames, resampling is performed independently for each column. Rows are interpreted as time samples and columns as individual signal channels.

The function uses resample internally. The resampling process includes interpolation, low-pass filtering, and decimation.

Value

A numeric vector, matrix, or data frame containing the resampled signal. The output type matches the input type. Column names are preserved for matrices and data frames.

Examples

# Numeric vector
signal <- sin(2 * pi * 5 * seq(0, 1, length.out = 400))
signal_resampled <- resample_signal(signal, p = 1, q = 4)

old.par <- par("mfrow", "mai")
par(mfrow = c(2, 1))
par(mai = c(0.9, 0.5, 0.3, 0.4))

plot(signal, type = "o")
plot(signal_resampled, type = "o")

par(old.par)

# Numeric matrix: samples in rows, channels in columns (256Hz, 10sec., 5 channels)
signal <- matrix(rnorm(2560 * 5), nrow = 2560, ncol = 5)
colnames(signal) <- paste0("channel_", seq_len(ncol(signal)))

# Resample to 64Hz
signal_64 <- resample_signal(signal, p = 1, q = 4)
dim(signal_64)

# Data frame
signal_df <- as.data.frame(signal)
signal_df_64 <- resample_signal(signal_df, p = 1, q = 4)
names(signal_df_64)

Methods for Signal Objects

Description

Methods for printing, summarizing, and plotting objects of class "sig", together with a print method for objects returned by summary().

Usage

## S3 method for class 'sig'
print(x, ...)

## S3 method for class 'sig'
summary(object, ...)

## S3 method for class 'summary.sig'
print(x, ...)

## S3 method for class 'sig'
plot(
  x,
  begin = NULL,
  end = NULL,
  panel_height = NULL,
  zero_line = TRUE,
  main = NULL,
  mar = c(4, 7, 2, 1),
  col = "black",
  lwd = 1,
  lty = 1,
  ...
)

Arguments

Details

Objects of class "sig" represent a single- or multi-channel signal together with its sampling frequency and corresponding time vector.

Value

print.sig() and print.summary.sig() return their input object invisibly.

summary.sig() returns an object of class "summary.sig" containing signal dimensions, sampling information, channel names, and basic descriptive statistics for each channel.

plot.sig() is called for its side effect and returns x invisibly.

See Also

read_csv_signals, as_sig

Examples

file <- system.file("extdata", "sample3.csv", package = "MatchingPursuit")

x <- read_csv_signals(
  file,
  col_names_in_csv = TRUE
)

print(x)
summary(x)

plot(x)
plot(x, mar = c(4, 12, 2, 1))
plot(x, mar = c(4, 12, 2, 1), col = "blue")
plot(x, col = "red", lwd = 2, lty = 4)
plot(x, begin = 0, end = 1, main = "Signal")

Convert multichannel signals to binary format

Description

Converts a numeric matrix or data frame containing one or more signals to the binary format required by EMPI. Rows correspond to samples and columns to channels. Values are stored as 4-byte floating-point numbers using little-endian byte order.

For multichannel signals, samples are written in time order, with all channel values for a given time point stored consecutively: first all channels at t = 0, then all channels at t = \Delta t, and so on.

Usage

signal_to_bin(data, write_to_file = FALSE, path = NULL, file_name = NULL)

Arguments

Value

A raw vector containing the binary representation of the signal. If write_to_file = TRUE, a .bin file is additionally created.

Note

The .bin files generated by this function are not intended for direct user manipulation. They are used internally by empi_execute(). The external program Enhanced Matching Pursuit Implementation (EMPI) requires binary input data. This conversion utility may also be useful for users who wish to run EMPI outside of the R environment.

Examples

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

signal_bin <- signal_to_bin(data = out$signal, write_to_file = FALSE)

# We have 3 channels. The first 4 time points.
head(out$signal, 4)

# The same elements of the signal in binary (floats are stored in 4 bytes).
head(signal_bin, 48)

# After decoding to numeric.
# Of course we get the same values as in out$signal.
readBin(signal_bin[1:4], what = "numeric", size = 4, endian = "little")
readBin(signal_bin[5:8], what = "numeric", size = 4, endian = "little")
readBin(signal_bin[41:44], what = "numeric", size = 4, endian = "little")
readBin(signal_bin[45:48], what = "numeric", size = 4, endian = "little")

Creates a time-frequency map using atoms from the Matching Pursuit algorithm

Description

Creates a time-frequency map using atoms from the Matching Pursuit algorithm. The resulting map can be: 1) displayed on the screen, 2) saved as a .png file, or 3) saved as an .RData object.

Usage

tf_map(
  x = NULL,
  channel,
  mode = "sqrt",
  freq_divide = NULL,
  increase_factor = 4,
  shortening_factor_x = 2,
  shortening_factor_y = 2,
  atom_centers = "crosses",
  display_grid = FALSE,
  color = "white",
  palette = "my custom palette",
  reverse_palette = TRUE,
  out_mode = "plot",
  path = NULL,
  file_name = NULL,
  size = c(512, 512),
  draw_ellipses = FALSE,
  plot_signals = TRUE,
  write_atoms = FALSE,
  verbose = FALSE
)

Arguments

Details

The function also computes basic reconstruction-quality measures for the selected channel, including signal energy, residual energy, and normalized reconstruction error (NRE).

Value

Depending on the out_mode parameter, the function:

• displays the time-frequency map on the screen

• saves the time-frequency map as a .png file

• saves the time-frequency map as a .RData file

Regardless of the output mode, the function also returns:

Examples

if (interactive()) {
file <- system.file("extdata", "sample1.csv", package = "MatchingPursuit")
signal <- read_csv_signals(file)

sample1_empi_out <- empi_execute (
  signal = signal,
  empi_options = "-o local --gabor -i 25",
)

# 'freq_divide' is set arbitrarily
out <- tf_map(
  x = sample1_empi_out,
  channel = 1,
  mode = "sqrt",
  freq_divide = 4,
  increase_factor= 4,
  atom_centers = "crosses",
  out_mode = "plot",
)

# 'freq_divide' is determined based on the atom with the highest frequency
out <- tf_map(
  x = sample1_empi_out,
  channel = 1,
  mode = "sqrt",
  increase_factor= 4,
  atom_centers = "numbers",
  out_mode = "plot",
)
}

Methods for Top-k Objects

Description

Print and summarize objects returned by topk_gabor_atoms().

Usage

## S3 method for class 'topk'
print(x, ...)

## S3 method for class 'topk'
summary(object, ...)

## S3 method for class 'summary.topk'
print(x, ...)

Arguments

Value

print.topk() returns x invisibly.

summary.topk() returns an object of class "summary.topk" containing basic information about the selected atoms and ranges of their Gabor parameters.

print.summary.topk() returns x invisibly.

See Also

topk_gabor_atoms

Examples

signal <- read_csv_signals(system.file("extdata", "sample1.csv", package = "MatchingPursuit"))
xml_file <- system.file("extdata", "sample1.xml", package = "MatchingPursuit")

dictionary <- read_gabor_dict(
  xml_file = xml_file,
  sampling_frequency = signal$sampling_frequency,
  duration = max(signal$time)
)

out_topk <- topk_gabor_atoms(
  atoms_dict = dictionary,
  signal = signal,
  topk = 100
)

print(out_topk)
summary(out_topk)

Select the most relevant Gabor atoms using phase-invariant similarity

Description

Constructs a sparse, signal-dependent Gabor dictionary by selecting the atoms with the largest phase-invariant projections onto the input signal.

Usage

topk_gabor_atoms(
  atoms_dict,
  signal,
  topk = NULL,
  sigma_divisor = NULL,
  verbose = FALSE
)

Arguments

Details

In the first step, complex projection coefficients are computed for all candidate atoms using gabor_projection_fft(), and atoms are ranked according to the magnitudes of these coefficients. In the second step, the top-ranked atoms are reconstructed as real-valued time-domain atoms using their optimal phases.

Candidate atoms are ranked separately for each signal channel according to the magnitudes of their complex projection coefficients. This makes the ranking independent of the phase of the corresponding real-valued Gabor atom. The phase of each selected atom is subsequently obtained from the complex coefficient and used to construct its real-valued representation.

Atom supports may extend beyond the observed signal boundaries when the input dictionary was generated with full_atoms_in_signal = FALSE in read_gabor_dict(). In this case, signal values outside the observed interval are implicitly treated as zero.

The complete real-valued Gabor atom is normalized before boundary truncation. If only part of an atom overlaps the observed signal, the retained fragment is not renormalized. Consequently, boundary-crossing atoms stored in the returned object may have an L2 norm smaller than one, whereas atoms whose complete support lies within the signal have unit L2 norm.

The returned "topk" object contains the channel-specific atom matrices that can be passed to mp_core() or omp_core() for sparse signal decomposition.

Value

An object of class "topk", containing:

See Also

read_gabor_dict, gabor_projection_fft, mp_core, omp_core, mp_omp_execute

Examples

# +-------------------------------------------------------------+
# | Step 1: Read signal                                         |
# +-------------------------------------------------------------+
sig_file <- system.file(
  "extdata",
  "sample3.csv",
  package = "MatchingPursuit"
)

signal <- read_csv_signals(
  sig_file,
  col_names_in_csv = TRUE
)

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

# +-------------------------------------------------------------+
# | Step 2: Read dictionary                                     |
# +-------------------------------------------------------------+
xml_file <- system.file(
  "extdata",
  "sample3.xml",
  package = "MatchingPursuit"
)

atoms_dict <- read_gabor_dict(
  xml_file = xml_file,
  sampling_frequency = sampling_frequency,
  duration = duration,
  verbose = FALSE,
  full_atoms_in_signal = TRUE
)

# +-------------------------------------------------------------+
# | Step 3: Select top-k atoms most similar to the signal       |
# +-------------------------------------------------------------+
out_topk <- topk_gabor_atoms(
  atoms_dict = atoms_dict,
  signal = signal,
  topk = 500,
  verbose = TRUE
)

class(out_topk)
dim(out_topk$atoms[[1]])
head(out_topk$frequency[, 1])

Methods for WFDB Objects

Description

Methods for printing, summarizing, and plotting objects of class "wfdb", together with a print method for objects returned by summary().

Usage

## S3 method for class 'wfdb'
plot(
  x,
  begin = NULL,
  end = NULL,
  panel_height = 3,
  small_squares = TRUE,
  zero_line = FALSE,
  ...
)

## S3 method for class 'wfdb'
print(x, ...)

## S3 method for class 'wfdb'
summary(object, ...)

## S3 method for class 'summary.wfdb'
print(x, ...)

Arguments

Details

Objects of class "wfdb" represent multi-channel physiological signals imported from WFDB records, together with their sampling frequency, time vector, lead names, and record name.

WFDB (WaveForm DataBase) is a widely used format and software framework for storing, reading, and analyzing physiological time-series signals. It is widely used for signals such as ECG, EEG, blood pressure, respiration, and other biomedical waveforms. It is the file format used by the PhysioNet project and is commonly used in research datasets.

A WFDB record typically consists of two main files: .dat, containing binary signal samples, and .hea, a header file describing how to interpret the data. Additional annotation files, such as .atr, may also be present and may contain beat labels or rhythm annotations.

The plotting method is designed primarily for ECG signals and displays individual leads in a layout resembling standard ECG paper. The small grid corresponds to 0.04 s by 0.1 mV and the large grid to 0.20 s by 0.5 mV.

Value

print.wfdb() and print.summary.wfdb() return their input object invisibly.

summary.wfdb() returns an object of class "summary.wfdb" containing record information, signal dimensions, sampling information, lead names, and basic descriptive statistics for each lead.

plot.wfdb() is called for its side effect and returns no value.

See Also

read_wfdb_signals

Examples

# ECG data from the PTB-XL database
file <- system.file("extdata", "00001_lr.hea", package = "MatchingPursuit")
x <- read_wfdb_signals(file)

print(x)
summary(x)

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