mathtools
Description¶
This module contains all the mathematical functions that are necessary for the library.
min_max_scaling(series_or_df)
¶
Minmax scaling of pd.series or pd.dataframes.
Minmax feature scaling is often simply referred to as normalisation, which rescales the dataset feature to a range of 0  1. It's calculated by subtracting the feature's minimum value from the value and then dividing it by the difference between the maximum and minimum value.
The formula looks like this: xnorm = x  xmin / xmax  xmin.
PARAMETER  DESCRIPTION 

series_or_df 
The minmax scaling is performed for the entire series, or for single columns in a pd.DataFrame.
TYPE:

RETURNS  DESCRIPTION 

object

The normalised pd.Series or pd.DataFrame (normalised by column).
TYPE:

norm_xcorr(sig1, sig2, out='both')
¶
Normalized crosscorrelation of 2 signals.
PARAMETER  DESCRIPTION 

sig1 
The two signals to correlate. These signals must be 1dimensional and of same length.
TYPE:

sig2 
The two signals to correlate. These signals must be 1dimensional and of same length.
TYPE:

out 
A string indicating the output value:
TYPE:

RETURNS  DESCRIPTION 

xcc

The crosscorrelation value depending on "out".
TYPE:

See also
 norm_twod_xcorr : Normalised 2dimensional crosscorrelation of STAs of two MUS.
norm_twod_xcorr(df1, df2, mode='full')
¶
Normalised 2dimensional crosscorrelation of 2.
The two inputs must have same shape. When this function is used to crosscorrelate MUAPs obtained via STA, df1 and df2 should contain the unpacked STA of the first and second MU, respectively, without np.nan columns.
PARAMETER  DESCRIPTION 

df1 
A pd.DataFrame containing the first 2dimensional signal.
TYPE:

df2 
A pd.DataFrame containing the second 2dimensional signal.
TYPE:

mode 
A string indicating the size of the output:
TYPE:

RETURNS  DESCRIPTION 

normxcorr_df

The results of the normalised 2d crosscorrelation.
TYPE:

normxcorr_max

The maximum value of the 2d crosscorrelation.
TYPE:

See also
 align_by_xcorr : to align the two STAs before calling norm_twod_xcorr.
 unpack_sta : for unpacking the sta dict in a pd.DataFrame before passing it to norm_twod_xcorr.
 pack_sta : for packing the sta pd.DataFrame in a dict where each matrix column corresponds to a dict key.
Examples:
Full steps to pass two dataframes to norm_twod_xcorr from the same EMG file.
 Load the EMG file and bandpass filter the raw EMG signal
 Sort the matrix channels and compute the spiketriggered average
 Extract the STA of the MUs of interest from all the STAs
 Unpack the STAs of single MUs and remove np.nan to pass them to norm_twod_xcorr
 Compute 2dxcorr to identify a common lag/delay
>>> import openhdemg.library as emg
>>> emgfile = emg.askopenfile(filesource="OTB", otb_ext_factor=8)
>>> emgfile = emg.filter_rawemg(emgfile, order=2, lowcut=20, highcut=500)
>>> sorted_rawemg = emg.sort_rawemg(
... emgfile,
... code="GR08MM1305",
... orientation=180,
... )
>>> sta = emg.sta(emgfile, sorted_rawemg, firings=[0, 50], timewindow=100)
>>> mu0 = 0
>>> mu1 = 1
>>> sta_mu1 = sta[mu0]
>>> sta_mu2 = sta[mu1]
>>> df1 = emg.unpack_sta(sta_mu1)
>>> no_nan_sta1 = df1.dropna(axis=1)
>>> df2 = emg.unpack_sta(sta_mu2)
>>> no_nan_sta2 = df2.dropna(axis=1)
>>> normxcorr_df, normxcorr_max = emg.norm_twod_xcorr(
... no_nan_sta1,
... no_nan_sta2,
... )
>>> normxcorr_max
0.7241553627564273
>>> normxcorr_df
0 1 2 125 126
0 0.000002 1.467778e05 3.013564e05 ... 1.052780e06 0.000001
1 0.000004 2.818055e05 6.024427e05 ... 4.452469e06 0.000001
2 0.000007 4.192479e05 9.223725e05 ... 1.549197e05 0.000002
3 0.000009 5.071660e05 1.174545e04 ... 3.078518e05 0.000007
4 0.000007 4.841255e05 1.239106e04 ... 4.232094e05 0.000012
.. ... ... ... ... ... ...
402 0.000005 1.641773e05 3.994943e05 ... 8.170792e07 0.000006
403 0.000001 4.535878e06 1.858700e05 ... 2.087135e06 0.000003
404 0.000004 1.241530e06 5.704194e06 ... 1.027966e05 0.000002
405 0.000004 1.693078e06 1.054646e06 ... 1.811828e05 0.000007
406 0.000002 2.473282e07 6.006046e07 ... 1.605406e05 0.000007
compute_sil(ipts, mupulses, ignore_negative_ipts=False)
¶
Calculate the Silhouette score for a single MU.
The SIL is defined as the difference between the withincluster sums of pointtocentroid distances and the same measure calculated between clusters. The output measure is normalised in a range between 0 and 1.
PARAMETER  DESCRIPTION 

ipts 
The decomposed source (or pulse train, IPTS[mu]) of the MU of interest.
TYPE:

mupulses 
The time of firing (MUPULSES[mu]) of the MU of interest.
TYPE:

ignore_negative_ipts 
If True, only use positive ipts during peak and noise clustering. This is particularly important for sources with large negative components.
TYPE:

RETURNS  DESCRIPTION 

sil

The SIL score.
TYPE:

See also
 compute_pnr : to calculate the Pulse to Noise ratio of a single MU.
Examples:
Calculate the SIL score for the third MU (MU number 2) ignoring the negative component of the decomposed source.
compute_pnr(ipts, mupulses, fsamp, constrain_pulses=[True, 3], separate_paired_firings=True)
¶
Calculate the pulse to noise ratio for a single MU.
PARAMETER  DESCRIPTION 

ipts 
The decomposed source (or pulse train, IPTS[mu]) of the MU of interest.
TYPE:

mupulses 
The time of firing (MUPULSES[mu]) of the MU of interest.
TYPE:

constrain_pulses 
If constrain_pulses[0] == True, the times of firing are considered those in mupulses + the number of samples specified in constrain_pulses[1]. If constrain_pulses[0] == False, the times of firing are estimated via a heuristic penalty funtion (see Notes). constrain_pulses[1] must be an integer (see Notes for instructions on how to set the appropriate value).
TYPE:

separate_paired_firings 
Whether to treat differently paired and nonpaired firings during the estimation of the signal/noise threshold (heuristic penalty funtion, see Notes). This is relevant only if constrain_pulses[0] == False. Otherwise, this argument is ignored.
TYPE:

RETURNS  DESCRIPTION 

pnr

The PNR in decibels.
TYPE:

See also
 compute_sil : to calculate the Silhouette score for a single MU.
Notes
The behaviour of the compute_pnr() function is determined by the argument constrain_pulses.
If constrain_pulses[0] == True, the times of firing are considered those in mupulses + a number of samples specified in constrain_pulses[1]. The inclusion of the samples around the mupulses values allows to capture the full ipts corresponding to the time of firing (e.g., including also the raising and falling wedges). The appropriate value of constrain_pulses[1] must be determined by the user and depends on the sampling frequency. It is suggested to use 3 when the sampling frequency is 2000 or 2048 Hz and increase it if the sampling frequency is higher (e.g. use 6 at 4000 or 4096 Hz). With this approach, the PNR estimation is not related to the variability of the firings.
If constrain_pulses[0] == False, the ipts values are classified as firings
or noise based on a threshold value (named "Pi" or "r") estimated from the
euristic penalty funtion described in Holobar 2012, as proposed in
Holobar 2014. If the variability of the firings is relevant, this
apoproach should be preferred. Specifically:
Pi = D ยท ฯ3,50 + CoVIDI + CoVpIDI
Where:
D is the median of the lowpass filtered instantaneous motor unit discharge
rate (firstorder Butterworth filter, cutoff frequency 3 Hz).
ฯ3,50 stands for an indicator function that penalizes motor units
with filtered discharge rate D below 3 pulses per second (pps) or above
50 pps:
ฯ3,50 = 0 if D is between 3 and 50 or D if D is not between 3 and 50.
Two separate coefficients of variation for interdischarge interval (IDI)
calculated as standard deviation (SD) of IDI divided by the mean IDI,
are used. CoVIDI is the coefficient of variation for IDI of nonpaired
MUs discharges only, whereas CoVpIDI is the coefficient of variation for
IDI of paired MUs discharges.
Holobar 2012 considered MUs discharges paired whenever the second
discharge was within 50 ms of the first.
Paired discharges are typical in pathological tremor and the use of both
CoVIDI and CoVpIDI accounts for this condition.
However, this heuristic penalty function penalizes MUs firing during
specific types of contractions like explosive contractions
(MUs discharge up to 200 pps).
Therefore, in this implementation of the PNR, we did not
include the
penalty based on MUs discharge.
Additionally, the user can decide whether to adopt the two coefficients
of variations to estimate Pi or not.
If both are used, Pi would be calculated as:
Pi = CoVIDI + CoVpIDI
Otherwise, Pi would be calculated as:
Pi = CoV_all_IDI
Examples:
Calculate the PNR value for the third MU (MU number 2) forcing the selction of the times of firing.
>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> mu_of_interest = 2
>>> pnr_value = emg.compute_pnr(
... ipts=emgfile["IPTS"][mu_of_interest],
... mupulses=emgfile["MUPULSES"][mu_of_interest],
... fsamp=emgfile["FSAMP"],
... constrain_pulses=[True, 3],
... )
Calculate the PNR value for the third MU (MU number 2) selecting the times of firing based on the euristic penalty funtion described in Holobar 2012 and considering, separately, the paired and the nonpaired firings.
>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> mu_of_interest = 2
>>> pnr_value = emg.compute_pnr(
... ipts=emgfile["IPTS"][mu_of_interest],
... mupulses=emgfile["MUPULSES"][mu_of_interest],
... fsamp=emgfile["FSAMP"],
... constrain_pulses=[False],
... separate_paired_firings=True,
... )