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:
Series
|
mupulses
|
The time of firing (MUPULSES[mu]) of the MU of interest.
TYPE:
ndarray
|
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:
list
DEFAULT:
[True, 3]
|
separate_paired_firings
|
Whether to treat differently paired and non-paired 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:
bool
DEFAULT:
False
|
| RETURNS |
DESCRIPTION |
pnr
|
TYPE:
float
|
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 low-pass filtered instantaneous motor unit discharge
rate (first-order Butterworth filter, cut-off 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 inter-discharge 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 non-paired
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 non-paired 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,
... )