HarmonICA: Quasi-Time-Invariant ICA for Motor Unit Decomposition during Dynamic Contractions

Alexander Kenneth Clarke; Dimitrios Halatsis; Agnese Grison; Irene Mendez Guerra; Noura Ezaz-Nikpay; Pranav Mamidanna; Shihan Ma; Silvia Muceli

cs.HCJun 28, 2024

HarmonICA: Quasi-Time-Invariant ICA for Motor Unit Decomposition during Dynamic Contractions

Authors: Alexander Kenneth Clarke, Dimitrios Halatsis, Agnese Grison, Irene Mendez Guerra, Noura Ezaz-Nikpay, Pranav Mamidanna, Shihan Ma, Silvia Muceli,

Organizations: Imperial College London, London, UK · Chalmers University of Technology, Gothenburg, Sweden

Abstract

Decomposing surface electromyography (EMG) into the spike trains of individual motor neurons is a long-standing inverse problem and a key step toward motor-neuron-driven neural interfaces such as prosthetics and exoskeletons. The standard approach, independent component analysis (ICA) of the multichannel signal, assumes that the mixing from neurons to electrodes is stationary in time. This assumption fails during movement, when volume-conductor deformation makes the mixing time-varying, and current decomposition algorithms are correspondingly restricted to isometric contractions. We introduce HarmonICA, a quasi-time-invariant ICA formulation in which a static linear separator is preceded by a learned, low-rank, time-varying invertible transformation. The separator is trained with an independence loss on the uncompensated projection, and the transformation with a stationarity loss on the recovered source. Gradients are not shared between the two, so source extraction remains a classical linear-ICA problem while non-stationary distortion is absorbed by the transformation. The closed-form inverse of the transformation enables per-spike subtraction with a time-varying template during sequential peel-off. On a public simulated benchmark of dynamic high-density EMG with ground-truth spike trains, the method outperforms five adaptive ICA baselines at every recall threshold, recovering more units at a higher accuracy.

Figures & tables

Figure 1: (a) A static ICA filter $\mathbf{w}_{j}$ recovers spikes only in the narrow local window where the mixing matches it; outside this window, the source approaches the noise floor and roughly half of the ground-truth spikes are missed. (b) HarmonICA’s time-varying compensator $A_{j}(t)$ rotates the filter (equivalently, the signal space) so that the same source extracts cleanly across the full recording. (c) The learned rotation tracks the angle-driven MUAP variation. (d) The two-component GMM stationarity loss propagates this rotation reliably; classical FastICA non-Gaussianity contrasts are less reliable on this benchmark (see supplementary § H.1 ).

Panel of Figure 2 — Figure 2: HarmonICA single-source extraction. (a) Forward pass: time encoding $\bm{\gamma}(t)$ , per-source network $\mathrm{NN}_{j}$ , low-rank invertible compensator $A_{j}(t)=I+U_{j}(t)V_{j}^{\top}$ applied to the whitened signal, and static separator $\mathbf{w}_{j}$ . (b) Decoupled training: independence step on $\mathbf{w}_{j}$ ( $L_{I}$ ) with the compensator off, and after a $T_{w}$ -step warmup a compensation step on $\phi_{j}$ ( $L_{N}$ ) with $\mathbf{w}_{j}$ detached and the spike/null sets re-split every $K_{g}$ steps. (c) Per-spike peel-off: subtract the spike-triggered average mapped back through $A_{j}(t_{j,i}+\tau)^{-1}$ , then restart on the residual.

Figure 3: Yield

\times

recall Pareto curve (ch320 benchmark, 168 recordings

\times

5 seeds).

x

-axis: average MUs per recording-seed at recall

\geq R

;

y

-axis: recall threshold

R

. Higher and further-right is better. Dashed line: unrefined ICA / cBSS reference; filled circles: each method’s coverage cap. Mendez Guerra 2024 is run with hyperparameters tuned on the validation split (§ 4.3 ). Per-grid Pareto-crossover yields (including the 70-channel grid) are in supplementary Table 13 .

Method	ch070 (70-channel)			ch320 (320-channel)
Method	$n_{\text{gt}}<40$	$n_{\text{gt}}\,40\!-\!69$	$n_{\text{gt}}\geq 70$	$n_{\text{gt}}<40$	$n_{\text{gt}}\,40\!-\!69$	$n_{\text{gt}}\geq 70$
ICA / cBSS [ 6 ]	$36\pm 2.5$	$151\pm 3.0$	$62\pm 1.9$	$40\pm 2.9$	$142\pm 2.9$	$67\pm 4.2$
HarmonICA (PE, ours)	$90\pm 4.6$	$252\pm 24$	$134\pm 5.5$	$96\pm 2.2$	$263\pm 4.4$	$156\pm 5.3$
HarmonICA ( $\theta$ , ours)	$\mathbf{91\pm 4.7}$	$\mathbf{265\pm 11}$	$\mathbf{146\pm 13}$	$\mathbf{102\pm 4.6}$	$\mathbf{290\pm 5.3}$	$\mathbf{179\pm 5.5}$
Chen 2020 [ 14 ]	$36\pm 2.9$	$133\pm 4.3$	$59\pm 3.4$	$39\pm 4.0$	$149\pm 8.7$	$81\pm 5.4$
Kramberger 2021 [ 15 ]	$20$ —	$59$ —	$25$ —	$25$ —	$73$ —	$35$ —

Table 1: Total motor units recovered at F

{}_{1}\geq 0.95

, stratified by per-recording motor-unit count

n_{\text{gt}}

. Cells are bin-level yield sums, mean

\pm

s.d. across

5

FastICA-initialisation seeds (

10

for HarmonICA (

\theta

); Kramberger is deterministic and is reported on a single seed; “ —” in the s.d. position). Mendez Guerra 2024 is run with hyperparameters tuned on the validation split (§ 4.3 ). Higher is better; bold marks the best method per cell. Per-recording mean yields with across-recording dispersion are in supplementary Table 7 ; all HarmonICA

-

baseline differences are significant at

p<10^{-19}

(paired Wilcoxon signed-rank, supplementary § D.1 ).

Appendix figures & tables36 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 4: The EMG decomposition problem. A motor-unit pool sends spike trains down its motor neurons to muscle fibres; the fibre depolarisations are picked up at the skin surface as a multichannel signal that mixes the contributions of every active motor unit through MUAP-shaped templates. The decomposition target is to recover the underlying per-neuron spike trains from the multichannel surface signal. Adapted from [ 24 ] .

Figure 5: Stationary vs. non-stationary mixing. Under isometric conditions a single static ICA filter $\mathbf{w}_{j}$ extracts the full spike train of a motor unit cleanly. As joint angle and limb position vary, each motor unit’s MUAP rotates and deforms; a fixed filter then matches the projection only inside a local window and the recovered source decays outside it.

Figure 6: MUAP variation across joint angles. The MUAP of a single motor unit traced across the angular sweep of the simulated joint. The waveform shape rotates and rescales smoothly with angle; HarmonICA’s compensator $A_{j}(t)$ targets this low-rank, slowly varying part of the mixing.

$T$	# pairs	pair-less MUs	paired MUs	mean cluster size	max cluster size
0.75	110 377	61 (0.7 %)	8 456 (99.3 %)	158.5	224
0.85	47 893	312 (3.7 %)	8 205 (96.3 %)	105.0	224
0.90	23 132	988 (11.6 %)	7 529 (88.4 %)	36.1	210
0.95	6 171	3 428 (40.3 %)	5 089 (59.8 %)	5.5	100
0.99	156	8 228 (96.6 %)	289 (3.4 %)	2.1	4

Appendix

Table 2: Pair- and cluster-level structure of within-pool MUAP similarity, aggregated across all 45 MUniverse motor-unit pools (

N_{\text{full}}=8{,}517

motor units). At any decomposability-relevant threshold the typical pool collapses into a single connected component that absorbs nearly every motor unit; the pair-less / paired transition only occurs above

T=0.95

	% pair-less motor units			spread
$T$	early (step 10)	mid (step 65)	late (step 120)	(pp)
0.75	0.58	0.72	0.80	0.22
0.85	3.35	3.66	3.82	0.47
0.90	9.51	11.60	10.83	2.09
0.95	36.48	40.25	39.10	3.77
0.99	95.26	96.61	95.02	1.59

Appendix

Table 3: Cluster topology under early / mid / late joint poses. Pose-spread is at most 3.8 percentage points at any threshold, and the maximum cluster size at

T=0.75

is invariant (224 across all three poses), so the same motor units sit in the same dominant component regardless of joint angle.

threshold $T$	kept MUs / pool (mean $\pm$ std)	retention
0.75	$25.5\pm 7.5$	13.5 %
0.80	$35.0\pm 10.6$	18.5 %
0.85	$50.5\pm 14.9$	26.7 %
0.90	$77.2\pm 20.9$	40.8 %

Appendix

Table 4: Decomposable-subset retention under a greedy Henneman-ordered crosstalk filter. Even at a permissive

T=0.90

, which admits MUAP pairs

90\,\%

identical, only

40.8\,\%

of nominal motor units survive; at the standard

T=0.75

only

13.5\,\%

do. Reporting raw

F_{1}

against the nominal pool size therefore under-estimates every method’s ceiling by an

\sim 3

–

7\times

factor that has nothing to do with method quality.

Tier	Pool size (MUs)	Threshold $T$	# pools	# recordings
small	26–45	$0.85$	6	96
medium	47–74	$0.85$	8	128
large	71–113	$0.85$ / $0.90$	7	112
total	26–113	—	21	336

Appendix

Table 5: Structure of the muniverse-benchmark-dynamic test split. Pools are stratified by decomposable size after the tier-dependent crosstalk filter; each pool is rendered under 8 angle/SNR conditions

\times

2 electrode grids (

16

recordings/pool). An additional held-out val split (5 pools, 80 recordings) is included in the public release for groups that wish to use it; all main-paper numbers and the supplementary tables that mirror §5 are on the

168

-per-grid test split.

Method	$F_{1}\geq 0.99$		$F_{1}\geq 0.95$		$F_{1}\geq 0.90$		$F_{1}\geq 0.80$
Method	ch070	ch320	ch070	ch320	ch070	ch320	ch070	ch320
ICA / cBSS [ 6 ]	0.97	0.93	1.48	1.49	1.82	1.84	2.64	2.81
HarmonICA (ours)	1.94	2.15	2.83	3.06	3.44	3.67	4.59	5.01
Chen 2020 [ 14 ]	0.50	0.50	1.36	1.60	1.87	2.30	2.94	3.41
Kramberger 2021 [ 15 ]	0.39	0.44	0.72	0.91	0.92	1.08	1.19	1.36
Yeung 2024 [ 16 ]	0.24	0.28	0.72	0.65	1.03	1.05	1.67	1.70

Appendix

Table 6: Yield at four

F_{1}

quality floors . Mean MUs per recording at the indicated

F_{1}

floor, mirroring the main-paper convention (per-recording mean over the

5

ICA-initialisation seeds, then mean across recordings within a density). Higher is better; bold marks the best method in each column. Kramberger 2021 is averaged over the recordings on which it returns

\geq 1

MU only (§ D ).

Method	ch070 (70-channel)			ch320 (320-channel)
Method	$n_{\text{gt}}<40$	$n_{\text{gt}}\,40\!-\!69$	$n_{\text{gt}}\geq 70$	$n_{\text{gt}}<40$	$n_{\text{gt}}\,40\!-\!69$	$n_{\text{gt}}\geq 70$
ICA / cBSS [ 6 ]	0.90 $\pm$ 1.31	1.88 $\pm$ 2.39	1.29 $\pm$ 2.38	1.01 $\pm$ 1.42	1.78 $\pm$ 2.30	1.40 $\pm$ 2.19
HarmonICA (ours)	2.25 $\pm$ 1.69	3.15 $\pm$ 2.45	2.78 $\pm$ 2.31	2.39 $\pm$ 1.75	3.28 $\pm$ 2.45	3.25 $\pm$ 2.64
Chen 2020 [ 14 ]	0.91 $\pm$ 1.09	1.66 $\pm$ 1.66	1.23 $\pm$ 1.57	0.98 $\pm$ 1.31	1.87 $\pm$ 1.94	1.68 $\pm$ 1.91
Kramberger 2021 [ 15 ]	0.51 $\pm$ 0.88	0.83 $\pm$ 0.89	0.74 $\pm$ 0.71	0.64 $\pm$ 0.87	0.99 $\pm$ 1.03	1.06 $\pm$ 0.93
Yeung 2024 [ 16 ]	0.56 $\pm$ 0.46	0.83 $\pm$ 0.84	0.67 $\pm$ 0.66	0.50 $\pm$ 0.44	0.78 $\pm$ 0.89	0.56 $\pm$ 0.50

Appendix

Table 7: Per-recording mean yield at F

{}_{1}\geq 0.95

stratified by per-recording motor-unit count

n_{\text{gt}}

. Each cell is mean MUs/rec

\pm

standard deviation across recordings within the cell, after each recording is collapsed to its mean over the

5

FastICA-initialisation seeds. Bin recording counts are

40

80

48

for the three tiers. Standard deviations are large because the test split is heterogeneous (per-recording pool sizes range

23

–

101

MUs); the corresponding across-seed standard deviation per cell is given in Table 8 . Kramberger 2021 is averaged over the recordings on which it returns

\geq 1

MU only (§ D ).

Method	ch070 (70-channel)			ch320 (320-channel)
Method	$n_{\text{gt}}<40$	$n_{\text{gt}}\,40\!-\!69$	$n_{\text{gt}}\geq 70$	$n_{\text{gt}}<40$	$n_{\text{gt}}\,40\!-\!69$	$n_{\text{gt}}\geq 70$
ICA / cBSS [ 6 ]	0.90 $\pm$ 0.06	1.88 $\pm$ 0.04	1.29 $\pm$ 0.04	1.01 $\pm$ 0.07	1.78 $\pm$ 0.04	1.40 $\pm$ 0.09
HarmonICA (ours)	2.25 $\pm$ 0.11	3.15 $\pm$ 0.30	2.78 $\pm$ 0.12	2.39 $\pm$ 0.05	3.28 $\pm$ 0.06	3.25 $\pm$ 0.11
Chen 2020 [ 14 ]	0.91 $\pm$ 0.07	1.66 $\pm$ 0.05	1.23 $\pm$ 0.07	0.98 $\pm$ 0.10	1.87 $\pm$ 0.11	1.68 $\pm$ 0.11
Kramberger 2021 [ 15 ]	0.51 $\pm$ —	0.83 $\pm$ —	0.74 $\pm$ —	0.64 $\pm$ —	0.99 $\pm$ —	1.06 $\pm$ —
Yeung 2024 [ 16 ]	0.56 $\pm$ 0.04	0.83 $\pm$ 0.04	0.67 $\pm$ 0.04	0.50 $\pm$ 0.06	0.78 $\pm$ 0.05	0.56 $\pm$ 0.04

Appendix

Table 8: Seed reproducibility of the per-recording mean yield at

F_{1}\geq 0.95

. Same shape and units as Table 7 (mean MUs/rec) but with the

\pm

taken across the

5

FastICA-initialisation seeds rather than across recordings. Kramberger 2021 is averaged over the recordings on which it returns

\geq 1

MU only (§ D ).

Figure 7: Pairwise MU-set overlap at $F_{1}\geq 0.95$ . Left: Jaccard similarity (symmetric). Right: % of row method’s accepted MUs also accepted by the column method (asymmetric; diagonal $=100$ ). adapt : Mendez Guerra 2024 with hyperparameters tuned by the method’s author on the validation split.

Figure 8: Per-method MU-discovery breakdown vs. the unrefined FastICA seed at

F_{1}\geq 0.95

. Each cluster of bars is one method; within each cluster the three bars correspond to GT motor-unit count bins (

<40

40

–

69

\geq 70

). Stack components: both = GT MU found by both method and ICA, method only = GT MU found by method but not by ICA (novel), ICA only = GT MU found by ICA but not by method (missed). adapt : Mendez Guerra 2024 with hyperparameters tuned by the method’s author on the validation split.

	ch070			ch320
Method	mean	95% CI lo	95% CI hi	mean	95% CI lo	95% CI hi
HarmonICA (ours)	$\mathbf{2.83}$	$2.51$	$3.14$	$\mathbf{3.06}$	$2.71$	$3.41$
ICA / cBSS reference	$1.48$	$1.16$	$1.82$	$1.49$	$1.18$	$1.82$
Chen 2020	$1.36$	$1.13$	$1.60$	$1.60$	$1.34$	$1.88$
Kramberger 2021	$0.72$	$0.59$	$0.86$	$0.91$	$0.75$	$1.06$
Yeung 2024	$0.72$	$0.61$	$0.83$	$0.65$	$0.54$	$0.76$

Appendix

Table 9: Bootstrap

95\%

confidence interval on per-recording mean yield at

F_{1}\geq 0.95

. The HarmonICA interval does not overlap with any baseline interval on either grid. Kramberger 2021 is averaged over the recordings on which it returns

\geq 1

MU only (§ D ).

HarmonICA vs.	grid	$n_{\text{rec}}$	$\Delta$ (MUs/rec)	$p$ -value	wins	losses (ties)
ICA / cBSS reference	ch070	$168$	$+1.35$	$2.4\!\times\!10^{-20}$	$150$	$\phantom{0}10$ ( $\phantom{0}8$ )
ICA / cBSS reference	ch320	$168$	$+1.57$	$2.3\!\times\!10^{-22}$	$142$	$\phantom{0}14$ ( $12$ )
Chen 2020	ch070	$168$	$+1.47$	$2.0\!\times\!10^{-26}$	$148$	$\phantom{00}6$ ( $14$ )
Chen 2020	ch320	$168$	$+1.46$	$6.0\!\times\!10^{-26}$	$148$	$\phantom{00}9$ ( $11$ )
Kramberger 2021	ch070	$144$	$+2.23$	$2.8\!\times\!10^{-24}$	$136$	$\phantom{00}2$ ( $\phantom{0}6$ )
Kramberger 2021	ch320	$146$	$+2.35$	$2.2\!\times\!10^{-25}$	$143$	$\phantom{00}1$ ( $\phantom{0}2$ )

Appendix

Table 10: Paired Wilcoxon signed-rank tests on per-recording yield at

F_{1}\geq 0.95

. “Wins” is the count of recordings where HarmonICA exceeds the baseline; “losses” the reverse; “ties” where they are equal.

\Delta

is the mean per-recording difference in MUs at

F_{1}\geq 0.95

. All p-values are far below the Bonferroni-corrected threshold for the

12

tests reported here.

HarmonICA pipeline
Hardware	1 $\times$ NVIDIA A40 GPU
Framework	PyTorch 2.x, CUDA 12.x
Sampling rate $f_{s}$	$2048$ Hz
Toeplitz extension factor $L$	$8$
PCA dimensionality $K$	$128$
Whitening	ZCA

Appendix

Table 11: Implementation details for HarmonICA and the baselines.

Grid	Side info	$n_{\text{gt}}<40$	$n_{\text{gt}}\,40\!-\!69$	$n_{\text{gt}}\geq 70$
ch070	PE	$90\pm 4.6$	$252\pm 24$	$134\pm 5.5$
ch070	angle	$\mathbf{91\pm 4.7}$	$\mathbf{265\pm 11}$	$\mathbf{146\pm 13}$
ch320	PE	$96\pm 2.2$	$263\pm 4.4$	$156\pm 5.3$
ch320	angle	$\mathbf{102\pm 4.6}$	$\mathbf{290\pm 5.3}$	$\mathbf{179\pm 5.5}$

Appendix

Table 12: HarmonICA bin-level total motor units recovered at

F_{1}\geq 0.95

under PE (sinusoidal positional encoding, blind) versus angle (real joint-angle trace) compensator inputs. Mean

\pm

s.d. over the

5

FastICA initialisation seeds (

10

for angle). Same test split as main Table 1 . Bin recording counts are

40

80

48

Grid	Method	$R\geq 0.99$	$R\geq 0.95$	$R\geq 0.90$	$R\geq 0.80$
ch070	HarmonICA (PE)	$2.23$	$3.16$	$3.80$	$4.71$
	HarmonICA (angle)	$2.36$	$3.34$	$3.96$	$4.85$
	ICA / cBSS	$0.94$	$1.42$	$1.67$	$2.11$
	Mendez Guerra 2024	$0.57$	$1.39$	$1.82$	$2.61$
ch320	HarmonICA (PE)	$2.17$	$3.24$	$3.84$	$4.70$
	HarmonICA (angle)	$2.49$	$3.66$	$4.21$	$5.02$

Appendix

Table 13: Pareto-crossover yield (MUs per recording-seed at recall

\geq R

) for HarmonICA under PE vs. angle, alongside the unrefined ICA / cBSS reference (which is identical under both regimes since it does not see the side input) and Mendez Guerra 2024, the strongest adaptive baseline at

R\geq 0.95

Loss	Mean $F_{1}$	$n$ at $F_{1}\geq 0.85$	$n\geq 0.90$	$n\geq 0.95$
ch070
GMM (ours)	$\mathbf{0.706\pm 0.004}$	$\mathbf{565.8\pm 8.1}$	$\mathbf{491.6\pm 5.1}$	$\mathbf{404.6\pm 2.5}$
Negentropy ( $g_{1}$ )	$0.616\pm 0.008$	$354.0\pm 14.0$	$282.8\pm 5.0$	$226.8\pm 6.8$
Negentropy ( $g_{2}$ )	$0.608\pm 0.009$	$341.8\pm 8.1$	$275.0\pm 6.0$	$221.8\pm 7.0$
Power- $p$ ( $\|\hat{s}\|^{4}$ )	$0.399\pm 0.004$	$132.6\pm 6.1$	$117.0\pm 5.6$	$102.2\pm 6.7$
Kurtosis	$0.296\pm 0.005$	$159.8\pm 6.4$	$136.6\pm 8.4$	$104.6\pm 5.4$

Appendix

Table 14: Stationarity-loss substitution on the per-source iso-ablation. Yield is per seed, summed over

168

recordings. Reference: GMM, our headline choice.

Figure 9: Rescue. Only the GMM loss and the power- $p$ loss work cleanly; kurtosis maximises a few spikes, the negentropy losses fail to discover the harder spikes, and power- $p$ generally succeeds in the easier cases.

Figure 10: Partial success. GMM yields the best result. The negentropy losses still fail to discover a substantial fraction of spikes; kurtosis still amplifies the central part of the source; power- $p$ regresses.

Figure 11: Regression. GMM and the negentropy losses regress by losing precision, while kurtosis and power- $p$ lose recall. GMM regresses slightly more than the negentropy losses (which are more conservative in general); power- $p$ still amplifies the bigger spikes and loses accuracy as a result; kurtosis amplifies only a few spikes.

Rank	Mean $F_{1}$	$n$ at $F_{1}\geq 0.85$	$n\geq 0.90$	$n\geq 0.95$
ch070
$r=3$	$0.736\pm 0.003$	$694.4\pm 7.3$	$602.0\pm 11.3$	$501.8\pm 8.6$
$r=5$	$0.740\pm 0.002$	$710.6\pm 10.5$	$618.8\pm 10.2$	$510.0\pm 9.4$
$r=7$	$\mathbf{0.741\pm 0.002}$	$\mathbf{716.6\pm 11.0}$	$\mathbf{624.6\pm 11.1}$	$\mathbf{512.6\pm 9.9}$
ch320
$r=3$	$0.759\pm 0.003$	$775.0\pm 13.1$	$677.8\pm 10.0$	$566.0\pm 6.7$

Appendix

Table 15: Compensator rank sweep on the angle-conditioned pipeline.

r=3

is the operating point used in §5 of the main paper for both the PE and angle variants. Yield is per seed, summed over

168

recordings; all three ranks come from the same sweep, run separately from Table 12 , so the

r=3

row differs from it by run-to-run variation.

Method	ch070	ch320
Linear ICA (reference)	$10.0$	$8.4$
TCL (angle-conditioned)	$1$ – $4$	$2.4$
TCL (time segments)	$0.6$	–
iVAE (angle-conditioned)	$0.6$	$0.4$

Appendix

Table 16: Nonlinear ICA baselines on muniverse-benchmark-dynamic . Mean motor units per recording at

F_{1}>0.5

. The ch070 TCL (angle-conditioned) result varied across runs; we report the observed range. Evaluated on

5

recordings per electrode grid with a single seed, not the full test split; iVAE’s latent dimension was set to the ground-truth number of units.

	ch070		ch320
Method	alone	best-of	alone	best-of
Linear ICA	$2.16$	–	$0.75$	–
Chen 2020	$2.96$	$3.41$	$3.43$	$3.52$
Yeung 2024	$1.67$	$2.79$	$1.77$	$2.06$
Kramberger 2021	$1.02$	$2.59$	$1.21$	$1.84$
Glaser 2018	$0.88$	$2.25$	$0.83$	$1.17$

Appendix

Table 17: Oracle upper bound. Motor units per recording at

F_{1}\geq 0.8

; “best-of” scores each unit at

\max(F_{1}^{\mathrm{ICA}},F_{1}^{\mathrm{method}})

. This analysis uses its own recording panel and extension factor, so absolute values are not comparable with Table 6 .

Method	baseline only	HarmonICA only	ratio
Chen 2020	$0.49$	$1.53$	$3.1\times$
Yeung 2024	$0.27$	$3.06$	$11\times$
Glaser 2018	$0.07$	$3.89$	$56\times$
Kramberger 2021	$0.07$	$3.93$	$56\times$

Appendix

Table 18: Per-unit coverage on ch070 (same panel as Table 17 ). Motor units per recording at

F_{1}\geq 0.7

recovered by one method and missed by the other.

Drift	ICA	HarmonICA	Chen	Yeung	Kramberger	Glaser
$0.0$ (static)	$29.7$	$29.0$	$17.7$	$15.3$	$25.0$	$18.7$
$0.3$	$22.3$	$24.7$	$16.7$	$14.0$	$14.3$	$15.7$
$0.6$	$8.3$	$13.3$	$9.3$	$8.3$	$8.7$	$7.0$
$0.9$	$3.7$	$8.0$	$4.0$	$2.3$	$4.0$	$2.3$

Appendix

Table 19: Drift-controlled generic simulation. Motor units per recording at

F_{1}\geq 0.7

; drift

0

is a static mixing.

Method	MU3	MU4	MU1+MU14	MU11	MU15	Mean	Worst run
Chen 2020 [ 14 ]	$0.986$	$0.979$	$0.993$	$0.932$	$0.700$	$0.918$	$0.102$
HarmonICA (ours)	$0.999$	$0.969$	$0.989$	$0.874$	$0.757$	$0.918$	$0.732$
Kramberger 2021 [ 15 ]	$0.979$	$0.979$	$0.989$	$0.809$	$0.773$	$0.906$	$0.773$
Yeung 2024 [ 16 ]	$0.980$	$0.980$	$0.989$	$0.910$	$0.574$	$0.887$	$0.079$
Linear ICA (reference)	$0.986$	$0.832$	$0.886$	$0.837$	$0.727$	$0.854$	$0.706$
Glaser 2018 [ 13 ]	$0.956$	$0.988$	$0.861$	$0.597$	$0.622$	$0.805$	$0.438$

Appendix

Table 20: Real dynamic HD-sEMG:

F_{1}

per validated unit against the self-consistency reference, mean over

5

seeds, for

6

methods. “Worst run” is the lowest single-run

F_{1}

over all

25

unit–seed runs.

Figure 12: Where linear ICA fails on real dynamic HD-sEMG. The two units whose accuracy moves materially under compensation; the other three validated units shift by $\leq 0.04$ . Left: the motor-unit action potential on its best channel (black) and the eight next-best (grey). Right: the first $20$ s of the recording through the static linear-ICA filter (lower) and through HarmonICA’s compensated filter (upper), each normalised by its own baseline, with ticks marking detected discharges. On MU1 $+$ MU14 linear ICA detects nothing for the first $8$ s and $60$ discharges over the $20$ s shown, against HarmonICA’s $163$ . Across the whole recording, $96\%$ of the $109$ discharges HarmonICA recovers on MU1 $+$ MU14, and $80\%$ of the $46$ on MU4, fall in the first third. HarmonICA seed $1909$ .

Method	$0^{\circ}$	$-10^{\circ}$	$-20^{\circ}$	$-30^{\circ}$	$-40^{\circ}$
Chen 2020	$1.000$	$0.980$	$0.979$	$1.000$	$0.972$
Yeung 2024	$1.000$	$0.980$	$1.000$	$1.000$	$0.982$
Kramberger 2021	$1.000$	$1.000$	$1.000$	$1.000$	$0.943$

Appendix

Table 21: Tracking accuracy of the most-drifting unit with calibration cut to a single

5.3

s hold at

0^{\circ}

Figure 13: Easy success. ICA cleanly extracts the spikes with a single filter; the non-stationarity is trivial. These are the straightforward ICA successes within the benchmark: the angle modulation is small and the filter changes too little to matter.

Figure 14: ICA complete failure. A scenario where ICA would fail even in the absence of non-stationarity; HarmonICA cannot refine pure noise, so both methods fail completely. This accounts for most of the unrecovered motor units in the benchmark, with the failure driven by some combination of amplitude, crosstalk, and non-stationarity.

Figure 15: Stuck source rescued. ICA yields spikes around a slow-moving non-stationarity; the spikes are clean but weaken as $t$ moves away from the correct filter window. HarmonICA successfully rotates the filter to discover the missing spikes, lifting recall from $\sim\!0.5$ to $1$ .

Figure 16: Partial success. HarmonICA discovers part of the spike train but cannot extend coverage to the entire signal. This happens when the filter rotates too much or the local degradation of the source is too abrupt for the compensator to track gracefully.

Figure 17: Noise amplification. ICA yields a low-precision source; HarmonICA has no way of distinguishing correct spikes from incorrect ones in this regime, since no auxiliary information breaks the tie. The compensator ends up amplifying false-positive spikes and the source regresses.

Figure 18: Zoomed-in view of the noise-amplification failure mode of the previous figure.

Symbol	Meaning	Value / shape
Signal and preprocessing
$f_{s}$	sampling rate	$2048$ Hz
$N$	number of recording channels	$70$ or $320$
$L$	Toeplitz extension factor (delayed copies / channel)	$8$
$Q$	MUAP length	samples
$K$	PCA-whitened dimensionality	$128$

Appendix

Table 22: Symbols used in §3 of the main paper and in supplementary § E .

Independent Component Analysis

Independent component analysis separates observed mixtures into statistically independent latent signals for source recovery and data analysis.

See all 23 papers

Oct 6, 2026

Neuromotor Hierarchy Network: Physiological Inductive Biases for Robust Generalization in sEMG Decoding

Surface electromyography (sEMG) provides a wearable, noninvasive interface to neuromuscular activity for movement decoding and human-computer interaction. Population-scale decoding remains difficult because the relationship between sEMG and neuromuscular activity varies across users and sessions, while task-relevant dynamics span channels and multiple timescales. Learning waveform-to-output mappings from task labels leaves the distinction between recording variability and coordinated motor activity implicit. We introduce the Neuromotor Hierarchy Network (NHN), which learns a compact latent neuromotor state from task supervision to represent task-relevant neuromuscular coordination. NHN constructs this latent state through a hierarchy inspired by neuromotor organization. It adapts recording statistics while preserving relative intensity. Its spatiotemporal encoder uses parameter-efficient channel interactions and modulates features with multi-timescale history. The resulting features yield candidate activations of learned motor primitives, which are temporally integrated and continuously weighted to form the state. Theoretical analysis characterizes the efficiency, temporal behavior, and optimization of NHN's core mechanisms. We evaluate the architecture for both continuous hand-pose estimation on emg2pose and touch-typing recognition on emg2qwerty. On emg2pose, NHN reduces user-averaged angular error by 0.52% to 2.84% across all three generalization splits in both Regression and Tracking relative to Hadidi et al.'s best task-specific variants, using 48.42% to 48.51% fewer parameters. On emg2qwerty, NHN reduces beam-search character error rate by 19.40% zero-shot and 30.42% after fine-tuning relative to SplashNet-Upscale, using 65.86% fewer parameters. Physiology-guided inference of a latent neuromotor state supports parameter-efficient sEMG decoding.

He Wang, Hongyuan Qi, Zhaoxian Zhang +4

National University of Singapore

Neural Decoding State Representation Learning Surface Electromyography

Ask PDF

May 8, 2026

Estimation of Motor Unit Parameters from Surface Electromyograms using an Informed Autoencoder

Motor unit parameters such as the innervation zone centre or the conduction velocity of the electrical potential harbour the potential to improve the fidelity of neuromechanical models used for movement and force prediction. Determining these parameters in a non-invasive way is challenging, as they are subject-specific and may vary with muscle contraction. Existing work on the estimation of motor unit parameters mainly relies on white-box modelling and therefore requires substantial manual modelling effort. This work targets the simultaneous estimation of multiple subject-specific motor unit parameters from electromyography (EMG) recordings measured non-invasively at the skin surface. This results in an inverse problem with a nonlinear loss function. To address this problem, an informed autoencoder is developed. This autoencoder reconstructs the surface EMG recordings while learning the parameters in its latent space and adhering to physical laws that relate the parameters to the EMG signals. In experiments on synthetic data, innervation zone centres are estimated with a mean absolute error of 2.5989

\mathrm{mm}

, and conduction velocities of the electric potential are estimated with a mean absolute error of 0.1697

\mathrm{m}\mathrm{s}^{-1}

. These results demonstrate the plausibility of this novel approach, which enables the simultaneous estimation of several motor unit parameters while reducing manual modelling effort through the integration of data-driven machine learning.

Kaja Balzereit, Malte Mechtenberg, Axel Schneider

Hochschule Bielefeld, University of Applied Sciences and Arts, Institute for System Dynamics

Parameter Estimation Surface Electromyography

Ask PDF

Jun 12, 2026

Test-Time Adaptation of Spiking Neural Networks for Intracortical Neural Decoding using Membrane Potential Alignment

Intracortical brain-computer interfaces suffer from day-to-day neural signal shifts that degrade pretrained decoders. Existing unsupervised adaptation methods rely on deep recurrent or adversarial architectures that are too computationally expensive for implantable hardware. We propose Membrane Potential Alignment (MPA), a test-time adaptation method for spiking neural networks that realigns a pretrained decoder to shifted recordings by only matching membrane potential distributions via KL divergence. By restricting updates to low-rank (LoRA) weights, MPA adapts fewer than 9% of parameters. On a non-human primate reaching task spanning over one month, MPA achieves performance competitive with the state-of-the-art NoMAD method, while using a simpler architecture and finer temporal resolution (4 ms vs. 20 ms). These results show that efficient SNN-based test-time adaptation is a practical path toward long-term, recalibration-free brain-computer interfaces.

Guangzhi Tang

Department of Advanced Computing Sciences, Maastricht University Maastricht, The Netherlands

Test-Time Adaptation Brain-Computer Interfaces Neural Decoding Spiking Neural Networks

Ask PDF

Abstract

Independent Component Analysis

Explore similar work

Neuromotor Hierarchy Network: Physiological Inductive Biases for Robust Generalization in sEMG Decoding

Estimation of Motor Unit Parameters from Surface Electromyograms using an Informed Autoencoder

Test-Time Adaptation of Spiking Neural Networks for Intracortical Neural Decoding using Membrane Potential Alignment