Package 'PhysioHDEMG'

Decompose HD-sEMG into motor-unit discharge patterns

Description

Decomposes a high-density surface EMG electrode-grid recording into its constituent motor-unit pulse trains by convolutive blind source separation: the channels are extended by delayed copies and whitened, each candidate source is seeded by a fixed-point contrast (kurtosis) maximisation and then refined by convolution-kernel compensation (the separation vector is iteratively set to the average whitened observation at the estimated discharge times), and discharge times are read off by peak detection with a two-class (spike / noise) split. Sources are accepted on a pulse-to-noise ratio and silhouette quality gate, and duplicate units are merged.

Usage

hdEMGDecompose(
  x,
  n_units = 8,
  extension_factor = 16,
  sampling_rate = NULL,
  assay_name = NULL,
  refractory_sec = 0.02,
  pnr_threshold = 26,
  sil_threshold = 0.9,
  max_iter = 100,
  tol = 1e-04,
  bandpass = NULL,
  seed = NULL
)

Arguments

x	HD-sEMG data: an ⁠n_time x n_channels⁠ numeric matrix, an hdemg_sim object from make_hdemg_sim(), or a PhysioExperiment (requires the PhysioCore package).
n_units	Maximum number of separation vectors to attempt (default 8).
extension_factor	Number of delayed copies per channel used to build the extended observations (default 16).
sampling_rate	Sampling rate in Hz; required when x is a bare matrix (taken from the object otherwise).
assay_name	Assay to use when x is a PhysioExperiment (default: the first assay).
refractory_sec	Minimum inter-discharge interval in seconds used for peak separation (default 0.02).
pnr_threshold	Minimum pulse-to-noise ratio in dB for an accepted unit (default 26).
sil_threshold	Minimum silhouette for an accepted unit (default 0.9).
max_iter	Maximum fixed-point iterations per source (default 100).
tol	Convergence tolerance on the separation vector (default 1e-4).
bandpass	Optional length-2 numeric c(low, high) in Hz to band-pass the signal (zero-phase) before decomposition. Recommended for real HD-sEMG (e.g. c(20, 500)) to remove baseline drift and high-frequency noise; NULL (default) leaves the signal unchanged.
seed	Optional integer random seed for reproducibility.

Details

The default quality gate (pnr_threshold = 26, sil_threshold = 0.9) is tuned for clean data (e.g. make_hdemg_sim()). Real HD-sEMG motor units typically have lower quality (PNR ~20-23 dB, silhouette ~0.85-0.92); for real recordings, band-pass the signal (bandpass = c(20, 500)) and relax the thresholds (pnr_threshold ~ 20, sil_threshold ~ 0.85). Validated against openhdemg on a real 64-channel recording: PhysioHDEMG recovered 4 of 5 reference motor units at rate-of-agreement 0.79-0.99 (lag-aligned).

Value

An hdemg_decomposition object with pulse_trains (list of discharge sample indices), a quality data frame (unit, n_discharges, discharge_rate, pnr_db, sil), the separation filters, and the retained signal/grid/sampling_rate for downstream firing statistics and action-potential maps.

References

Holobar, A. & Zazula, D. (2007). Multichannel blind source separation using convolution kernel compensation. IEEE TSP 55(9):4487-4496. Negro, F. et al. (2016). J Neural Eng 13(2):026027.

See Also

make_hdemg_sim(), muFiringStats(), muActionPotentialMaps()

Examples

sim <- make_hdemg_sim(n_units = 2, duration_sec = 3, grid = c(4, 4),
                      noise_sd = 0.005, seed = 1)
dec <- hdEMGDecompose(sim, n_units = 4, extension_factor = 10, seed = 1)
dec$quality

---

Wrap an HD-sEMG electrode grid in a PhysioExperiment

Description

Builds a PhysioExperiment (from PhysioCore) holding a single HD-sEMG assay (⁠n_time x n_channels⁠) at the given sampling rate, recording the electrode grid layout (rows, cols) in the channel metadata and object metadata so that hdEMGDecompose() can recover the grid for action-potential maps.

Usage

hdemgToPhysioExperiment(
  signal,
  sampling_rate,
  grid,
  unit = "uV",
  metadata = list()
)

Arguments

signal	An ⁠n_time x n_channels⁠ matrix of HD-sEMG (channels ordered column-major over the grid: row varies fastest).
sampling_rate	Sampling rate in Hz.
grid	Integer c(rows, cols) electrode-grid layout; rows * cols must equal the number of channels.
unit	Signal unit stored in the channel metadata (default "uV").
metadata	Optional named list stored in the object metadata.

Value

A PhysioExperiment object.

See Also

hdEMGDecompose(), make_hdemg_sim()

Examples

sim <- make_hdemg_sim(n_units = 2, duration_sec = 2, grid = c(4, 4), seed = 1)
# PhysioCore supplies the shared PhysioExperiment data model.
if (requireNamespace("PhysioCore", quietly = TRUE) &&
    requireNamespace("S4Vectors", quietly = TRUE)) {
  pe <- hdemgToPhysioExperiment(sim$signal, sim$sampling_rate, grid = c(4, 4))
  pe
}

---

Simulate a high-density surface EMG recording with known motor units

Description

Generates a synthetic HD-sEMG electrode-grid signal as the linear superposition of a known set of motor units, each firing at a target discharge rate with inter-spike-interval jitter and each with its own spatial action-potential map over the grid. The ground-truth discharge times are returned so that a decomposition can be validated (rate of agreement, discharge-rate error).

Usage

make_hdemg_sim(
  n_units = 4,
  duration_sec = 10,
  sampling_rate = 2048,
  grid = c(5, 5),
  firing_rates = NULL,
  muap_length = 24,
  cov_isi = 0.1,
  noise_sd = 0.02,
  seed = NULL
)

Arguments

n_units	Number of motor units to simulate (default 4).
duration_sec	Recording duration in seconds (default 10).
sampling_rate	Sampling rate in Hz (default 2048).
grid	Integer c(rows, cols) electrode-grid layout (default c(5, 5)); rows * cols channels are generated.
firing_rates	Optional numeric vector of per-unit mean discharge rates (Hz); defaults to an evenly spaced spread in ⁠[8, 18]⁠ Hz.
muap_length	MUAP template length in samples (default 24).
cov_isi	Coefficient of variation of the inter-spike interval (firing-time jitter, default 0.1).
noise_sd	Standard deviation of additive white noise, in the same units as the unit-peak MUAPs (default 0.02; smaller = higher SNR).
seed	Optional integer random seed for reproducibility.

Value

An hdemg_sim object (a list) with:

signal

An ⁠n_time x n_channels⁠ matrix of simulated HD-sEMG.

sampling_rate

Sampling rate in Hz.

grid

The c(rows, cols) grid layout.

spikes

A length-n_units list of integer ground-truth discharge sample indices.

discharge_rates

The realised per-unit mean discharge rate (Hz).

muap

A length-n_units list of ⁠muap_length x n_channels⁠ ground-truth action-potential templates.

n_units

Number of motor units.

References

Farina, D. et al. (2010). Decoding the neural drive to muscles from the surface electromyogram. J Physiol 588(11):1969-1975.

See Also

hdEMGDecompose()

Examples

sim <- make_hdemg_sim(n_units = 3, duration_sec = 2, grid = c(4, 4),
                      seed = 1)
dim(sim$signal)
lengths(sim$spikes)

---

Match decomposed motor units to a reference decomposition

Description

Scores a decomposition against a reference set of motor-unit discharge patterns with a lag-aware rate-of-agreement (RoA): for each reference unit the best-matching decomposed unit is found (allowing a constant latency of up to max_lag samples), and the RoA – the fraction of discharges the two share within tol samples, matched / (n_a + n_b - matched) – is reported. This scores cross-tool recovery ("did we recover the units a reference finds, and which do we miss") rather than bit-identical agreement, so a constant latency between two decompositions does not penalise a correct match.

Usage

matchMotorUnits(decomposition, reference, tol = 1L, max_lag = 0L)

Arguments

decomposition	An hdemg_decomposition from hdEMGDecompose().
reference	The reference decomposition: a list of integer discharge- sample vectors (one per reference unit), or a PhysioExperiment carrying them in metadata reference_mupulses (e.g. from readHDEMG()).
tol	Matching tolerance in samples (default 1).
max_lag	Maximum constant latency in samples to search when matching (default 0 = exact alignment; use a few milliseconds' worth for cross-tool scoring).

Value

A data frame (class mu_match), one row per reference unit, with columns ref_mu, n_discharges, matched_unit, roa_lag_aware, roa_no_lag, lag, pnr_db, sil.

References

Holobar A, Minetto MA, Farina D (2014). Accurate identification of motor unit discharge patterns from high-density surface EMG and validation with a novel signal-based performance metric. J Neural Eng 11:016008.

See Also

hdEMGDecompose(), readHDEMG(), muPulseTrains()

Examples

# Score a decomposition against the simulator's known discharge times.
sim <- make_hdemg_sim(n_units = 2, duration_sec = 3, grid = c(4, 4),
                      noise_sd = 0.005, seed = 1)
dec <- hdEMGDecompose(sim, n_units = 4, extension_factor = 10, seed = 1)
matchMotorUnits(dec, sim$spikes, tol = 2, max_lag = 5)

---

Motor-unit action-potential maps

Description

Spike-triggered average of the raw HD-sEMG for each motor unit, yielding the multi-channel action-potential waveform and its peak-to-peak amplitude across the electrode grid (the spatial "map" used to localise the unit).

Usage

muActionPotentialMaps(decomposition, window_ms = 15)

Arguments

decomposition	An hdemg_decomposition from hdEMGDecompose().
window_ms	Half-window in milliseconds averaged around each discharge (default 15).

Value

A length-n_units list; each element has waveform (⁠(2*w+1) x n_channels⁠ spike-triggered average), amplitude (per-channel peak-to-peak), map (the amplitude reshaped to the ⁠rows x cols⁠ grid, or NULL if the grid is unknown), and peak_channel (the strongest channel).

References

Farina, D. et al. (2016). Reflections on the evolution of surface EMG decomposition. J Electromyogr Kinesiol 26:39-46.

See Also

hdEMGDecompose()

Examples

sim <- make_hdemg_sim(n_units = 2, duration_sec = 3, grid = c(4, 4),
                      noise_sd = 0.005, seed = 1)
dec <- hdEMGDecompose(sim, n_units = 4, extension_factor = 10, seed = 1)
maps <- muActionPotentialMaps(dec)
maps[[1]]$peak_channel

---

Motor-unit firing statistics

Description

Per-unit discharge-rate and recruitment statistics: the mean discharge rate, the mean and coefficient of variation of the inter-spike interval, the recruitment (first discharge) and derecruitment (last discharge) times, and the decomposition quality (pulse-to-noise ratio and silhouette).

Usage

muFiringStats(decomposition)

Arguments

decomposition

An hdemg_decomposition from hdEMGDecompose().

Value

A data frame with one row per motor unit and columns unit, n_discharges, discharge_rate (Hz), mean_isi_sec, cov_isi, recruitment_sec, derecruitment_sec, pnr_db, sil.

References

De Luca, C.J. et al. (1982). Behaviour of human motor units in different muscles during linearly varying contractions. J Physiol 329:113-128.

See Also

hdEMGDecompose(), muPulseTrains()

Examples

sim <- make_hdemg_sim(n_units = 2, duration_sec = 3, grid = c(4, 4),
                      noise_sd = 0.005, seed = 1)
dec <- hdEMGDecompose(sim, n_units = 4, extension_factor = 10, seed = 1)
muFiringStats(dec)

---

Motor-unit discharge pulse trains

Description

Returns the discharge instants of every decomposed motor unit in tidy form.

Usage

muPulseTrains(decomposition)

Arguments

decomposition

An hdemg_decomposition from hdEMGDecompose().

Value

A data frame with one row per discharge: unit (integer), discharge (sample index), and time_sec (discharge time in seconds).

See Also

hdEMGDecompose(), muFiringStats()

Examples

sim <- make_hdemg_sim(n_units = 2, duration_sec = 3, grid = c(4, 4),
                      noise_sd = 0.005, seed = 1)
dec <- hdEMGDecompose(sim, n_units = 4, extension_factor = 10, seed = 1)
head(muPulseTrains(dec))

---

Read a high-density surface EMG recording from an HDF5 file

Description

Reads an HD-sEMG electrode-grid recording from an HDF5 (.h5) file into a PhysioExperiment, the ingest counterpart of hdemgToPhysioExperiment() for on-disk grids. The file must contain a signal dataset (the raw grid signal) and a sampling_rate scalar; it may also contain a grid dataset (c(rows, cols)) and a reference_mus group of per-unit integer discharge-sample vectors – a bundled reference decomposition (e.g. exported from the openhdemg sample file) – which are stored in the returned object's metadata as reference_mupulses for scoring a decomposition with matchMotorUnits().

Usage

readHDEMG(
  path,
  grid = NULL,
  signal_name = "signal",
  rate_name = "sampling_rate",
  unit = "uV"
)

Arguments

path	Path to the HDF5 file.
grid	Optional integer c(rows, cols) electrode-grid layout override; taken from the file's grid dataset when present, otherwise c(1, n_channels).
signal_name, rate_name	Dataset names for the signal and the sampling rate (defaults "signal" and "sampling_rate").
unit	Signal unit stored in the channel metadata (default "uV").

Details

The signal dataset may be stored ⁠n_time x n_channels⁠ or ⁠n_channels x n_time⁠; the longer axis is taken as time.

Value

A PhysioExperiment holding the HD-sEMG grid, with any bundled reference decomposition in metadata reference_mupulses.

See Also

hdEMGDecompose(), matchMotorUnits(), hdemgToPhysioExperiment()

Examples

# Round-trip a tiny grid through a temporary HDF5 file (needs rhdf5).
if (requireNamespace("rhdf5", quietly = TRUE) &&
    requireNamespace("PhysioCore", quietly = TRUE)) {
  sim <- make_hdemg_sim(n_units = 2, duration_sec = 2, grid = c(4, 4), seed = 1)
  f <- tempfile(fileext = ".h5")
  rhdf5::h5createFile(f)
  rhdf5::h5write(sim$signal, f, "signal")
  rhdf5::h5write(sim$sampling_rate, f, "sampling_rate")
  rhdf5::h5write(as.integer(sim$grid), f, "grid")
  rhdf5::H5close()
  pe <- readHDEMG(f)
  unlink(f)
  pe
}

---