emgdecompy.decomposition

Module Contents

Functions

initial_w_matrix(z, l=31)

Find highest activity regions of z to use as initializations of w.

deflate(w, B)

w = w - BB^{T} * w

gram_schmidt(w, B)

Gram-Schmidt orthogonalization.

orthogonalize(w, B, fun=gram_schmidt)

Performs orthogonalization using selected orthogonalization function.

normalize(w)

Normalize the input vector (scale the elements of the vector so its length is 1).

separation(z, w_init, B, Tolx=0.001, contrast_fun=skew, ortho_fun=gram_schmidt, max_iter=10, verbose=False)

Finds the separation vector for the i-th source using latent component analysis

silhouette_score(s_i, peak_indices)

Calculates silhouette score on the estimated source.

pnr(s_i, peak_indices)

Returns pulse-to-noise ratio of an estimated source.

refinement(w_i, z, i, l=31, sil_pnr=True, thresh=0.9, max_iter=10, random_seed=None, verbose=False)

Refines the estimated separation vectors determined by the separation function

decomposition(x, discard=None, R=16, M=64, bandpass=True, lowcut=10, highcut=900, fs=2048, order=6, Tolx=0.001, contrast_fun=skew, ortho_fun=gram_schmidt, max_iter_sep=10, l=31, sil_pnr=True, thresh=0.9, max_iter_ref=10, random_seed=None, verbose=False)

Blind source separation algorithm that utilizes the functions

emgdecompy.decomposition.initial_w_matrix(z, l=31)[source]

Find highest activity regions of z to use as initializations of w. Highest activity regions of z refers to the time instances corresponding to the highest values in the squared summation of all the whitened and extended observation vectors. Used for step 1 in Negro et al. 2016.

Parameters
  • z (numpy.ndarray) – The whitened extended observation matrix. shape = M*(R+1) x K M = number of channels R = extension factor K = number of time points

  • l (int) – Required minimal horizontal distance between peaks. Default value of 31 samples is approximately equivalent to 15 ms at a 2048 Hz sampling rate.

Returns

  • numpy.ndarray – Peak indices for high activity columns of z.

  • numpy.ndarray – Corresponding peak heights for each peak index.

Examples

>>> initial_w_matrix(z)
emgdecompy.decomposition.deflate(w, B)[source]

w = w - BB^{T} * w Note: this is not true orthogonalization, such as the Gram–Schmidt process. This is dubbed in Negro et al. (2016) as the “source deflation procedure.”

Parameters
  • w (numpy.ndarray) – Vector we are “orthogonalizing” against columns of B.

  • B (numpy.ndarray) – Matrix of vectors to “orthogonalize” w by. Should contain float dtype.

Returns

‘Deflated’ w.

Return type

numpy.ndarray

Examples

>>> w = np.array([7, 4, 6])
>>> B = np.array([[ 1. ,  1.2,  0. ],
                  [ 2. , -0.6,  0. ],
                  [ 0. ,  0. ,  0. ]])
>>> deflate(w, B)
array([-28. ,  -3.2,   6. ])
emgdecompy.decomposition.gram_schmidt(w, B)[source]

Gram-Schmidt orthogonalization.

Parameters
  • w (numpy.ndarray) – Vector we are orthogonalizing against columns of B.

  • B (numpy.ndarray) – Matrix of vectors to orthogonalize w by. Should contain float dtype.

Returns

Orthogonalized w.

Return type

numpy.ndarray

Examples

>>> w = np.array([7, 4, 6])
>>> B = np.array([[ 1. ,  1.2,  0. ],
                  [ 2. , -0.6,  0. ],
                  [ 0. ,  0. ,  0. ]])
>>> gram_schmidt(w, B)
array([0., 0., 6.])
emgdecompy.decomposition.orthogonalize(w, B, fun=gram_schmidt)[source]

Performs orthogonalization using selected orthogonalization function.

Parameters

w: numpy.ndarray

Vector we are orthogonalizing against columns of B.

B: numpy.ndarray

Matrix of vectors to orthogonalize w by. Should contain float dtype.

fun: function

What function to use for orthogonalizing process. Current options are:

  • gram_schmidt (default)

  • deflate

Returns

Orthogonalized w.

Return type

numpy.ndarray

Examples

>>> w = np.array([7, 4, 6])
>>> B = np.array([[ 1. ,  1.2,  0. ],
                  [ 2. , -0.6,  0. ],
                  [ 0. ,  0. ,  0. ]])
>>> orthogonalize(w, B)
array([0., 0., 6.])
emgdecompy.decomposition.normalize(w)[source]

Normalize the input vector (scale the elements of the vector so its length is 1). This is done using the formula w/||w||.

Parameters

w (numpy.ndarray) – Vector to normalize.

Returns

Normalized vector.

Return type

numpy.ndarray

Examples

>>> w = np.array([5, 6, 23, 29])
>>> normalize(w)
array([0.13217526, 0.15861032, 0.60800622, 0.76661653])
emgdecompy.decomposition.separation(z, w_init, B, Tolx=0.001, contrast_fun=skew, ortho_fun=gram_schmidt, max_iter=10, verbose=False)[source]

Finds the separation vector for the i-th source using latent component analysis that maximizes for sparsity. Implemented with a fixed point algorithm. Step 2 in Negro et al.(2016).

Parameters
  • z (numpy.ndarray) – Extended and whitened observation matrix.

  • w_init (numpy.ndarray) – Initial separation vector.

  • B (numpy.ndarray) – Current separation matrix.

  • Tolx (numpy.ndarray) – Tolx for element-wise comparison.

  • contrast_fun (function) – Contrast function to use. skew, log_cosh or exp_sq

  • ortho_fun (function) – Orthogonalization function to use. gram_schmidt or deflate or None

  • max_iter (int > 0) – Maximum iterations for fixed point algorithm. When to stop if it doesn’t converge.

  • verbose (bool) – If true, print fixed-point algorithm iterations.

Returns

Estimated separation vector for the i-th source.

Return type

numpy.ndarray

Examples

>>> w_i = separation(z, w_init, B)
emgdecompy.decomposition.silhouette_score(s_i, peak_indices)[source]

Calculates silhouette score on the estimated source.

Defined as the difference between within-cluster sums of point-to-centroid distances and between-cluster sums of point-to-centroid distances. Measure is normalized by dividing by the maximum of these two values (Negro et al. 2016).

Parameters
  • s_i (numpy.ndarray) – Estimated source. 1D array containing K elements, where K is the number of samples.

  • peak_indices_a (numpy.ndarray) – 1D array containing the peak indices.

Returns

Silhouette score.

Return type

float

Examples

>>> s_i = np.array([0.80749775, 10, 0.49259282, 0.88726069, 5,
                    0.86282998, 3, 0.79388539, 0.29092294, 2])
>>> peak_indices = np.array([1, 4, 6, 9])
>>> silhouette_score(s_i, peak_indices)
0.740430148513959
emgdecompy.decomposition.pnr(s_i, peak_indices)[source]

Returns pulse-to-noise ratio of an estimated source.

Parameters
  • s_i (numpy.ndarray) – Square of estimated source. 1D array containing K elements, where K is the number of samples.

  • peak_indices (numpy.ndarray) – 1D array containing the peak indices.

Returns

Pulse-to-noise ratio.

Return type

float

Examples

>>> s_i = np.array([0.80749775, 10, 0.49259282, 0.88726069, 5,
                    0.86282998, 3, 0.79388539, 0.29092294, 2])
>>> peak_indices = np.array([1, 4, 6, 9])
>>> pnr(s_i, peak_indices)
8.606468362838562
emgdecompy.decomposition.refinement(w_i, z, i, l=31, sil_pnr=True, thresh=0.9, max_iter=10, random_seed=None, verbose=False)[source]

Refines the estimated separation vectors determined by the separation function as described in Negro et al. (2016). Uses a peak-finding algorithm combined with K-Means clustering to determine the motor unit spike train. Updates the estimated separation vector accordingly until regularity of the spike train is maximized. Steps 4, 5, and 6 in Negro et al. (2016).

Parameters
  • w_i (numpy.ndarray) – Current separation vector to refine.

  • z (numpy.ndarray) – Centred, extended, and whitened EMG data.

  • i (int) – Decomposition iteration number.

  • l (int) – Required minimal horizontal distance between peaks in peak-finding algorithm. Default value of 31 samples is approximately equivalent to 15 ms at a 2048 Hz sampling rate.

  • sil_pnr (bool) – Whether to use SIL or PNR as acceptance criterion. Default value of True uses SIL.

  • thresh (float) – SIL/PNR threshold for accepting a separation vector.

  • max_iter (int > 0) – Maximum iterations for refinement.

  • random_seed (int) – Used to initialize the pseudo-random processes in the function.

  • verbose (bool) – If true, refinement information is printed.

Returns

  • numpy.ndarray – Separation vector if SIL/PNR is above threshold. Otherwise return empty vector.

  • numpy.ndarray – Estimated source obtained from dot product of separation vector and z. Empty array if separation vector not accepted.

  • numpy.ndarray – Peak indices for peaks in cluster “a” of the squared estimated source. Empty array if separation vector not accepted.

  • float – Silhouette score if SIL/PNR is above threshold. Otherwise return 0.

  • float – Pulse-to-noise ratio if SIL/PNR is above threshold. Otherwise return 0.

Examples

>>> w_i = refinement(w_i, z, i)
emgdecompy.decomposition.decomposition(x, discard=None, R=16, M=64, bandpass=True, lowcut=10, highcut=900, fs=2048, order=6, Tolx=0.001, contrast_fun=skew, ortho_fun=gram_schmidt, max_iter_sep=10, l=31, sil_pnr=True, thresh=0.9, max_iter_ref=10, random_seed=None, verbose=False)[source]

Blind source separation algorithm that utilizes the functions in EMGdecomPy to decompose raw EMG data. Runs data pre-processing, separation, and refinement steps to extract individual motor unit activity from EMG data. Runs steps 1 through 6 in Negro et al. (2016).

Parameters
  • x (numpy.ndarray) – Raw EMG signal.

  • discard (slice, int, or array of ints) – Indices of channels to discard.

  • R (int) – How far to extend x.

  • M (int) – Number of iterations to run decomposition for.

  • bandpass (bool) – Whether to band-pass filter the raw EMG signal or not.

  • lowcut (float) – Lower range of band-pass filter.

  • highcut (float) – Upper range of band-pass filter.

  • fs (float) – Sampling frequency in Hz.

  • order (int) – Order of band-pass filter.

  • Tolx (float) – Tolerance for element-wise comparison in separation.

  • contrast_fun (function) – Contrast function to use. skew, og_cosh or exp_sq

  • ortho_fun (function) – Orthogonalization function to use. gram_schmidt or deflate

  • max_iter_sep (int > 0) – Maximum iterations for fixed point algorithm.

  • l (int) – Required minimal horizontal distance between peaks in peak-finding algorithm. Default value of 31 samples is approximately equivalent to 15 ms at a 2048 Hz sampling rate.

  • sil_pnr (bool) – Whether to use SIL or PNR as acceptance criterion. Default value of True uses SIL.

  • thresh (float) – SIL/PNR threshold for accepting a separation vector.

  • max_iter_ref (int > 0) – Maximum iterations for refinement.

  • random_seed (int) – Used to initialize the pseudo-random processes in the function.

  • verbose (bool) – If true, decomposition information is printed.

Returns

Dictionary containing:
B: numpy.ndarray

Matrix whose columns contain the accepted separation vectors.

MUPulses: numpy.ndarray

Firing indices for each motor unit.

SIL: numpy.ndarray

Corresponding silhouette scores for each accepted source.

PNR: numpy.ndarray

Corresponding pulse-to-noise ratio for each accepted source.

Return type

dict

Examples

>>> gl_10 = loadmat('../data/raw/gl_10.mat')
>>> x = gl_10['SIG']
>>> decomposition(x)