decomposition

Description

This module contains the functions needed to extract MU discharge times using convolutive blind source separation techniques.

extend_emg_signal(sig, ext_fact)

Extend an HDsEMG signal by interleaving shifted copies of the channels.

The function generates ext_fact shifted copies of the input signal and interleaves them into a larger channel matrix.

PARAMETER	DESCRIPTION
sig	The EMG signal to extend of shape (n_channels, n_samples). TYPE: ndarray
ext_fact	The extension factor. Determines how many shifted copies of the signal are created. Must be ≥ 1. TYPE: int

RETURNS	DESCRIPTION
e_sig	The extended signal of shape (n_channels * ext_fact, n_samples - ext_fact). Each group of n_channels rows contains a copy of the original signal shifted by 0, 1, ..., ext_fact - 1 samples. TYPE: ndarray

Notes

To re-align the extended signal to the original signal, consider adding a number of zeros equal to 'ext_fact' a the beginning of the signal, or 'ext_fact/2' at the start and end depending on use cases.

Examples:

>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> sig = np.transpose(emgfile["RAW_SIGNAL"].to_numpy())
>>> print(f"sig shape: {np.shape(sig)}")
>>> e_sig = emg.extension(sig, 16)
>>> print(f"e_sig shape: {np.shape(e_sig)}")
sig shape: (64, 66560)
e_sig shape: (1024, 66544)

svd_whitening(e_sig, eigenvalue_percentile)

Perform SVD-based whitening of the extended signal.

This function whitens the input signal by: 1. Computing its covariance matrix. 2. Performing an eigenvalue decomposition. 3. Retaining only those eigenvectors whose eigenvalues exceed a specified percentile threshold. 4. Scaling the retained eigenvectors by the inverse square root of their eigenvalues. 5. Projecting the extended signal onto this whitening basis.

Whitening decorrelates the channels and normalises their variance, producing a transformed signal in which: - Rows correspond to whitened components. - Columns correspond to samples.

PARAMETER	DESCRIPTION
e_sig	Extended signal of shape (n_channels, n_samples). Typically obtained by spatial extension to improve decomposition robustness. TYPE: ndarray
eigenvalue_percentile	Percentile (0-100) used to threshold the eigenvalues of the covariance matrix. Only eigenvectors whose eigenvalues exceed this percentile are retained for whitening. Higher values keep fewer components. If set to 0, all strictly positive eigenvalues are retained. TYPE: float

RETURNS	DESCRIPTION
e_w_sig	Whitened signal array of shape (n_components, n_samples), where n_components equals the number of eigenvalues retained after percentile thresholding. TYPE: ndarray

normc(data)

Column-wise L2 normalisation of a matrix.

PARAMETER	DESCRIPTION
data	Input array of shape (n_rows, n_columns). TYPE: ndarray

RETURNS	DESCRIPTION
ndarray	Array with the same shape as data, where each column has unit L2 norm. Columns with zero norm are left unchanged.

symdecor(data)

Symmetric decorrelation of a square or rectangular matrix.

This orthogonalises the rows of data using the polar decomposition (via SVD), returning the closest orthogonal matrix.

PARAMETER	DESCRIPTION
data	Input array of shape (n_rows, n_columns). TYPE: ndarray

RETURNS	DESCRIPTION
ndarray	Orthogonalised matrix of the same shape as data.

get_cost_function(cost_function_id)

Return a FastICA nonlinearity g(x) and its derivative dg(x) based on the selected identifier.

PARAMETER	DESCRIPTION
cost_function_id	Identifier of the cost function to use: skew g(x) = x² dg(x) = 2x kurtosis g(x) = x³ dg(x) = 3x² log_cosh g(x) = tanh(x) dg(x) = 1 - tanh²(x) TYPE: str

RETURNS	DESCRIPTION
g	FastICA nonlinearity g(x). TYPE: callable
dg	First derivative g'(x). TYPE: callable

RAISES	DESCRIPTION
ValueError	If an invalid cost_function_id is provided.

Notes

The functions returned by this function are the FastICA nonlinearities g(x) and their first derivatives g'(x), rather than the underlying contrast functions G(x), where:

g(x) = G'(x)

The corresponding contrast functions are:

• skew: G(x) = x³ / 3, g(x) = x², and g'(x) = 2x.

• kurtosis: G(x) = x⁴ / 4, g(x) = x³, and g'(x) = 3x².

• log_cosh: G(x) = log(cosh(x)), g(x) = tanh(x), and g'(x) = 1 - tanh²(x).

The term cost function is retained because it is commonly used in MU blind source separation to refer to the nonlinearity employed in the fixed-point optimisation.

source_peaks_classification(ipts, fsamp, init='k-means++', max_drate=None, random_state=42)

Identify MU discharge times from a 1-D source signal.

The function performs: 1. Normalisation of the input source. 2. Peak detection on the normalised signal using a minimum inter-spike interval derived from the maximum physiological discharge rate. 3. K-means clustering (k=2) on the peak amplitudes to separate likely MU discharges (high-amplitude cluster) from noise. 4. Cluster-quality estimation using the silhouette score.

PARAMETER	DESCRIPTION
ipts	One-dimensional source signal (component extracted by the decomposition). TYPE: ndarray
fsamp	Sampling frequency in Hz. TYPE: float
init	Initialisation strategy for k-means. k-means++ selects initial cluster centroids using sampling based on an empirical probability distribution of the points' contribution to the overall inertia. This technique speeds up convergence. The algorithm implemented is "greedy k-means++". It differs from the vanilla k-means++ by making several trials at each sampling step and choosing the best centroid among them. random choose n_clusters observations (rows) at random from data for the initial centroids. TYPE: str {"k-means++", "random"} DEFAULT: "k-means++"
max_drate	Maximum physiological discharge rate in Hz. Used to compute the minimum allowed distance between successive peaks: min_distance_samples = fsamp / max_drate. If None, no minimum spacing constraint is applied. TYPE: float or None DEFAULT: None
random_state	Random seed for k-means reproducibility. TYPE: int DEFAULT: 42

RETURNS	DESCRIPTION
mu_spike_idx	Sample indices of the accepted spikes, corresponding to motor unit discharge times. TYPE: ndarray
sil	The silhouette score indicating how well the "high" cluster is separated from the "low" cluster. Higher values imply better separability. TYPE: float

compute_mu_filter(e_w_sig, mu_mupulses)

Calculate the MU filter.

PARAMETER	DESCRIPTION
e_w_sig	Extended and whitened signal array of shape (n_components, n_samples), where n_components equals the number of eigenvalues retained after percentile thresholding. TYPE: ndarray
mu_mupulses	The times of firing (in samples) of each MU. TYPE: ndarray

RETURNS	DESCRIPTION
mu_filter	Array of shape (n_components,). TYPE: ndarray

ConvolutiveBSSParams dataclass

Parameters for Convolutive Blind Source Separation (BSS).

To be used in the convolutive_bss function.

ATTRIBUTE	DESCRIPTION
n_iterations	Number of decomposition iterations (maximum number of sources to attempt to extract). TYPE: int, default 100
silhouette_threshold	Minimum silhouette threshold for classifying a source as a valid MU. TYPE: float, default 0.9
extension_factor	Factor to extend the signal by. TYPE: int, default 16
cost_function_id	Cost function used by FastICA. Defaults to a skew cost function. See emg.get_cost_function for possible parameters. TYPE: str, default "skew"
rem_activity_index	If True, removes activity around identified spikes from the activity index to reduce re-detection of the same source. TYPE: bool, default False
max_discharge_rate	Maximum allowed discharge rate (in Hz) for spike train classification. If 0, no upper bound is enforced. TYPE: float, default 0
min_spike_count	Minimum number of detected spikes required to accept a motor unit. TYPE: int, default 0
eigenvalue_percentile	Percentile (0-100) of the smallest eigenvalues discarded during whitening. Only eigenvectors whose eigenvalues exceed this percentile are retained for whitening. Higher values keep fewer components. If set to 0, all strictly positive eigenvalues are retained. TYPE: float, default 10
artifact_win_sec	Length (in seconds) of the window removed around high-energy peaks in the activity index. TYPE: float, default 10e-3
tolerance_fastica	Convergence tolerance for the FastICA iteration. TYPE: float, default 1e-4
max_iter_fastica	Maximum number of iterations for FastICA updates for each source. TYPE: int, default 100
cluster_dist_type	Initialisation strategy for spike clustering. TYPE: str {"k-means++", "k-means", "random"}, default "k-means++"
spike_clip_max	Maximum value used to clip the projected source before applying the contrast function. If 0, no clipping is performed. TYPE: float, default 0
random_state	Random seed used for clustering reproducibility. TYPE: int, default 42

Examples:

Use default parameters

>>> import openhdemg.library as emg
>>> params = emg.ConvolutiveBSSParams()
>>> print(params)

Modify the number of iterations and silhouette threshold at params creation

>>> import openhdemg.library as emg
>>> params = emg.ConvolutiveBSSParams(
...     n_iterations=200,
...     silhouette_threshold=0.95,
... )
>>> print(params.n_iterations)

Modify the number of iterations after params creation

>>> import openhdemg.library as emg
>>> params = emg.ConvolutiveBSSParams()
>>> params.n_iterations = 500

convolutive_bss(emgsig, fsamp, decomposition_params=None, **plotting_kwargs)

Decompose HDsEMG signals into motor unit spike trains using convolutive blind source separation (Negro et al., 2016).

PARAMETER	DESCRIPTION
emgsig	Multi-channel EMG signal to decompose of shape (n_channels, n_samples). TYPE: ndarray
fsamp	Sampling frequency of the EMG recording, in Hz. TYPE: float
decomposition_params	Parameters for Convolutive Blind Source Separation (BSS). Please see emg.ConvolutiveBSSParams and the examples for more details and default values. TYPE: dataclass DEFAULT: None

Plotting Parameters

signals : class, optional A class providing Qt Signals for external plotting and progress monitoring obtained from DecompositionSignals. send_only_process : Bool, optional If true, only the signal marking the iteration is sent. This can save data transfer and plotting time. stop_object : class, optional A class providing a stop flag obtained from StopDecompositionObject.

Returns:

mupulses : list of np.ndarray List of arrays, each containing the sample indices (spike times) of one decomposed motor unit. ipts : list of np.ndarray List of arrays, each containing the innervation pulse train of one motor unit. sil : list of float List of silhouette scores for each motor unit.

EMGDecomposer

__init__()

This class provides a convenient interface to run a full decomposition pipeline on an emgfile structure, including optional preprocessing (bandpass and power-line removal), optional removal of bad channels, decomposition using a selected algorithm, and optional post-processing (duplicate MU removal).

Defaults

Decomposition - method: emg.convolutive_bss - parameters: emg.ConvolutiveBSSParams()

Filtering - bandpass: enabled by default (order=2, 20-500 Hz) - power-line removal: disabled by default

Channels - bad-channel exclusion: enabled by default (exclude_bad_channels=True)

Duplicate removal - enabled by default

Examples:

Common part

>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()

Minimal usage (defaults: bandpass enabled, notch disabled, duplicate removal enabled)

>>> decomposer = emg.EMGDecomposer()
>>> decomposed_emgfile = decomposer.run_decomposition(emgfile)

Full customisation of decomposition parameters

>>> decomposer = emg.EMGDecomposer()
>>> params = emg.ConvolutiveBSSParams()
>>> params.n_iterations = 5
>>> params.silhouette_threshold = 0.88
>>> decomposer.set_decomposition_parameters(params=params)
>>> decomposed_emgfile = decomposer.run_decomposition(emgfile)

Quick customisation of decomposition parameters

>>> decomposer = emg.EMGDecomposer()
>>> decomposer.decomposition_parameters.n_iterations = 5
>>> decomposed_emgfile = decomposer.run_decomposition(emgfile)

Disable filtering entirely

>>> decomposer = emg.EMGDecomposer()
>>> decomposer.change_filtering_parameters(
...     bandpass_enabled=False,
...     notch_enabled=False,
... )
>>> decomposed_emgfile = decomposer.run_decomposition(emgfile)

Disable duplicate removal

>>> decomposer = emg.EMGDecomposer()
>>> decomposer.change_duplicate_removal_parameters(
...     duplicate_removal_enabled=False,
... )
>>> decomposed_emgfile = decomposer.run_decomposition(emgfile)

set_decomposition_function(func)

Set the decomposition function to be used by the pipeline.

PARAMETER	DESCRIPTION
func	Decomposition callable. See emg.convolutive_bss for the expected function signature. TYPE: callable

Notes

The function name is stored internally (via func.__name__) and written into the output metadata under emgfile["DECOMPOSITION_PARAMETERS"].

set_decomposition_parameters(params)

Set decomposition parameters.

PARAMETER	DESCRIPTION
params	Dataclass instance containing method-specific parameters. TYPE: dataclass instance

Notes

Parameters are stored and later serialised (via dataclasses.asdict) into emgfile["DECOMPOSITION_PARAMETERS"]. Ensure parameters are JSON-serialisable if you intend to save them to disk.

change_filtering_parameters(bandpass_enabled=True, bandpass_order=2, bandpass_lowcut=20, bandpass_highcut=500, notch_enabled=False, notch_freq=50.0, notch_width=5.0)

Configure preprocessing filters.

PARAMETER	DESCRIPTION
bandpass_enabled	Whether to bandpass filter the signal. TYPE: bool DEFAULT: True
bandpass_order	The filter order. Note that a band-pass transformation doubles the order, so order=2 produces a 4th-order band-pass filter. TYPE: int DEFAULT: 2
bandpass_lowcut	The lower cut-off frequency in Hz. TYPE: float DEFAULT: 20
bandpass_highcut	The higher cut-off frequency in Hz. TYPE: float DEFAULT: 500
notch_enabled	Whether to apply notch filtering. TYPE: bool DEFAULT: False
notch_freq	Fundamental power-line frequency (e.g., 50 or 60 Hz). TYPE: float DEFAULT: 50.0
notch_width	Width of each notch (± notch_width/2), in Hz. TYPE: float DEFAULT: 5.0

get_filtering_parameters()

Return the current filtering parameters.

RETURNS	DESCRIPTION
params	Dictionary with keys: - "bandpass_enabled" - "bandpass_order" - "bandpass_lowcut" - "bandpass_highcut" - "notch_enabled" - "notch_freq" - "notch_width" TYPE: dict

Notes

This method is intended for inspection. To modify filtering parameters, use change_filtering_parameters().

use_good_channels_only(ans=True)

Enable or disable exclusion of bad channels.

PARAMETER	DESCRIPTION
ans	If True, channels marked as False in emgfile["GOOD_CHANNELS"] are excluded prior to decomposition. TYPE: bool DEFAULT: True

Notes

If enabled but emgfile["GOOD_CHANNELS"] is missing, a warning is issued and no channels are removed.

change_duplicate_removal_parameters(duplicate_removal_enabled=True, correlation_max_lag=0.05, peak_window_half_width=0.0025, duplicate_threshold=30, which='accuracy')

Configure post-processing removal of duplicate MUs.

PARAMETER	DESCRIPTION
duplicate_removal_enabled	Whether to remove duplicate units. TYPE: bool DEFAULT: True
correlation_max_lag	Maximum lag (in seconds) used when computing the cross-correlation between MU spike trains. Defines the full search range for possible synchronisation peaks. Larger values allow detection of synchrony over wider time shifts. This must be < 1. TYPE: float DEFAULT: 50e-3
peak_window_half_width	Half-width (in seconds) of the window used around the cross-correlation peak to compute the duplication sensitivity metric. This window should capture the narrow temporal jitter expected for duplicate MUs. This must be < 1 and < correlation_max_lag. TYPE: float DEFAULT: 2.5e-3
duplicate_threshold	Threshold (in percent) for classifying two MUs as duplicates. The sensitivity metric is computed as the sum of the cross-correlation values within a ±peak_window_half_width window around the correlation peak, normalised by the size of the larger spike train. TYPE: float DEFAULT: 30
which	How to remove the duplicated MUs. accuracy The MU with the lowest accuracy is removed. The emgfile must already contain the 'ACCURACY' dataframe. covisi The MU with the highest CoV of interspike interval is removed. TYPE: str {"accuracy", "covisi"} DEFAULT: "accuracy"

get_duplicate_removal_parameters()

Return the current duplicate-removal parameters.

RETURNS	DESCRIPTION
params	Dictionary with keys: - "duplicate_removal_enabled" - "correlation_max_lag" - "peak_window_half_width" - "duplicate_threshold" - "which" TYPE: dict

Notes

This method is intended for inspection. To modify duplicate-removal parameters, use change_duplicate_removal_parameters().

run_decomposition(emgfile, **plotting_kwargs)

Run the full decomposition pipeline on an emgfile.

The pipeline performs: 1) Optional bandpass filtering. 2) Optional power-line harmonics removal. 3) Optional exclusion of bad channels (if GOOD_CHANNELS is present). 4) Decomposition using the configured method. 5) Reconstruction of the EMG file with decomposition outputs. 6) Optional removal of duplicate MUs.

PARAMETER	DESCRIPTION
emgfile	emgfile structure. Must provide at least: - "RAW_SIGNAL": pd.DataFrame of shape (n_samples, n_channels) - "FSAMP": sampling frequency (Hz) Optionally: - "GOOD_CHANNELS": dict mapping channel identifiers to booleans - "EMG_LENGTH": expected number of samples (validated if present) TYPE: dict - like

Plotting Parameters

signals : class, optional A class providing Qt Signals for external plotting and progress monitoring obtained from DecompositionSignals. send_only_process : bool, optional If true, only the signal marking the iteration is sent. This can save data transfer and plotting time. stop_object : class, optional A class providing a stop flag obtained from StopDecompositionObject.

RETURNS	DESCRIPTION
decomposed_emgfile	A deep copy of the input EMG file, updated with decomposition results. Keys added/updated include: - "DECOMPOSITION_PARAMETERS" - "NUMBER_OF_MUS" - "MUPULSES" (when MUs are found) - "IPTS" (when MUs are found) - "ACCURACY" (silhouette scores, when MUs are found) - "BINARY_MUS_FIRING" (when MUs are found) TYPE: dict - like

RAISES	DESCRIPTION
ValueError	If EMG_LENGTH is present and inconsistent with the length of "RAW_SIGNAL".