Assessment of the non-Gaussianity and non-linearity levels of simulated sEMG signals on stationary segments

Accepted Manuscript Assessment of the Non-Gaussianity and Non-Linearity Levels of Simulated sEMG Signals on Stationary Segments Noureddine Messaoudi, Raïs El'hadi Bekka, Philippe Ravier, Rachid Harba PII: DOI: Reference:

S1050-6411(16)30324-8 http://dx.doi.org/10.1016/j.jelekin.2016.12.006 JJEK 2040

To appear in:

Journal of Electromyography and Kinesiology

Received Date: Revised Date: Accepted Date:

27 May 2016 23 November 2016 22 December 2016

Please cite this article as: N. Messaoudi, R.E. Bekka, P. Ravier, R. Harba, Assessment of the Non-Gaussianity and Non-Linearity Levels of Simulated sEMG Signals on Stationary Segments, Journal of Electromyography and Kinesiology (2016), doi: http://dx.doi.org/10.1016/j.jelekin.2016.12.006

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Assessment of the Non-Gaussianity and Non-Linearity Levels of Simulated sEMG Signals on Stationary Segments Noureddine MESSAOUDI (1), (2), Raïs El’hadi BEKKA (2)*, Philippe RAVIER (3), Rachid HARBA(3) (1), (2)

Université de Boumerdès, Faculty of Sciences, Department of Physics, 35000 Boumerdès, Algeria

(2)

Université de Sétif 1, Faculty of Technology, Department of Electronics, LIS Laboratory, 19000 Sétif, Algeria

(3)

Image and Signal Teams of the PRISME Laboratory of Orleans’ University, 12 rue de Blois, 45067 Orleans, France [email protected]

Abstract The purpose of this paper was to evaluate the effects of the longitudinal single differential (LSD), the longitudinal double differential (LDD) and the normal double differential (NDD) spatial filters, the electrode shape, the inter-electrode distance (IED) on nonGaussianity and non-linearity levels of simulated surface EMG (sEMG) signals when the maximum voluntary contraction (MVC) varied from 10% to 100% by a step of 10%. The effects of recruitment range thresholds (RR), the firing rate (FR) strategy and the peak firing rate (PFR) of motor units were also considered. A cylindrical multilayer model of the volume conductor and a model of motor unit (MU) recruitment and firing rate were used to simulate sEMG signals in a pool of 120 MUs for 5 seconds. Firstly, the stationarity of sEMG signals was tested by the runs, the reverse arrangements (RA) and the modified reverse arrangements (MRA) tests. Then the nonGaussianity was characterised with bicoherence and kurtosis, and non-linearity levels was evaluated with linearity test. The kurtosis analysis showed that the sEMG signals detected by the LSD filter were the most Gaussian and those detected by the NDD filter were the least Gaussian. In addition, the sEMG signals detected by the LSD filter were the most linear. For a given filter, the sEMG signals detected by using rectangular electrodes were more Gaussian and more linear than that detected with circular electrodes. Moreover, the sEMG signals are less non-Gaussian and more linear with reverse onion-skin firing rate strategy than those with onion-skin strategy. The levels of sEMG signal Gaussianity and linearity increased with the increase of the IED, RR and PFR. Keywords: Detection systems, Firing rates, Gaussianity, Linearity, Stationarity, Surface EMG signal. 1

1. Introduction The surface techniques for electromyographic (EMG) signal recording have become increasingly important compared to the invasive techniques because they offer some advantages (avoiding discomfort and the risk of infection and providing an indication of the muscle activity in a broad region of space) (Disselhorst-Klug et al., 2000; Farina and Cescon, 2001a). Many spatial filters have been proposed for surface EMG signal detection (Reucher et al., 1987; Disselhorst-Klug et al., 1997; Farina and Cescon, 2001a; Dimitrov et al., 2003; Östlund et al., 2004; Farina et al., 2008; Zhou et al., 2009). The characteristics of the sEMG signals detected by these filters are different. A Gaussian signal can be completely described by its correlation function or its power spectrum. However, there are practical situations as most biomedical signals where the processes are non-Gaussian and non-linear. Therefore, it is necessary to deal them with higher order statistics (HOS) known as cumulants and their associated Fourier transforms of a stationary signal, known as polyspectra (Nikias and Mendel, 1993; Chua et al., 2010). Cumulants and their polyspectra extract information due to deviations from Gaussianity and characterize nonlinear properties in signals (Nikias and Mendel, 1993; Nazarpour et al., 2007; Chua et al., 2010; Orosco et al., 2013). The non-Gaussianity and non-linearity levels of the sEMG signal were addressed experimentally and by computer-simulation approach in some works (Kaplanis et al., 2000; Kaplanis et al., 2009; Zhao and Li, 2012). The test S g of non-Gaussianity is based on the mean bicoherence power. The sEMG linearity test S l is based on the measure of how the bicoherence index is constant in the bi-frequency domain (Kaplanis et al., 2000; Kaplanis et al., 2009; Zhao and Li, 2012). Kaplanis et al. (Kaplanis et al., 2000) used the S g test to characterize the Gaussianity of the sEMG signal recorded from the biceps brachii muscle. They showed that the sEMG signal distribution is highly non-Gaussian at low and high levels of force while having a maximum Gaussianity at 50 % of MVC. A measure of the linearity in (Kaplanis et al., 2000) showed that the EMG signal was less linear at 50% MVC and more linear at 100% MVC which is the reverse pattern with the measure of Gaussianity. In a simulation study on the relation between muscle motor unit numbers and the non-Gaussianity and non-linearity levels of sEMG signals, Zhao and Li (Zhao and Li, 2012) used S g and S l tests and showed that the sEMG signals were more Gaussian and linear when exertion and the 2

number of motor units increased. They concluded that the S g and S l values essentially depended on the number of motor units, the firing rate strategies and their interaction. It was shown in (Zhao and Li, 2012) that the sEMG signal tended towards a more Gaussian distribution when the force increased in contrary to what was reported in (Kaplanis et al., 2000) with increased force.

Kaplanis et al. (Kaplanis et al., 2009) examined the non-

Gaussianity and non-linearity of the sEMG data recorded from 94 subjects for 5 s at 10, 30, 50, 70, 100% MVC. They reported that the surface EMG signal presented a more Gaussian distribution with increase of force up to 70% MVC and was less linear at 70% of MVC compared to 10, 30, 50 and 100% MVC. The level of Non-Gaussianity of the surface EMG signal was measured by other parameters (Nazarpour et al., 2005; Naik et al., 2011; Nazarpour et al., 2013). By using negentropy, Nazarpour et al. (Nazarpour et al., 2005) showed that the probability distribution function (PDF) of sEMG signal was non-Gaussian at small MVC below 30% and it tended to a more Gaussian distribution as MVC level increased. By Kurtosis analysis, Nazarpour et al. (Nazarpour et al., 2013) proved that the sEMG signals had non-Gaussian properties at low contraction levels and tended to a more Gaussian distribution when the force level increased. Naik et al. (Naik et al., 2011) studied the Gaussianity for three MVC levels with kurtosis and concluded that sEMG signal shifted towards the Gaussian with increasing the level of MVC. Earlier, it was found experimentally that the distribution of the sEMG signal recorded during constant-force, constant-angle and non-fatiguing contractions fell between Gaussian and Laplacian densities, with a best fit to the Gaussian density (Clancy and Hogan, 1999). In all studies linked to the evaluation of the non-Gaussianity and non-linearity levels of sEMG signals, the detection system used was always the simple bipolar recording (single differential configuration) (Clancy and Hogan, 1999; Naik et al., 2011; Zhao and Li, 2012; Nazarpour et al. 2005; Nazarpour et al., 2013). Those results showed that the non-Gaussianity level of sEMG signal was significant at low MVC and that its non-linearity level also depended on MVC. Therefore, application of higher order statistics in sEMG signal processing is needed (Nazarpour et al. 2007). But in order to improve the detection quality of the sEMG signal, other more complex spatial filters have been proposed (Reucher et al., 1987; Disselhorst-Klug et al., 1997; Farina and Cescon, 2001a). The differences between them are the configuration of the electrodes within the mask of the spatial filter and the weights given to the electrodes. The effect of the electrode configuration of the different spatial filters on the levels of the 3

non-Gaussianity and non-linearity of the sEMG signals has not been investigated yet and then remain unclear. The objective of this study was to test the non-Gaussianity and non-linearity of stationary surface EMG signals recorded by the LSD, LDD and NDD detection systems when the MVC varied from 10% to 100% by a step of 10% for: i) two recruitment range thresholds (30% and 70% of maximum excitation), and ii) two firing rate strategies (onion skin and reverse onion skin) for each recruitment range thresholds, and iii) three values of the peak firing rate (20, 25 and 30 Hz) for each firing rate strategy. The effects of the electrodes shape, the inter-electrode distance on the non-Gaussianity and non-linearity levels of the sEMG signals were also examined. A suitable window size for analysis of sEMG during contractions was equally investigated using Runs test, Reverse Arrangements (RA) test and Modified Reverse Arrangements (MRA) test. The results of this study were based on the mean values of the bicoherence, kurtosis, and linearity test. The results of this study showed that the level of non-Gaussianity of sEMG signals was not only influenced by the number of motor units and firing rate strategies as in (Zhao and Li, 2012), but it also depended on the used detection system, the electrodes shape, the interelectrode distances, the recruitment range thresholds and the peak firing rates of motor units. 2. Methods 2.1. Surface EMG simulation The simulated surface EMG was implemented in MATLAB according to an EMG model that simulated the surface EMG detected above a elliptical muscle within a cylindrical limb containing bone, muscle, fat, and skin tissues (Farina et al., 2004b), and to model the motor unit (MU) recruitment and firing rate (Fuglevand et al., 1993). The duration of the simulated signal was 5 s at a sampling frequency of 2048 Hz. The model's basic parameters are shown in Table 2.1. 2.1.1. Motor-unit action potential The motor-unit action potential (MUAP) comprises the sum of the action potentials of the muscle fibers belonging to the MU. The stimulated muscle has an elliptical cross section

30 mm, 20 mm (Keenan

et al., 2007) and includes 120 motor units (MUs) uniformly

distributed (Fuglevand et al., 1993; Keenan et al., 2007). The number of muscle fibres (MFs) within MUs is linked to their sizes. The fibres belonging to one MU were uniformly

4

distributed in the MU territory with a density of 20 fibres / mm2 (Fuglevand et al., 1993; Keenan et al., 2007). The diameters of MUs were distributed according to the Poisson law in the range [2  8 mm] with an average of 6 mm (Stashuk, 1993; Wang et al., 2006). The smallest MU innervated 63 fibres and the largest MU innervated 1005 fibres (see table 2.2 for more numerical details). The current density source is proportional to the second derivative of the intracellular action potential (IAP). The simulation of the IAP was based on the analytical description of Rosenfalck (Rosenfalck, 1969). For each muscle fibre, two IAPs are generated at the neuromuscular junction (NMJ), and then they propagate in two opposite directions towards the tendons where they extinguish (Dimitrov and Dimitrova, 1998). The values of conduction velocity of all MUs varied between 2.5 m / s and 5.5 m / s. They were normally distributed in this interval with an average of 4 m / s and a standard deviation of 0.75 m / s (Fuglevand et al., 1993). The simulated signals were detected by LSD, LDD and NDD systems. The electrodes were located halfway between the innervations zone and the distal tendon. The masks of the three investigated spatial filters are given by (Farina et al., 2003): M

LSD



0 0 0

1 0  1 0 , M  LDD 0 0 

0 0 0

1 0  2 0, 1 0 

0 M NDD  1 0

1 0  4 1 1 0 

Two forms of electrodes were considered with LSD and LDD filters. One rectangular of 1 mm wide and 10 mm long and an inter-electrode (IED) distances of 5 and 10 mm. The other was circular of 1 mm diameter and an IED of 5 and 10mm. However, the electrodes configured in the NDD filter were circular of 1 mm electrode diameter and an IED of 5 and 8 mm (Reffad et al., 2014). Table 2.1. Anatomical and physiological parameters of the model. Parameter Muscle cross section Bone (radius and conductivity) Muscle (thick, longitudinal conductivity, transversal conductivity) Fat (thick and conductivity)

Value

30 mm, 20 mm  20 mm

0.02 S / m

Description Elliptic

Reference (Keenan et al., 2007)

Homogeneous and isotropic

26 mm

0.5 S / m

Homogeneous and anisotropic

Farina et al., 2004b

0.1 S / m 3 mm

0.05 S / m

Homogeneous and isotropic,

5

1 mm

Skin (thick and conductivity) Fibre diameter

1S / m 46 m

MFs lengths

[40  160 mm]

MUs diameters

[2  8 mm]

MUs firing rate strategies

FR1(onion skin) and FR2 (reverse onion skin)

Homogeneous and isotropic The same for all fibres Gaussian distribution (mean  80 mm, SD  1 mm) Poisson distribution (mean  6 mm)

(Stashuk, 1993) (Wang et al., 2006 ; Keenan et al., 2007) (Fuglevand et al., 1993; Wang et al., 2006) (Hu et al. 2014)

Gaussian distribution (mean  4 m / s, SD  0.75 m / s)

MUs conduction velocities Recording duration

5 s

MUs territory

-

Circular

MUs axis distribution within muscle

-

Uniform distribution

Number of MUs

120

-

Fibre distribution within the Mus

-

Uniform distribution

Endplates spread

5 mm

Uniform distribution

Tendons regions spread

5 mm

Uniform distribution

MFs density

20 fibres / mm 2

Of the muscle area

ISIs

CV  20 %

Gaussian distribution

The recruitment range thresholds Level of contraction

30 and 70 % of maximum excitation from 10% to 100%

Narrow and broad recruitment range conditions

[2.5  5.5 m / s]

(Fuglevand et al., 1993) (Bilodeau et al., 1997) (Stashuk, 1993 ; Fuglevand et al., 1993; Wang et al., 2006) (Stashuk, 1993 ; Fuglevand et al., 1993; Wang et al., 2006) (Fuglevand et al., 1993) (Farina and Merletti, 2001; Farina et al., 2004; Keenan et al., 2007) (Farina and Merletti, 2001; Keenan et al., 2007) (Farina and Merletti, 2001; Keenan et al., 2007) (Fuglevand et al., 1993; Keenan et al., 2007) (Fuglevand et al., 1993) (Fuglevand et al., 1993) (Zhao and Li, 2012)

Table 2.2. Numerical details of the model at MVC = 100%. (MFs = Muscle fibres). MUs sizes (mm) 2 3 4 5 6 7 8

The number of MUs corresponds to each diameter 6 17 13 31 19 20 14

The number of MFs in each MU 63 141 251 393 565 770 1005

The number of MFs in MUs with the same diameter 378 2397 3263 12183 10735 15400 14070

The number of MFs in the muscle

58426

2.1.2. Motor units recruitment 6

The motor unit is the basic functional element of the central nervous system and muscle that produces movement. It consists of one motor neuron, its axon, and all the muscle fibers that the axon innervates. The group of motor neurons that innervate a single muscle is referred to as a motor pool. The function of a motor neuron is to discharge action potentials and to bring muscle fibers to contract. The force generated by a muscle during contraction depends on the number of motor units that are activated and the firing rates at which the motor neurons discharge action potentials (Enoka, 2006). The recruitment thresholds within a motor unit pool were modeled as an exponential (Fuglevand et al., 1993):

RTE (i)  e

 ln RR   i  n 

(1)

Where i is an index identifying the MU, ln is the natural logarithm, RR is the range of RTE values desired and n is the total number of motor units in the pool. Two RR were selected. A narrow RR was set at 30% (narrow recruitment condition) and wide RR was set at 70 % in a pool of 120 motoneurons. A motoneuron remains inactive as long as excitatory drive is less than its recruitment threshold value (Fuglevand et al., 1993) 2.1.3. Firing rate strategies The train of a motor unit action potential (MUAPT) is the result of the repeated discharge of an alpha motoneuron (Stashuk, 1993). According to the Fuglevand et al. (Fuglevand et al., 1993) models, the repetitive discharge behaviour of motoneurons has four characteristics: 1) the minimum firing rate. 2) The relationship between the excitatory drive and the firing rate. 3) The peak firing rate, and 4) the variability in inter-spike intervals (ISIs) under steady-state excitation. Two types of motor unit firing rate strategies were considered in the simulation (Hu et al., 2014). In the first pattern (FR1, Fig. 1a) termed “onion skin”, the firing rates of earlier recruited MUs (low-threshold) are larger than those of later recruited MUs (high-threshold). In the second strategy (FR2, Fig. 1b) called “reverse onion skin”, the firing rates of earlier recruited MUs (low-threshold) are smaller than those of later recruited MUs. The relationship between the firing rate and the excitatory drive was expressed as:

FRi (t )  g e .( E (t )  RTEi )  MFR

E(t )  RTEi

(2)

7

where FRi is the firing rate and RTEi is the recruitment thresholds of the ith MU, E is the excitatory drive, MFR is the minimum firing rate, and g e is the gain between the excitatory drive and firing rates. The peak firing rate of a MU is inversely proportional to its recruitment threshold according to:

 RTEi  , PFR i  PFR1  PFRD .  RTE120 

i  2 : 120

(3)

where PFR1 is the peak firing rate of the first recruited motor unit, PFRD is the desired difference in peak firing rates between the first and the last recruited motor units, and RTE120 is the recruitment threshold of the last recruited motor unit. In these two strategies, the slope of the relation between the firing rate of the motor unit and the excitatory drive was set to be the same for all motor units. The gain was specified by the equation (Fuglevand et al., 1993):

ge 

PFR i  MFR E max i  RTEi

(4)

where PFR i is the peak firing rate of the ith MU, MFR is the minimum firing rate of a MU (it is the same for all motor units), RTEi is the recruitment threshold of the ith MU and E max i is the level of the excitatory drive for which the ith motor unit bring its PFR. Despite the gain between the motor units firing rates and the excitatory drive has the same mathematical equation with both strategies FR1 and FR2, but the numerical values taken by this gain with the two strategies are not the same. The difference between them resides in the value of the PFR of the last recruited motor unit. Such as, with FR1, the value of the PFR of the last recruited motor unit equal to the value of the PFR of the first recruited motor unit minus the PFRD i.e. it is equal to 20 Hz (Fig. 1a). However with FR2, the value of the PFR of the last recruited motor unit equal to the value of the PFR of the first recruited motor unit plus the PFRD i.e. it is equal to 30Hz (Fig. 1b). PFR varies from 10 Hz up to 35 Hz (De Luca and Hostage, 2010; Hu et al. 2014). In the simulation, PFR1 (case of FR1) and PFR120 (case of FR2) were set at 20, 25 and 30 Hz (Hu et al. 2013), PFRD was set at 10Hz, and RTE120 was set at 30% (narrow recruitment) of maximum excitation and 70% (broad recruitment).

8

The inter-spike intervals (ISI) of the MUs firing was modeled as a Gaussian distribution function using the original model of Fuglevand et al. (Fuglevand et al., 1993). The coefficient of variation (CV) in ISIs was assumed to remain constant. For human muscles during isometric contractions, the CV in ISIs generally ranged between 0.1 and 0.3. In our simulation, CV was set at 0.2. The combination of the three spatial filter configurations, the two electrodes shape, the two IED, seven values (two RR, two FR strategies and three PFR) of MU recruitment and the ten MVC levels allows obtaining 1200 sEMG signals (720 signals with the circular electrode and 480 signals with the rectangular electrode).

Fig.1. The relationship between motor units firing rates and the excitatory drive with the strategies FR1 (a) and FR2 (b).

2.2. Evaluation of stationarity, non-Gaussianity and non-linearity 9

The non-Gaussianity level of sEMG signals was characterized by the bicoherence and the kurtosis, and its non-linearity was evaluated by the S l test. Specifically, in order to apply bicoherence analysis, the stationarity of sEMG signals was first evaluated using the Runs, RA and MRA tests (Bendat et al., 1971; Beck et al., 2006). 2.2.1. sEMG stationarity Although, the long-term sEMG signals are non-stationary, it may be locally stationary. A strictly stationary process is one for which the mean, correlation function, and higher-order moments are invariant to time shifts. In practice, the wide sense stationary (WSS) process is adequate for analysis. In this case the mean and variance are constant and autocorrelation is invariant to time shifts. The study of WSS may be limited to the analysis of mean and the variance rather entire autocorrelation function (Bendat et al., 1971). Therefore, it is first necessary to determine the segment or segments where the sEMG signal is WSS. There are many tests in the literature that can be used to test stationarity of the signals. In this study, the stationarity of the sEMG signals was based on the Runs test, RA test and MRA test which are usually used to examine the stationarity of surface EMG signals (Beck et al., 2006). The WSS test of 1200 sEMG signals (120 signals by MVC) was conducted in accordance to the following protocol: 1) Since the level of the excitatory drive over the initial 1 s is not sufficient to recruit all motor units of the muscle (Fuglevand et al., 1993), the first second was skipped. The stationarity analysis was thus only applied to the last 4 s which is the first segment. 2) This first segment was segmented by one of the four window sizes including 2000, 1000, and 500 and 250 ms to test the local stationarity of each sEMG segment (to test the local stationarity on 2000 ms window, there are 2 adjacent segments of 2000 ms and to test the local stationarity on 1000 ms window there are 4 adjacent segments of 1000 ms, etc.). Consequently, we have thirty one segments. 3) Each of the thirty one segments was divided into 16 equal adjacent sub-segments of duration T ms, and the mean and variance of sub-segments were calculated. The trend of means and variances were tested by the three tests to examine the stationarity of the corresponding segment. For example, to test the stationarity of two segments of 2000 ms (two adjacent time intervals of the signal, each having a duration of 2000 ms), the two segments were divided in 16 sub-segments of 125 ms; a segment was considered stationary only when both parameters have satisfied the three tests for the ten levels of MVC. 10

4) For each window size, the stationarity rate (ST) in % of the sEMG signal was calculated using:

Ns (5) NT N s is the number of stationary segments with a test for each MVC level with respect to a STi % 

window size and N T is the total number of segments with respect to a window size. i is the index of the sEMG signal. 5) In order to compare the stationarity with respect to each MVC level, each spatial filter and each window size by varying the electrodes shape, the IED, the RR, the firing rate strategy and the PFR, the stationarity level was formulated by the mean stationarity rate (MST) in % defined as:

MST % 

N sEMG

ST j %

j 1

N sEMG



(6)

N sEMG is the number of generated signals by each spatial filter for each MVC level ( N sEMG  48 for LSD and LDD, and N sEMG  24 for NDD). 2.7.2. sEMG non-Gaussianity and non-linearity evaluation The non-Gaussianity of a zero-mean, stationary process is quantified by the normalized bispectrum, or bicoherence estimated as (Kaplanis et al., 2009): Bic (1 ,  2 ) 

B(1 ,  2 ) P(1 ) P( 2 ) P(1   2 )

(7)

where P() is the power spectrum. The linearity test ( S l ) involves whether or not the estimated bicoherence is constant in the bifrequency domain, employing a measure of the absolute difference (R) between a theoretical inter-quartile range, Rth and the estimated inter-quartile range, Restim. The non linearity hypothesis was adopted when (Kaplanis et al., 2000): Sl 

R 2 Rth

(8)

The kurtosis is another higher order statistic parameter used to measure the degree of peakedness of PDF. The kurtosis of a process x is defined by (Nazarpour et al., 2013):

Kurt ( x) 

E( x 4 )

E( x ) 2

2

3

(9)

11

where E () denotes the statistical expectation operator. The Kurtosis is 0 for a Gaussian distribution. Kurtosis can be either positive or negative. A process with positive kurtosis is called super-Gaussian or leptokurtic. A process with a negative kurtosis is called sub-Gaussian or platykurtic. 3. Results Fig. 2 shows examples of simulated sEMG signals generated in a cylindrical multilayer volume conductor model and detected by LSD filter when the MVC varied from 10% to 100% by a step of 10% with both firing rate strategies. It appeared clearly that when the level of MVC increases, the amplitude and the frequency of the detected sEMG signal increases. This increase in the EMG signal amplitude and frequency is justified as follows: when the level of MVC increased, the number of the recruited motor units increased which implies an increase in the number of the MUAPTs and as long as the sEMG signal is the sum of the MUAPTs, then an increase in the number of the MUAPTs generated from the recruited motor units implies an increase in the sEMG signal amplitude. At the same time, by the increase of the level of MVC, the firing rates of the recruited motor units increase also which implies an increase of the frequency of the resultant sEMG signal. Parameters of simulation of the surface EMG signals shown in Fig. 2 are given in table 2.1. and table 2.2. FR1 270

MVC = 100%

240

MVC = 90%

Amplitude (V)

210

MVC = 80%

180

MVC = 70%

150

MVC = 60%

120

MVC = 50%

90

MVC = 40%

60

MVC = 30% MVC = 20% MVC = 10%

30 0 0

1

2

3

Time (s)

4

5

(a)

12

FR2 MVC = 100%

300

MVC = 90%

Amplitude (V)

250

MVC = 80%

200

MVC = 70%

150

MVC = 60% MVC = 50%

100

MVC = 40% MVC = 30% MVC = 20% MVC = 10%

50 0 0

1

2

3

4

Time (s)

5

(b)

Fig. 2. Simulated sEMG signals generated in a cylindrical multilayer model and detected by the LSD filter with circular electrodes shape. The diameter of each electrode was set 1mm. The IED was set 5mm. The RR was fixed at 30%. The PFR was set at 25Hz.

To take into account the statistical relevance of the results of stationarity, Gaussianity and linearity, fifteen trials with different start parameters have been made with the combination of the assigned values of the six investigated parameters at each MVC level (1200 sEMG signals were generated in each trial). Then we divided these fifteen trials into three groups containing five trials. The results of stationarity are shown in Fig. 3 to Fig. 5. The results of Gaussianity and linearity are presented in Fig. 6 to Fig. 11. For each parameter of evaluation (stationarity rate in % and mean stationarity rate in % for the results of stationarity, and Kurtosis and bicoherence for results of Gaussianity, and

S l test for the results of

linearity), we computed the mean value of each of these parameters in the same group. However, given the large number of the large number of results to examine the effect of each parameter (three graphs for spatial filter effect (Fig. 6), four graphs for electrodes shape effect (Fig. 7), four graphs for IED effect (Fig. 8), six graphs for RR effect (Fig. 9), six graphs firing rate strategy effect (Fig. 10) and nine graphs for the effect of the PFR (Fig. 11)), only the mean value of the first group was shown. A statistical analysis based on a one-way analysis of variance (ANOVA) was conducted to determine if the results obtained from the three groups are statistically different or not. A one-way ANOVA was used to compare the data means of the three groups. We use hypothesis

H0

which states that all of the data means are equal and the alternative

hypothesis H 1 which states that the 3 data means are not all equal. We retain the null hypothesis

H0

when a p-value “p” is more than 0.05. 13

The results from Fig. 3 to Fig. 11 are presented as mean ± standard error (SE) of the mean and the p-value corresponds to each graph was shown on the legend of each figure. 3.1. sEMG Stationarity Using the Runs test, RA test and MRA test, we investigated the effects of window size and MVC level on the mean stationarity (MST %) of sEMG signals detected by the three spatial filters for two electrodes shape, two IED, two RR, two FR and three PFR. Fig. 3 and 4 present the effects of the window size on the MST % for each MVC level with the runs test (Fig. 3) and the MRA test (Fig. 4). The results show that the mean MST % depends on window size and the MVC level. On average over all MVC levels, 250 ms window produces a higher stationarity level than other window sizes. Other windows lead to good scores > 99% with MRA test. It appears that the MST % were larger with the MRA test than with the Runs test. It should be noted for all window sizes and all MVC levels, the RA test showed that the sEMG signals have higher stationarity level than the two other tests. A one-way ANOVA of the average MST % was not significant (p > 0.05) in the case of runs test (Fig. 3) and MRA test (Fig. 4). LSD filter, p = 0.78 91.84 %

85.51 % 80.21 % 90.93 % 91.41 %

MST %

100 95 90 85 80 75 70 65 60 55 50

250

500

1000

2000

4000

MVC = 10% MVC = 20% MVC = 30% MVC = 40% MVC = 50% MVC = 60% MVC = 70% MVC = 80% MVC = 90% MVC = 100% (a)

Window size (ms)

14

LDD filter, p = 0.59 91.75 %

85.96 % 80.97 % 90.66 % 90.91 %

MST %

100 95 90 85 80 75 70 65 60 55 50

250

500

1000

2000

4000

MVC = 10% MVC = 20% MVC = 30% MVC = 40% MVC = 50% MVC = 60% MVC = 70% MVC = 80% MVC = 90% MVC = 100% (b)

Window size (ms) NDD filter, p = 0.96 100

96.82 %

85.54 %

77 %

93.58 %

1000

2000

95.08 %

95 90 85

MST %

80 75 70 65 60 55 50

250

500

4000

MVC = 10% MVC = 20% MVC = 30% MVC = 40% MVC = 50% MVC = 60% MVC = 70% MVC = 80% MVC = 90% MVC = 100% (c)

Window size (ms)

Fig.3. Effect of the window sizes on the MST % (mean ± SE) with the Run test when the MVC varied from 10% to 100% by a step of 10%. The values shown at the top of the bars of each window size represent the mean of MST % of the ten MVC. LSD filter, p = 0.62 99.25 %

97.68 %

97.27 %

97.2 %

MVC = 10% MVC = 20% MVC = 30% MVC = 40% MVC = 50% MVC = 60% MVC = 70% MVC = 80% MVC = 90% MVC = 100%

MST %

100 99 98 97 96 95 94 93 92 91 90

98.95 %

250

500

1000

2000

4000

(a)

Window size (ms)

15

LDD filter, p = 0.75 100

99.21 %

98.90 %

98.18 %

97.18 %

96.54 %

MVC = 10% MVC = 20% MVC = 30% MVC = 40% MVC = 50% MVC = 60% MVC = 70% MVC = 80% MVC = 90% MVC = 100%

99 98 97

MST %

96 95 94 93 92 91 90

250

500

1000

2000

4000

(b)

Window size (ms)

NDD filter, p = 0.38 98.67 %

98.52 % 98.04 % 97.08 %

MST %

100 99 98 97 96 95 94 93 92 91 90

99.44 %

250

500

1000

2000

4000

MVC = 10% MVC = 20% MVC = 30% MVC = 40% MVC = 50% MVC = 60% MVC = 70% MVC = 80% MVC = 90% MVC = 100% (c)

Window size (ms)

Fig. 4. Effect of the windows sizes on the MST % (mean ± SE) with the MRA test when the MVC varied from 10% to 100% by a step of 10%. The values shown at the top of the bars of each window size represent the mean of MST % of the ten MVC.

3.2. Gaussianity test and linearity test To evaluate the effects of the six investigated parameters on the levels of non-Gaussianity and non-linearity, we first determined the stationarity rate (SR) that corresponds to each parameter and for each window size and this for each trial. For a given window size, a segment was considered to be stationary if and only if it satisfies the three stationarity tests for all parameters studied and all their values. This SR was defined as:

SR 

N SSP N TSP

(10)

For a given window size and for given parameter, N SSP is the number of stationary segments and N TSP is the total number of segments. Fig. 5 shows the effect of the window size on the SR in regards to each of the six investigated parameters. It can be seen that the SR is the largest with 250 ms window size and it is the 16

smallest with 1000 ms window size. Moreover, for a given window size, the highest SR was found with the effect of the IED (except 500 and 1000 ms) and the lowest SR was found with the effect of the PFR. This figure also shows that SR is zero in the following cases: - 1000 ms window size in comparison to the effects of the spatial filter, RR, FR and PFR - 4000 ms window size regarding the effects of the spatial filter and the PFR. A one way ANOVA of the SR corresponding to each one of the six investigated parameters was not significant (p > 0.05) (Fig. 5). On the basis of the results shown in Fig. 3-5, the evaluation of the levels of non-Gaussianity and non linearity of the sEMG signals detected by each spatial filter was based on 250 ms window size. With the same trail, to study the effect of each parameter on the Gaussianity and linearity levels, the analyzed parameter was varied and all other parameters were fixed. This analysis produced a large number of cases. Actually, the number of cases that corresponds to effects of the spatial filter, electrodes shape, IED, RR, firing rate strategy and PFR were 12, 48, 60, 60, 60 and 40, respectively For example, the result of 48 cases for the effect of the electrodes shape means that there are 48 pairs of signals (48 for circular electrodes and 48 for rectangular electrodes). Therefore, the effect of electrodes shape on the levels of sEMG signal Gaussianity and linearity was based on the analysis of 96 signals per MVC level. The same principle was applied to other parameters. Consequently, the study was based on the mean value. For each MVC level and each value of the studied parameter, we calculated the kurtosis, bicoherence and the value of S l test of the sEMG signals for all stationary segments of 250 ms. Then the average values of the kurtosis, bicoherence and S l test were deducted for the comparative study. It is emphasized that the same procedure was repeated to calculate the average values of the kurtosis, bicoherence and S l test for each of the fifteen trials.

17

Filter effect (p = 0.87) Electrodes shape (p = 0.75) IED effect (p = 0.89) RR effect (p = 0.91) FR effect (p = 0.84) PFR effect (p = 0.79)

0,36 0,32 0,28 0,24

SR

0,20 0,16 0,12 0,08 0,04 0,00 0

500 1000 1500 2000 2500 3000 3500 4000

Window size (ms)

Fig. 5. Stationarity rate (mean ± SE) corresponding to the effect of each parameter for each window size.

Fig. 6-11 generally show when the MVC level increases, the average values of the kurtosis (a), bicoherence (b) and S l (c) decrease. 3.2.1. Effects of spatial filter configurations Fig. 6a shows that the average values of the kurtosis of the sEMG signals detected by the NDD filter are larger than those of the sEMG signals detected by the LSD and LDD filters and those of signals detected by the LSD are the lowest. This result reveals that the sEMG signals detected by the LSD filter are more Gaussian than the sEMG signals detected by the LDD and NDD filters. However, the bicoherence analysis (Fig.6b) shows that the average values of the bicoherence are the largest with the LSD filter and they are the smallest with the LDD filter and are thus disagreed with those of kurtosis. This implies that with bicoherence, the sEMG signals detected by the LDD filter are more Gaussian than the sEMG signal detected by the LSD filter. The S l test values indicate that the linearity of the sEMG signals depends on detection system and the MVC level, as shown in Fig. 6c. The signals detected by the LDD and NDD filters are not linear with all MVC levels. The one-way ANOVA measures of the kurtosis (Fig. 6a), the bicoherence (Fig. 6b) and the S l (Fig. 6c) revealed that there was no significant effect (p > 0.05) when the simulation parameters are chosen at different way.

18

54 48 42

18

LSD filter, p = 0.76 LDD filter, p = 0.69 NDD filter, p = 0.66

14

Bicoherence

Kurtosis

36 30 24 18 12

12 10 8 6 4

6 0

LSD filter, p = 0.59 LDD filter, p = 0.66 NDD filter, p = 0.62

16

2

10 20 30 40 50 60 70 80 90 100

Sl

7,0 6,5 6,0 5,5 5,0 4,5 4,0 3,5 3,0 2,5 2,0 1,5 1,0

10 20 30 40 50 60 70 80 90 100

(a)

MVC %

MVC %

(b)

LSD filter, p = 0.77 LDD filter, p = 0.88 NDD filter, p = 0.92

10 20 30 40 50 60 70 80 90 100

MVC %

(c)

Fig. 6. Effect of the spatial filter configurations on the: (a) kurtosis (mean ± SE), (b) bicoherence (mean ± SE) and (c) S l ( mean ± SE) of the detected sEMG signals.

3.2.2. Effect of the electrodes shape Fig. 7a shows that with the LSD and LDD filters and for the ten levels of MVC, the average values of the kurtosis are larger with the circular electrodes shape than with rectangular electrodes. Therefore, the detected signals with rectangular electrodes are more Gaussian than those detected using circular electrodes. Fig.7b shows that the average values of the bicoherence of the sEMG signals detected by the LSD filter are larger than that detected by the LDD filter, regardless of electrodes shape. Moreover, the bicoherence analysis shows that the sEMG signals detected by the LDD filter with rectangular electrodes are less nonGaussian than those detected by the LSD filter with rectangular electrodes. The S l test values indicate that the linearity of the sEMG signals depends on detection system, the electrodes shape and MVC level (Fig. 7c). The signals detected by the LDD filter with rectangular electrode shape are the most linear (from 40 to 100% MVC). It is notable

19

also that the sEMG signals detected by the two spatial filters and both shapes of electrodes become linear from 80% MVC. The one-way ANOVA results of the kurtosis (Fig. 7a), the bicoherence (Fig. 7b) and the

S l test (Fig. 7c) showed that there was not a statistically significant difference between groups (p > 0.05). 27 24

Kurtosis

21

LSD (Circ, p = 0.71) LSD (Rect, p = 0.68) LDD (Circ, p = 0.69) LDD (Rect, p = 0.79)

14

LSD (Circ, P = 0.53) LSD (Rect, p = 0.48) LDD (Circ, p = 0.46) LDD (Rect, p = 0.50)

12

Bicoherence

30

18 15 12 9

10

6

8 6 4 2

3

0

10 20 30 40 50 60 70 80 90 100

MVC % 5,0 4,5 4,0 3,5

10 20 30 40 50 60 70 80 90 100

(a)

MVC %

(b)

LSD (Circ, p = 0.90) LSD (Rect, p = 0.70) LDD (Circ, p = 0.92) LDD (Rect, p = 0.78)

Sl

3,0 2,5 2,0 1,5 1,0 10 20 30 40 50 60 70 80 90 100

MVC %

(c)

Fig.7. Effect of the electrodes shape on the: (a) kurtosis (mean ± SE), (b) bicoherence (mean ± SE) and (c) S l (mean ± SE) of the detected sEMG signals. Here the NDD filter was not considered because it was analyzed only with circular electrodes shape.

3.2.3. Effects of the IED Regarding the average values of the kurtosis (Fig. 8a), the signals detected by LDD filter with

IED  10 mm are more Gaussian than those detected by LDD filter with IED  5 mm. The signals detected by LSD filter with IED  10 mm are slightly more Gaussian than those detected by LSD filter with IED  5 mm. This result implies that for a given filter the level of Gaussianity of the sEMG signals increases with the increase of the IED. In addition, the LSD filter with IED  10 mm allows the detection of the most Gaussian signals. The bicoherence 20

analysis also shows that the level of Gaussianity of the sEMG signals detected by LSD increases with the increase of the IED, but the sEMG signals detected by the LDD filter with both IED are the most Gaussian (Fig. 8b). Fig. 8c shows that the sEMG signals detected by the LDD filter with IED  10mm are not linear only with 10% MVC. For the case of LSD filter with IED  5 mm and IED  10mm and LDD with IED  5mm , the detected sEMG signals become linear from 70% of MVC. For the effect of the IED, the one-way ANOVA measures of the kurtosis (Fig. 8a), the bicoherence (Fig. 8b) and the S l (Fig. 8c) showed there was not significantly different (p > 0.05) between groups.

24 21

Kurtosis

18

LSD (IED = 5 mm, p = 0.63) LSD (IED = 10 mm, p = 0.60) LDD (IED = 5 mm, p = 0.81) LDD (IED = 10 mm, p = 0.81)

15 12 9 6

14 12 10

Bicoherence

27

3

LSD (IED = 5 mm, p = 0.43) LSD (IED = 10 mm, p = 0.47) LDD (IED = 5 mm, p = 0.62) LDD (IED = 10 mm, p = 0.51)

8 6 4 2

10 20 30 40 50 60 70 80 90 100

MVC % 5,0 4,5 4,0

10 20 30 40 50 60 70 80 90 100

(a)

MVC %

(b)

LSD (IED = 5 mm, p = 0.93) LSD (IED = 10 mm, p = 0.76) LDD (IED = 5 mm, p = 0.98) LDD (IED = 10 mm, p = 0.88)

3,5

Sl

3,0 2,5 2,0 1,5 1,0 0,5

10 20 30 40 50 60 70 80 90 100

MVC %

(c)

Fig. 8. Effect of the IED on the: (a) kurtosis (mean ± SE), (b) bicoherence (mean ± SE) and (c) S l (mean ± SE) of the detected sEMG signals. Here the NDD filter was not considered because the values given the IED in the case of LSD and LDD filters are not the same as given in the case of NDD filter.

3.2.4. Effects of the recruitment range thresholds In respect of each investigated filter, Fig. 9a and 9b show that the average values of the kurtosis (Fig. 9a) and bicoherence (Fig. 9b) are larger with RR  70% and smaller with 21

RR  30% . Therefore, the sEMG signals are more Gaussian at narrow MUs recruitment and less Gaussian at broad MUs recruitment. At high levels of MVC, the degree of nonGaussianity of the sEMG signals is similar for both MUs recruitments ranges. For each spatial filter, Fig.9c shows that the average values of S l were smaller with

RR  30% than with RR  70% for many levels of MVC. For the two RR, the sEMG signals detected by the LDD filter become linear from 50% MVC, while those detected by NDD filter are not linear for any level of MVC. For the recruitment range thresholds, the one-way ANOVA results of the kurtosis (Fig. 9a), the bicoherence (Fig. 9b) and the S l (Fig. 9c) were not significantly different with

55 50 45 40 35 30 25 20 15 10 5 0

LSD (RR = 30 %, p = 0.90) LSD (RR = 70 %, p = 0.90) LDD (RR = 30 %, p = 0.82) LDD (RR = 70 %, p = 0.77) NDD (RR = 30 %, p = 0.55) NDD (RR = 70 %, p = 0.72)

16 14 12

Bicoherence

Kurtosis

three groups (p > 0.05).

10

LSD (RR = 30 %, p = 0.35) LSD (RR = 70 %, p = 0.60) LDD (RR = 30 %, p = 0.48) LDD (RR = 70 %, p = 0.70) NDD (RR = 30 %, p = 0.38) NDD (RR = 70 %, p = 0.62)

8 6 4 2

10 20 30 40 50 60 70 80 90 100

10 20 30 40 50 60 70 80 90 100

MVC %

MVC %

(a)

(b)

6,0 5,5

LSD (RR = 30 %, p = 0.90) LSD (RR = 70 %, p = 0.90) LDD (RR = 30 %, p = 0.90) LDD (RR = 70 %, p = 0.77) NDD (RR = 30 %, p = 0.66) NDD (RR = 70 %, p = 0.97)

5,0 4,5

Sl

4,0 3,5 3,0 2,5 2,0 1,5 1,0

10 20 30 40 50 60 70 80 90 100

MVC %

(c)

Fig.9. Effect of the RR on the: (a) kurtosis (mean ± SE), (b) bicoherence (mean ± SE) and (c) S l (mean ± SE) of the detected sEMG signals.

3.2.5. Effects of the firing rate strategy

22

For a given detection system, Fig. 10a and 10b show that the average values of the kurtosis (Fig. 10a) and bicoherence (Fig. 10b) of the sEMG signals simulated under the onion skin firing rate condition are larger than those simulated under the reverse onion skin strategy. This result implies that the reverse onion skin scheme led to more Gaussian signals than the onion skin scheme. The average values of S l show that the reverse onion skin pattern led to more linear sEMG signals than the onion skin scheme for each spatial filter (Fig. 10c). For the firing rate strategy, there was not a statistically significant difference between groups as determined by one-way ANOVA of the kurtosis (Fig. 10a), the bicoherence (Fig. 10b) and the S l test (Fig. 10c) LSD (FR1, p = 0.70) LSD (FR2, p = 0.64) LDD (FR1, p = 0.68) LDD (FR2, p = 0.70) NDD (FR1, p = 0.58) NDD (FR2, p = 0.66)

40 35

Kurtosis

30 25

14

20 15 10 5 0

LSD (FR1, p = 0.50) LSD (FR2, p = 0.48) LDD (FR1, p = 0.43) LDD (FR2, p = 0.48) NDD (FR1, p = 0.46) NDD (FR2, p = 0.52)

12

Bicoherence

45

10 8 6 4 2

10 20 30 40 50 60 70 80 90 100

10 20 30 40 50 60 70 80 90 100

(a)

MVC %

MVC %

(b)

6,0 5,5 5,0 4,5

LSD (FR1, p = 0.92) LSD (FR2, p = 0.82) LDD (FR1, p = 0.93) LDD (FR2, p = 0.89) NDD (FR1, p = 0.91) NDD (FR2, p = 0.89)

Sl

4,0 3,5 3,0 2,5 2,0 1,5 1,0

10

20

30

40

50

60

MVC %

70

80

90 100

(c)

Fig.10. Effect of the firing rate strategy on the: (a) kurtosis (mean ± SE), (b) bicoherence (mean ± SE) and (c) S l (mean ± SE) of the detected sEMG signals.

3.2.6. Effects of the peak firing rate

23

When the PFR increases, the average values of kurtosis (Fig. 11a) and bicoherence (Fig. 11b) decrease whatever the investigated filters. Thus, sEMG signals are more Gaussian at higher PFR. Fig. 11c shows that the linearity of the sEMG signals simulated with three typical values assigned to PFR depend on the filter configuration and the MVC level. For example, the signals detected by the NDD filter under PFR = 20 Hz are not linear for any MVC level. For the effect of the peak firing rate, we accept the null hypothesis H 0 which states that the means of the kurtosis (Fig. 11a), the bicoherence (Fig. 11b) and the S l (Fig. 11c) of

LSD (PFR=20Hz, p=0.32) LSD (PFR=25Hz, p=0.48) LSD (PFR=30Hz, p=0.64) LDD (PFR=20Hz, p=0.29) LDD (PFR=25Hz, p=0.67) LDD (PFR=30Hz, p=0.74) NDD (PFR=20Hz, p=0.42) NDD (PFR=25Hz, p=0.69) NDD (PFR=30Hz, p=0.75)

50 45 40

Kurtosis

35 30 25 20 15 10 5 0

Bicoherence

the three groups are equal (p > 0.05).

15 14 13 12 11 10 9 8 7 6 5 4 3 2 1

LSD (PFR=20Hz, p=0.13) LSD (PFR=25Hz, p=0.40) LSD (PFR=30Hz, p=0.47) LDD (PFR=20Hz, p=0.16) LDD (PFR=25Hz, p=0.67) LDD (PFR=30Hz, p=0.59) NDD (PFR=20Hz, p=0.22) NDD (PFR=25Hz, p=0.62) NDD (PFR=30Hz, p=0.62)

10 20 30 40 50 60 70 80 90 100

10 20 30 40 50 60 70 80 90 100

MVC %

MVC %

(a)

(b)

5,5 5,0

LSD (PFR=20Hz, p=0.80) LSD (PFR=25Hz, p=0.89) LSD (PFR=30Hz, p=0.92) LDD (PFR=20Hz, p=0.34) LDD (PFR=25Hz, p=0.73) LDD (PFR=30Hz, p=0.82) NDD (PFR=20Hz, p=0.48) NDD (PFR=25Hz, p=0.93) NDD (PFR=30Hz, p=0.92)

4,5 4,0

Sl

3,5 3,0 2,5 2,0 1,5 1,0 10 20 30 40 50 60 70 80 90 100

MVC %

(c)

Fig.11. Effect of the peak firing rate on the: (a) kurtosis (mean ± SE), (b) bicoherence (mean ± SE) and (c) S l (mean ± SE) of the detected sEMG signals.

24

4. Discussion The purpose of this study was to test the hypothesis of the non-Gaussianity and non-linearity levels of sEMG signals in accordance with the detection system, the electrodes shape, the inter-electrode distances, the recruitment range thresholds, the firing rate strategy and the peak firing rates of motor units. The non Gaussianity and the non linearity of sEMG signals were evaluated only on stationary segments. The stationarity of the sEMG signals was evaluated using the Runs test, RA test and MRA test. Recent research has reported that surface sEMG signals have non-Gaussian properties and their PDF depends on the level of MVC (Yuan et al. 1997; Kaplanis et al., 2000; Nazarpour et al., 2007; Zhao and Li, 2012; Nazarpour et al., 2013). Moreover, a simulation study on the relation between muscle motor unit numbers and the non-Gaussianity and non-linearity levels of the sEMG signals was investigated in (Zhao and Li, 2012). In all the previously mentioned works, the non-Gaussianity and non-linearity levels of sEMG signals were evaluated without assessing the stationarity of signals. In this work, we began first by determining the segments where the sEMG signal was stationary relive to three tests for the ten levels of MVC. Then, we evaluated their levels of non-Gaussianity and non-linearity. Cho and Kim, (Cho and Kim, 2012) showed that the stationarity level of the EMG signals increased when the window size increased with a relatively low stationarity with window sizes of 250 and 500 ms in the RA test. They showed also that the level of stationary in the MRA test was maximum with a 750 ms window size for any load. The results of this study showed that the stationarity level of the EMG signals was the highest with 250 ms window size and the lowest with 1000 ms window size in the Runs (Fig. 3) and MRA (Fig. 4) tests. These figures also show the level of stationarity of the sEMG signals depended of investigated filter and for a given window size it depended of MVC level. In addition, it was shown that the stationarity rate used to analyze the effect of each one of the six investigated parameters on the non-Gaussianity and non-linearity was the highest with 250 ms window size (Fig. 5). The effects of the six investigated parameters on the levels of non-Gaussianity and nonlinearity have never been treated and were the subject of this study. HOS methods were used for testing the non-Gaussianity of the surface EMG signals (Kaplanis et al., 2000; Zhao and Li, 2012; Nazarpour et al., 2013). Kaplanis et al. (Kaplanis et al., 2000) used bicoherence to quantify the non-Gaussianity level and the linearity test S l to 25

evaluate the non-linearity level of the sEMG signal. They showed that the PDF of sEMG signal is highly non-Gaussian at low and high levels of MVC whereas it is nearly Gaussian at 50 % MVC, and it is more linear at 100% of MVC. Zhao and Li (Zhao and Li, 2012) showed that the non-Gaussianity level of the sEMG signal below 40% of MVC was significant and the linearity level of the sEMG signal above 60% MVC was significant for three firing rate strategies. Our results showed that the nonGaussianity and non-linearity levels of the sEMG signals don't only depend on the number of motor units and firing rate strategies (Zhao and Li, 2012), but also depend on the spatial filter configurations (Fig. 6), the electrodes shape (Fig. 7), the IED (Fig. 8), the recruitment range thresholds (Fig.9), the FR strategy (Fig.10) and the PFR (Fig.11). There is not as yet any work on the characterisation of the non-Gaussianity and non-linearity of sEMG signals detected by various detection systems. In the basis of kurtosis analysis, we showed that the PDF of sEMG signals detected by LSD filter was the most Gaussian and that the PDF of the sEMG detected by the NDD filter was the least Gaussian (Fig. 6a). However, with bicoherence analysis, we found that the sEMG signals detected by the LDD filter were the most Gaussian and that the sEMG signals detected by the LSD filter were the least Gaussian (Fig. 6b). This implies that for a given signal, its gaussianity level depended on the parameter of its evaluation i.e. for the same signal conditions; the evaluation of the Gaussianity and linearity levels of the sEMG signal was different from a one parameter of evaluation to another. It was showed that during isometric contractions, the sEMG signal distribution is more or less super-Gaussian. During low levels of MVC less motor units are recruited and with the increase of the MVC level more motor units begin to fire and making the sEMG signal less super-Gaussian (Naik et al., 2011). This result was confirmed in this study by the values of the kurtosis (from Fig. 6a to Fig. 11a) and the bicoherence (from Fig. 6b to Fig. 11b) of the sEMG signals detected by the LSD, LDD and NDD filters which decrease with the increase of the MVC level. Moreover, in the context of the sEMG signal generation, if the level of the excitatory drive is sufficient to recruit all motor units, at the end of the first second and the beginning of the second, all motor units have began to fire. The last recruited motor unit reaches its PFR after one second i.e. after two seconds of sEMG generation all recruited motor units reach his peak firing rates (Fuglevand et al., 1993). From Fig. 6 to Fig.11, we showed that with low MVC levels (when all motor units were not yet recruited), the values of the Kurtosis and the bicoherence of the sEMG signals detected by the LSD, LDD and NDD 26

filters are far from zero which implies that the sEMG signals were more super-Gaussian. However, when the MVC level was equal or greater than the recruitment range threshold, which was set at 30 and 70% (Fig. 9) of maximum excitation (all motor units were spatially recruited), the values of the kurtosis and the bicoherence became closer to zero, except for the NDD filter where the Kurtosis and bicoherence values were far from zero. This result shows that the increase in the level of MVC shifts the PDF of sEMG signals detected by LSD and LDD towards the Gaussian distribution. Nazarpour et al. (Nazarpour et al., 2013) explored the PDF of sEMG signal by investigating the bicoherence and kurtosis of the sEMG measurement at different isometric muscle contraction levels. They showed that the bicoherence analysis did not help to infer the PDF of measured EMG signals. However, they showed that with kurtosis analysis, the PDF of the EMG signal shift from a non-Gaussian distribution to a more Gaussian-like distribution when the contraction level increased (kurtosis decreased). With kurtosis analysis, we observed trend in the kurtosis values relative to MVC. Indeed, the kurtosis values decreased when the contraction level increased. There are many factors (anatomical, physical and detection system parameters) that may have an influence on the surface EMG signal (Farina et al., 2002) and these factors have automatically an influence on the levels of its non-Gaussianity and non-linearity. Here we concentrated the study only on six parameters (electrodes configuration and shape, IED, RR, FR strategy and PFR). We showed that the effect of each parameter on the non-Gaussianity and non-linearity levels was different of other parameters (Fig. 6 to Fig. 11). The results of stationarity, Gaussianity and linearity showed there was no statistically significant difference between the three groups when the simulations parameters are chosen at different way as determined by one-way ANOVA (p > 0.05). 5. Conclusion The results of this study showed that the levels of non-Gaussianity and nonlinearity of sEMG signals was not only influenced by the number of motor units and firing rate strategies, but also by the used detection system, the electrodes shape, the inter-electrode distances, the recruitment range thresholds, the firing rate strategy and the peak firing rates of motor units. Regarding the average values of the kurtosis, the main conclusion of this study is that the sEMG signal was more Gaussian and more linear with LSD filter, with rectangular electrodes shape, with large inter-electrode distances, with narrow recruitment range thresholds, with reverse onion-skin strategy and with large peak firing rate. 27

References Beck TW, Housh TJ, Weir JP, Cramer JT, Vardaxis V, Johnson GO, Coburn JW, Malek MH, Mielke M. An examination of the Runs Test, Reverse Arrangements Test, and modified Reverse Arrangements Test for assessing surface EMG signal stationarity. Journal of Neuroscience Methods 2006;156:242–48. Bendat J, Enochson LD, Piersol H. Tests for randomness, stationary, normality and comparison of spectra, Random Data: Analysis and Measurement Procedures. New York: Wiley Interscience, 1971. Bilodeau M, Cincera M, Arsenault A, Gravel D. Normality and stationarity of EMG signals of elbow flexor muscles during ramp and step isometric contractions. Journal of Electromyography and Kinesiology 1997;7:87-96. Chua. KC, Chandran V, Acharya UR, Min CL. Application of higher order statistics/spectra in biomedical signals-A review, Medical Engineering & Physics 2010;32:679-89. Clancy EA, Hogan N. Probability Density of the Surface Electromyogram and Its Relation to Amplitude Detectors. IEEE Transactions on Biomedical Engineering 1999;46:730-39. De Luca C J and Hostage E C. Relationship between firing rate and recruitment threshold of motoneurons in voluntary isometric contractions. J. Neurophysiol. 2010;104 :1034–46 Dimitrov GV, Dimitrova NA. Precise and fast calculation of the motor unit potentials detected by a point and rectangular plate electrode. Medical Engineering & Physics 1998;20:37481. Dimitrov GV, Disselhorst-Klug C, Dimitrova NA, Schulte E, Rau G. Simulation analysis of the ability of different types of multi-electrodes to increase selectivity of detection and to reduce cross-talk. J Electromyogr Kinesiol 2003;13:125–38. Disselhorst-Klug C, Silny J, Rau G. Improvement of spatial resolution in surface-EMG: A theoretical and experimental comparison of different spatial filters. IEEE Transactions on Biomedical Engineering 1997;44:567-74. Disselhorst-Klug C, Bahm J, Remaekers V, Trachterna A, Rau G. Non-invasive approach of motor unit recording during muscle contractions in humans. Eur J Appl Physiol 2000;83:144-50. Cho YJ, Kim JY. The effects of load, flexion, twisting and window size on the stationarity of trunk muscle EMG signals. International Journal of Industrial Ergonomics 2012 ; 42:287-92. Enoka RM. Motor Unit. Wiley Encyclopedia of Biomedical Engineering, 2006. 28

Farina D, Merletti R. A novel approach for precise simulation of the EMG signal detected by surface electrodes. IEEE Transactions on Biomedical Engineering 2001b;48:637-46. Farina D, Cescon C, Merletti R. Influence of anatomical, physical, and detection-system parameters on surface EMG. Biological Cybernetics 2002; 86(6):445-56. Farina D, Arendt-Nielsen L, Merletti R, Indino B, Graven-Nielsen T. Selectivity of Spatial Filters for Surface EMG Detection From the Tibialis Anterior Muscle. . IEEE Transactions on Biomedical Engineering 2003;50(3):354-64. Farina D, Mesin L, Martina S. Advances in surface EMG signal simulation with analytical and numerical descriptions of the volume conductor. Medical & Biological Engineering & Computing 2004a;42:467-76. Farina D, Mesin L, Simone M, Merletti R. A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE Transactions on Biomedical Engineering 2004b;51:415-26. Farina D, Kamavuako EN, Wu J, Naddeo F. Entropy-based Optimization of wavelet spatial filters. IEEE Transactions on Biomedical Engineering 2008;55:914-22. Fuglevand AJ, Winter DA, Patla AE. Models of recruitment and rate coding organization in motor-unit pools. J Neurophysiol 1993;70:2470–88. Hu X, Rymer WZ, Suresh NL. Motor unit pool organization examined via spike-triggered averaging of the surface Electromyogram. J. Neurophysiol 20013; 110:1205-1220. Hu X, Rymer WZ, Suresh N. Motor unit firing rate patterns during voluntary muscle force generation: a simulation study. J. Neural Eng 2014; 11:1-9. Kaplanis PA, Pattichis CS, Hadjileontiadis IJ, Panas SM. Bispectral analysis of surface EMG. In: Proceedings of the 10th Mediterranean electrochemical conference 2000;2:770-3. Kaplanis PA, Pattichis CS, Hadjileontiadis LJ, Roberts VC. Surface EMG analysis on normal subjects based on isometric voluntary contraction. Journal of Electromyography and Kinesiology 2009;19:157–71. Keenan KG, Farina D, Meyer FG, Merletti R, Enoka RM. Sensitivity of the cross-correlation between simulated surface EMGs for two muscles to detect motor unit synchronization. Journal of Applied Physiology 2007;102:1193-1201. Keenan KG, Valero-Cuevas. Experimentally valid predictions of muscle force and EMG in models of motor-unit function are most sensitive to neural properties. Journal of Neurophysiology 2007;98:1581-1590.

29

Naik G, Kumar D, Arjunan S. Kurtosis and negentropy investigation of myoelectric signals during different MVCs, In: Proceedings of the BRC, Vitoria, Brazil (2011). Nazarpour K, Sharafat AR, Firoozabadi SM. Negentropy analysis of surface electromyogram signal, in Proc. IEEE SSP, 13th Statistical Signal Processing Workshop, Bordeaux, France, 2005; 974–977 Nazarpour K, Sharafat AR, Firoozabadi SMP. Application of higher order statistics to surface Electromyogram signal classification. IEEE Transactions on Biomedical Engineering 2007;54:1762-9. Nazarpour K, Ali ATH, Guido B, Andrew J. A note on the probability distribution function of the surface electromyogram signal. Brain Research Bulletin 2013;90:88-91. Nikias L and Mendel JM. Signal processing with higher order spectra. IEEE Signal processing magazine 1993; 10 (3): 10-37 Orosco EC, Lopez NM, Sciascio Fdi. Bispectrum-based features classification for myoelectric control. Biomedical Signal Processing and Control 2013;8:153-68. Östlund N, Yu J, Roeleveld K, Karlsson JS. Adaptive spatial filtering of multichannel surface electromyogram signals. Medical & Biological Engineering & Computing 2004;42:82531. Reffad A, Mebarkia K, Vieira TMM, Disselhorst-Klug C. Effect of contraction force and knee joint angle on the spatial representation of soleus activity using high-density surface EMG. Biomed Tech. 2014;59(5):399-411. Reucher H, Rau G, Silny J. Spatial filtering of noninvasive multielectrode EMG: Part IIntroduction to measuring technique and applications. IEEE Transactions on Biomedical Engineering 1987;34:98-05. Rosenfalck P. Intra and extracellular fields of active nerve and muscle fibres: A physicomathematical analysis of different models. Acta. Physiol. Scand 1969;321:1-49. Stashuk DW. Simulation of electromyographic signals. Journal of Electromyography & Kinesiology 1993;3:157-73. Wang W, Stefano Ade, Allen R. A simulation model of the surface EMG signal for analysis of muscle activity during the gait cycle. Computers in Biology and Medicine 2006;36:601-18. Yuan DF, Zhang YT, Herzog W. The probability density function of vibromyographic and electromyographic signals for different levels of contraction of human quadriceps muscles. Proc. Of IEEE 1994;351-52. 30

Zhao Y, LI DX. A simulation study on the relation between muscle motor unit numbers and the non-Gaussianity/non-linearity levels of surface electromyography, Sci China Life Sci, 2012, 55: 958–967, doi: 10.1007/s11427-012-4400-1. Zhou P, Rymer WZ. Factors governing the form of the relation between muscle force and the EMG: A simulation study. Journal of Neurophysiology 2004;92:2878-86. Zhou P, Suresh NL, Lowery MM, Rymer WZ. Nonlinear spatial filtering of multichannel surface electromyogram signals during low force contractions. IEEE Transactions on Biomedical Engineering 2009;56:1871-78.

31

Raïs El’hadi Bekka was born in M’sila, Algeria. He graduated in Electronics Engineering from the Ecole Polytechnique, Alger, 1980. He obtained the Magister degree and the Doctorat d’Etat degree in Electronics from the University of Sétif 1, Algeria, in 1987 and 1994, respectively. Since May 2000, he is Professor at University of Sétif 1. His main research interests are in the areas of signal processing of biomedical signals and multiple-description image coding.

32