Motor-unit decomposition of HD-sEMG

PhysioHDEMG decomposes high-density surface EMG (HD-sEMG) electrode-grid recordings into the discharge patterns of the individual motor units that produced them. This vignette runs the whole pipeline offline on a synthetic recording with known motor units, so every number can be checked against ground truth — no device or data file is needed.

library(PhysioHDEMG)
#> PhysioHDEMG v0.2.1 - HD-sEMG motor-unit decomposition

1. Simulate a recording with known motor units

make_hdemg_sim() superimposes a set of motor units — each with its own discharge train and its own spatial action-potential map over the electrode grid — and returns the ground-truth discharge times alongside the signal.

sim <- make_hdemg_sim(
  n_units      = 3,
  duration_sec = 4,
  grid         = c(4, 4),   # 16-channel electrode grid
  noise_sd     = 0.005,     # smaller = higher SNR
  seed         = 1
)
sim
#> <hdemg_sim> 3 units, 4 x 4 grid, 4.0 s @ 2048 Hz
#>   discharge rates: 7.2, 12.5, 17.5 Hz
dim(sim$signal)             # n_time x n_channels
#> [1] 8192   16
lengths(sim$spikes)         # ground-truth discharges per unit
#> [1] 29 50 70

2. Decompose

hdEMGDecompose() extends and whitens the channels, seeds each source by a fixed-point contrast maximisation, refines it by convolution-kernel compensation, and accepts units that clear a pulse-to-noise-ratio and silhouette quality gate.

dec <- hdEMGDecompose(sim, n_units = 6, extension_factor = 10, seed = 1)
dec
#> <hdemg_decomposition> 1 motor unit(s) from 16 channels @ 2048 Hz
#>   MU 1: 29 discharges, 7.2 Hz, PNR 30.0 dB, SIL 0.98
dec$quality
#>   unit n_discharges discharge_rate   pnr_db       sil
#> 1    1           29           7.25 29.97549 0.9815063

3. Inspect the recovered units

The accessors turn a decomposition into tidy tables. muPulseTrains() lists every discharge, muFiringStats() summarises discharge rate, recruitment, and quality per unit, and muActionPotentialMaps() recovers each unit’s spatial footprint on the grid.

head(muPulseTrains(dec))
#>   unit discharge  time_sec
#> 1    1       754 0.3676758
#> 2    1       989 0.4824219
#> 3    1      1258 0.6137695
#> 4    1      1533 0.7480469
#> 5    1      1803 0.8798828
#> 6    1      2052 1.0014648
muFiringStats(dec)
#>   unit n_discharges discharge_rate mean_isi_sec    cov_isi recruitment_sec
#> 1    1           29           7.25     0.126552 0.09035235       0.3676758
#>   derecruitment_sec   pnr_db       sil
#> 1          3.911133 29.97549 0.9815063

maps <- muActionPotentialMaps(dec)
maps[[1]]$peak_channel     # strongest electrode for the first unit
#> [1] 1

4. Validate against the ground truth

Because the simulator recorded the true discharge times, we can score the decomposition with a lag-aware rate of agreement. roa_lag_aware near 1 means a reference unit was recovered; lag is the constant latency between the two.

matchMotorUnits(dec, sim$spikes, tol = 2, max_lag = 5)
#>   ref_mu n_discharges matched_unit roa_lag_aware roa_no_lag lag   pnr_db
#> 1      1           29            1    0.00000000 0.00000000  -5 29.97549
#> 2      2           50            1    0.03947368 0.02597403  -2 29.97549
#> 3      3           70            1    0.04210526 0.04210526   0 29.97549
#>         sil
#> 1 0.9815063
#> 2 0.9815063
#> 3 0.9815063

Working with your own recordings

On real data, read a grid from an HDF5 file with readHDEMG() (which can also carry a bundled reference decomposition for matchMotorUnits()), or wrap an in-memory matrix with hdemgToPhysioExperiment() to join the shared PhysioExperiment data model. Real motor units have lower quality than this clean simulation, so band-pass the signal (bandpass = c(20, 500)) and relax the gates (pnr_threshold, sil_threshold) — see ?hdEMGDecompose.

Conventional (non-decomposed) EMG features live in the sibling PhysioEMG package.