2 Holobar A, Zazula D: Surface EMG Decomposition Research Paper ■ Surface EMG Decomposition Using a Novel Approach for Blind Source Separation Aleš Holobar, Damjan Zazula Abstract. We introduce a new method to perform a blind deconvolution of the surface electromyogram (EMG) signals generated by isometric muscle contractions. The method extracts the information from the raw EMG signals detected only on the skin surface, enabling longtime noninvasive monitoring of the electromuscular properties. Its preliminary results show that surface EMG signals can be used to determine the number of active motor units, the motor unit firing rate and the shape of the average action potential in each motor unit. ■ Infor Med Slov 2003; 8(1): 2-14 Author's institution: Faculty of Electrical Engineering and Computer Science, University of Maribor. Contact person: Aleš Holobar, Faculty of Electrical Engineering and Computer Science, University of Maribor, Smetanova 17, 2000 Maribor. email: ales.holobar@uni-mb.si. Informatica Medica Slovenica 2003; 8(1) 3 Introduction The activity of muscles has been the subject of many studies of bioengineers, physiologists, neurophysiologists, and clinicians for more than 100 years. Many different methods of gathering and interpreting the physiological data and information have been developed. In the past few decades, the assessment of the electrical activity of muscles has proved to be very important. Using the computer based digital processing, many valuable knowledge has been extracted from the electromyographic signals enabling precise medical diagnosing and prevention of possible neurological and muscle disorders1,2. The muscle movement in all human beings is controlled by the central nervous system. It generates the electrical pulses that travel through the motor nerves to different muscles. The neuromuscular junction is called innervation zone and is usually situated about the middle of the muscle body. Each muscle is composed of a large number of tiny muscle fibres, which are organized into so called motor units (MU). Each MU gathers all the fibres that are innervated by the same nerve, i.e. axon. When electrically excited, fibres produce a measurable electrical potential, called action potential (AP), which propagates along the fibres to both directions towards muscle tendons and causes the contraction of the fibres. The electromyograms (EMGs) are taken with different kinds of electrodes whose uptake areas vary according to the electrode size. The majority of the EMG recordings is based on the invasive methods with the needle electrodes3, whose invasive character prevents the long-time monitoring of the electromuscular parameters. The non-invasive method of the surface EMG (SEMG) has been developed recently and has numerous advantages. There is significantly less discomfort, no tissue damage and therefore no subsequent tissue scarring. This allows for unlimited repetition of tests in exactly the same place. Furthermore, recording of SEMG is inexpensive and gives global information about muscle activity1. Finally, linear surface electrode arrays can be used providing additional information about innervation zone location, fibre length and conduction velocity4. The main disadvantage of SEMG is poor morphological information about the MU action potentials (MUAP), caused by different filtering effects of several tissue layers (skin, fat, muscle, etc.)5. In the case of needle electrodes we can selectively observe the action potentials of only a few active motor units, or even of a single muscle fiber3. In the SEMG case, on the other hand, we deal with several tens of active motor units and the measurable contributions from other muscles not being under the clinical investigation, what is often referred to as muscle cross-talk. Many attempts to enhance the MUAP information and to suppress the cross-talk in SEMG were maid in the past as different surface recording technique, such as double differential, Laplacian, etc., were investigated6. Nevertheless, the manual decomposition of SEMG to separate MUAPs is, even with lower muscle contractions, virtually impossible and computer assisted decomposition is required. The non-invasive assessment of muscle properties through the information extracted from the surface EMG signals has introduced new issues also in the field of the EMG signal processing. Despite the numerous efforts7-11, no final technique for SEMG decomposition has been proposed yet. The main problem, as mentioned above, is a very high complexity of the SEMG signals. They are composed of a high number of individual, filtered MUAPs being superimposed into the surface signal. In addition, no a priori information on the number of active motor units and the nature of their mixture is available, either; hence the SEMG signals should be decomposed blindly. The blind separation of the sources has been widely studied in the past and many solutions exist for both instantaneous and convolutive mixtures12. The problem of instantaneous mixtures is most often addressed by exploiting the Actually, this 4 Holobar A, Zazula D: Surface EMG Decomposition possible non-Gaussianity of the sources. is the only possible route when the source signals are independently and identically distributed (i.i.d)13'14. When the first 'i' of 'i.i.d.' is not valid, i.e. when the source signals are correlated in time, another route is to exploit these correlations. Identifiability in such an approach is granted even when the signals are normally distributed, provided the source signals have different spectra15,16. The contributions from other authors17-19 have explored the case where the second 'i' of 'i.i.d.' is failing, that is the non-stationary case. The latter can be successfully applied to the problem of muscular cross-talk20. The problem of convolutive mixtures is much more complex and although many different attempts (wavelets and neural networks classifications, time-frequency decomposition, etc.) 12 were made there is no general solutions for their complete deconvolution. However, the theoretical model of SEMG signals is, under the assumption of stable measurements and isometric muscle contractions, usually based on convolutive mixtures of the nerve train pulses and MUAPs detected on different electrodes. Considering their time-varying nature, the SEMG signals can also be modelled as non-stationary signals. A. Belouchrani and M. Amin18, 19 have addressed general convolutive mixtures of non-stationary signals and exploited the differences of energy locations of sources in time-frequency (t-f) domain. They have proposed to deal with the convolution in the form of mixing matrix by adding delayed repetitions of sources. New sources then form the block diagonal source spatial t-f distribution (STFD) matrices. Hence, with joint block diagonalization24 of observation spatial time-frequency distribution matrices several versions of each source can be retrieved, but only up to a filtering effect19. In this paper, we present preliminary results of a new method for full deconvolution of surface EMG signals. Our approach is based on the work of A. Belouchrani and M. Amin, however, it additionally suggests the construction of diagonal source STFD matrices. These latter matrices are processed into a joint-diagonalization scheme (instead of joint block diagonalization), which provides an estimation of the transfer functions (MUAPs detected on different electrodes) and sources up to the scale factor and the phase shift. Section 2 briefly reviews the algorithms and the problems of joint block diagonalization of spatial t-f distribution matrices for the instantaneous (convolutive) mixtures. In Section 3 we propose a new method for the construction of diagonal spatial t-f distribution matrices in the convolutive case. Section 4 reveals results of the proposed method on the synthetic SEMG signals. We end our paper with the conclusions and discussion in Section 5. Separation of convolutive mixtures in t-f domain Consider a general discrete convolutive multiple-input multiple-output (MIMO), linear, time invariant model given by N L (Os, (t -i)+n (t). (1) j=1 l=0 xi (i = 1,...,M)is one ofMobservations and s, (j = 1,...,N) is one of Nsources (M>N) that are mutually independent for every time/lag and have different structures and localization properties in time-frequency domain. hij is the transfer function between the j-th source and the i-th sensor with the overall extent of (L+1) taps. ni(t) (i=1,...,M ) is additive i.i.d noise, independent from the sources defined by. E[n(t + r)n(t)] = 8(z)<21 m, (2) where E is mathematical expectation, I M the M X M identity matrix, 8(t) the Dirac impulse and N(L+K). The covariance of vectors s(t) and x(t) is then19 Rss (t, t) = E[s (t + t)s (t)] = = diag[R ~i~i(t ,t),..„ r ~n ~n (t, (8) Rxx(t, t) = ARss(t, t)Ah + 5(t) (20) where S(-) is the Dirac impulse, Tj are deterministic, and 0 jn random variables with Gaussian distribution. Again, the sources are assumed to have different localization properties in time-frequency domain, i.e. their pulses should not overlap in time. We will also assume that the average interpulse interval is longer then the length of the impulse response hj, i.e. L in (1). Under these assumptions, the sources correlation matrix R— (0,0) is diagonal. Moreover, the sources have well localized energy in time, which will become advantageous in the process of deconvolution, where we will try to make the STFD matrices D( (t, f) diagonal, too. Diagonal source STFD matrices To preserve a high time resolution, the Wigner-Ville spectra defined by22 should be used as the time-frequency distribution in (13). Using any other distribution (kernel () would spread the energy in time and made the source STFD matrices D( (t, f) block-diagonal. The main disadvantage of using the Wigner-Ville spectra is its high sensitivity to the crossterms (non-zero cross Wigner-Ville spectra WVss (t, f)), which, again, make the source STFD matrices D( (t, f ) block-diagonal. The cross Wigner-Ville spectra in an optional (tk,fk) point is a summation of all the pulses from sources Sj and Sj that fulfil the following relation: 1, t, = -(t +t ), 2 jm/> (22) where tin is the time of appearance of the n-th pulse in source Si, tjm is the time of appearance of the m-th pulse in source Sj, and n, m = As a consequence, the source STFD matrix D^s (tk, fk ) is not diagonal in such (t,, f) point. We can reduce the number of pulse contributions in WVSS(t, f) by calculating pseudo Wigner-Ville spectra22, that is, by limiting T in (21) to the finite interval [-a, a]: PWV (t,f) = (t+1T)Sj (t -T J27fT (23) The number of pulses that contribute to crossterms in PWVsS(t,f) depends on the size of the interval [-a, a] and is finite. Making the limit a small enough, all the crossterms in PWVSS(t, f) are left out, and the source STFD matrices D( (t k , fk ) begin to show their diagonal structure. r=-w z=-a 8 Holobar A, Zazula D: Surface EMG Decomposition The time-frequency plane has now shrunk, because we decreased the resolution of the frequency content in the Wigner-Ville spectra. This made the blind source separation less noise resistant, because the energy of the noise is now less spread along the frequency axis (the time-frequency plane has the property of spreading the noise energy along the frequency axis while localizing the energy of the signal). As a consequence, we have to process longer signals to compensate the noise. However, the source STFD matrices D^ (tk, fk ) are diagonal and we can now use the joint diagonalization 18,21,25-28 instead of joint block diagonalization24 and, hence, retrieve the unfiltered version of the sources up to a scale factor. Selection of diagonal STFD matrices in the case with overlapped pulses Assume certain pulses of the sources overlap. Denote by {li,..., lp} the positions ofp sources in vector of sources S(t) = ),S2(t),...,SN(L+K)(t)]r whose pulses overlap at certain time tk . The source STFD matrix D^ (tk, fk ) will then have p2 non-zero elements at the positions (ki,k2) where k\, k2 g {l\,..., lp } as follows: DÎ (tk, fk ) = Mp, (24) id if k\, k2 g {l\,..., I } MS(k\,k2) = \n u/ p \, (25) 0 where d denotes the energy of one pulse, and the assumption that all pulses have the same energy, like in (20), has been considered. Only p ofp non-zero elements in (25) lie on the diagonal and Mp is far from being diagonal. Selecting the observation STFD matrices at time tk strongly affects the performance of joint diagonalization18,21,25-28. Because U is a unitary matrix, the following relation is valid: eig (Dt (t, f)) = eig (UDt (t, f )U H ) = = eig (Dt (t, f )), (26) where eig(M) denotes the eigenvalues of the matrix M. Noting that the only non-zero eigenvalue of matrix M dp is eig (M dp ) = pd, we can exclude the STFD source matrices that correspond to overlapped source pulses from the process of the joint diagonalization by simply comparing the maximum eigenvalues of the whitened observation STFD matrices D^ (t, f ) . Criteria and algorithms for the selection of (t,f) points in which D^ (t, f ) is diagonal are more in detail described by Holobar et al.28 and by Nguyen et al.29. Results on synthetic SEMG signals In this section, the preliminary results, as investigated via computer simulations, are reported. SEMG signals were generated by SEMG simulator31, that allows to simulate the main features of the surface EMG signal, including the generation and extinction phenomena of the action potentials at the end-plate and tendon regions and the size and shape of the recording electrodes without approximation of the current density source. The simulator models the volume conductor as an anisotropic layered medium with muscle (anisotropic), fat (isotropic) and skin (isotropic) layers. The model allows simulation of multi-channel spatially filtered surface EMG signals and is based on efficient numerical algorithms, which implement the simulation of signals generated during voluntary contractions by the activity of a large number of MUs. The detection systems can be placed either along the fibres direction (usual linear electrode array configuration) or transversally with respect to the Informatica Medica Slovenica 2003; 8(1) 9 muscle fibres. Motor units are placed randomly in the detection volume and are active at the selectable firing rates. In our experiment the surface EMG signals from the biceps brachii muscle during isometric voluntary contraction at low force level were simulated. The main parameters applied in our simulation were the following: 1. MA(4,4) model was chosen, so 4 active MUs were assumed and 6 surface electrodes using the double differential recording technique were simulated what resulted in 4 measurements of SEMG. 2. The length of muscle fibres in all MUs was set of 70 mm. The first MU with 9 fibres was assumed 6 mm, the second with 12 fibres 4 mm, the third with 15 fibres 5 mm, and the fourth with 10 fibres 3 mm deep in the muscle layer. 3. Fibres of the MUs were inclined by 2, 8, 4 and 7 degrees and shifted 7, -5, -2 and 6 mm in the transversal (x in Fig. 1) direction from the centre of the electrode array, respectively. 4. The spread of the innervation zone was taken 6 mm for the first, 6 mm for the second, 5 mm for the third and 5 mm for the fourth MU. 5. The conduction velocity was assumed to be normally distributed with the mean of 4 m/s and standard deviation of 1 m/s (4.54 m/s for the first, 3.83 m/s for the second, 4.33 m/s for the third, and 3.56 m/s for the fourth MU. 6. The thickness of the skin layer was set to 2 mm and thickness of the fat layer to 7 mm. 7. Mean firing rate of the first, second and the third MU were set to 13 Hz, 14 Hz, 13 Hz and 15 Hz with the variance of 3 Hz, respectively. 8. Rectangular electrodes of 5 by 1 mm were simulated with the 5 mm interelectrode distance. 9. The electrode array was assumed placed between the innervation zone and the tendons of fibres. 10. Transversal orientation of detecting system with respect to the muscle fibres was simulated in order to emphasize the differences in contributions (impulse responses) of different motor units to observations, that is to say, to improve the rank of the mixing matrix A. 11. The sampling frequency of 1024 Hz was used for the generated EMG signals. 12. The length of synthetic surface EMG signals was set to 5 s (5120 samples). 13. Signal-to-noise ratio (SNR) was set to 15 dB. The position and orientation of the detection system and the MUs is schematically depicted in Fig. 1, respectively. The generated SEMG signals are partially depicted in Fig. 2. 20 12.2 mm ï9 ir ~ 0 —i-Çïjmm .6 mm 20 1 l I I -50 0 50 100 z axis 4> ¡nervation zone i electrode o tendon f MU radius J number of MU fibers Figure 1 Position and orientation of the detecting system with respect to the simulated active motor units. The MUs are depicted by grey lines ending with circles (tendons), innervation zones by striped rectangular, electrodes by black rectangular. The depth, radius, inclination and the number of fibres in each MU is also depicted. 10 Holobar A, Zazula D: Surface EMG Decomposition Figure 2 Parts of synthetic SEMG signals at four different channels. The detection system was placed transversally with respect to the muscle fibres. The interelectrode distance was set to 5 mm and double differential recording technique was selected. Setting the length of impulse response hj, to 26 samples, 26 estimations of each source were calculated. The estimations were then classified, aligned, normalized, and summed together. They showed almost perfect match with the original sources. Almost all the triggering pulses were successfully recovered. The mean normalized energy of recovered pulses was 0.57 with the variance of 0.21 and the minimum value of 0.14. The ground jitter stayed bellow 0.18, with the mean value of -0.03 and the variance of 0.11. Note the ground jitter is proportional to the nose and exceeds the recovered pulses at the SNR of 8 dB. The results for all 4 train pulses are partly depicted in Fig. 3. 1 0.8 0.6 0.4 0.2 0 ■ jki f5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Time [s] a) ¡0.8 10.6 "O a> ^ 0.4 CO ¡0.2 z 0 0.2 0.4 0.6 Time [s] 0.8 1 b) «,0.8 =3 "SO.6 10.4 CO E0.2 0.2 c) 0.4 0.6 Time [s] d) Figure 3 Comparison of the original innervation pulse trains (black) and corresponding retrieved sources (grey) over a 1s time interval for the first (a), second (b), third (c) and fourth (d) source, respectively. Notice the exact matching of the pulse trains. The recovered MUAPs showed a good match with the original ones, too. The average absolute sample difference between the original MUAPs, detected by different electrodes, and decomposed MUAPs, expressed in percentage to the MAUP Informatica Medica Slovenica 2003; 8(1) 11 amplitude span, was 7.3 %. The recovered MUAPs are partially depicted in Fig. 4. a) c) d) 0.01 0.02 0.03 Time [s] 0.04 0.05 0.01 0.02 0.03 0.04 Time [s] 0.01 0.02 0.03 Time [s] 0.04 0.05 Figure 4 Comparison of the original MUAPs (black traces) and retrieved MUAPs (grey traces); the first MUAP in the third channel (a), the second MUAP in the fourth channel (b), the third MUAP in the second channel (c) and the fourth MUAP in the first channel (d). The impulse responses can be retrieved only up to a scaling factor; the amplitudes of the depicted impulse responses are normalized. Conclusions In this paper, a new approach for blind decovolution of surface EMG signals using pseudo Wigner-Ville time-frequency distributions is introduced. It takes the advantage of non-stationary characteristics of sources (the localization of energy in time) and is based on the joint diagonalization of a combined set of spatial time-frequency (STFD) matrices. A diagonal structure of the source STFD matrices is essential for the proposed approach and is enforced by incorporating only the (t,f) points corresponding to the autoterms (diagonal elements in STFD matrices) of one particular source. The offdiagonal elements in source STFD matrices are crossterms that become zero when calculating the Wigner-Ville spectra on the finite (short enough) interval. The proposed method shows a number of attractive features. The expansion of convolutions to instantaneous mixture makes the STFD matrices very large. As a consequence, the separation can be time and space consuming. Our approach simplifies and fastens the calculation of auto (cross) time-frequency distributions. Moreover, since the time-frequency distributions are shrunk along the frequency axis, the approach is also space efficient. As a result, longer signals can be processed, which compensates the potential noise cancellation due to the effect of spreading the noise power while localizing the source energy in time-frequency domain as a whole. Finally, the K estimations of each source (train of pulses) and its transfer function (MUAPs detected by different electrodes) are retrieved up to a scaling factor by our approach. Hence, the calculated sources and impulse responses can further be improved by averaging the corresponding estimations, which makes the approach more noise resistant. What are the limitations of our approach? Due to the uniqueness property of joint diagonalization18, 21,25-28 and the structure of the searched source STFD matrices (only one non-zero diagonal 12 Holobar A, Zazula D: Surface EMG Decomposition element), we have to find at least one source STFD matrix per source, i.e., each source (including the added delayed repetitions of sources) has to have at least one non-overlapped pulse. This situation is very probable when processing long enough signals and uncorrelated sources with low firing frequencies, i.e. surface EMG signals at low levels of muscle contraction. The performance drops at high contraction forces, due to the effects of motor unit synchronization32. A similar drop in performance is noticed at high firing rates (30Hz and more) when the interpulse intervals of sources become small in comparison to the length of impulse responses. The impulse responses in our model are modelled with constant coefficients. As a consequence, the changes of the shapes of MUAPs in time, such as caused by fatigue, are not recognized by our method. Moreover, the number of active motor units is assumed to be constant. As a consequence, motor unit recruitment during the constant force and muscle contraction of long duration 33, and increasing force, respectively, is not recognized by our method, either. Longer SEMG signals must be broken into subsequent epochs and processed separately. The recommended length of each epoch depends on the muscle type and is generally inversely proportional to the level of muscle contraction. Processing the SEMG signals detected in biceps brachii at 30 % of its maximum voluntary contraction, for example, the recommended length of each epoch is approximately 10 s. Due to possible permutations of source indices, caused by the indeterminacy of blind source separation, the reconstructed train pulses from each epoch may appear in the arbitrary order (two pulse trains, which are reconstructed from two different epochs and share the same index may correspond to different original pulse trains). In order to properly reconstruct the pulse sequences over all epochs, the subsequent epochs must share some common samples (we recommend each pair of subsequent epochs to share a half of common samples, i.e. 10 s long subsequent epochs should share 5 s of common SEMG signal). Aligning the common pulses in recovered innervation trains we easily identify the permutations of source indices and form the whole sequence of triggering pulses for each active MU. The MUAPs must retain constant only through the corresponding epoch, hence, the changes in the shape of MUAPs and the variation of the number of active MUs in time (subsequent epochs) can be monitored, respectively. The analysis of individual motor units is quite important in clinical electromyography. The amplitude and duration of the motor unit action potential provide information on the type of muscle disorder incurred by the peripheral nervous system, the length of time since the disorder's onset, and the evidence for recovery. Furthermore, the reconstruction of the MUAP trains provides information on the firing rate of the individual MUs and on the change of this rate during an increase or decrease of the muscle contraction level. No SEMG decomposition technique has so far provided this information, which makes our approach unique. Acknowledgement This work was supported by the Slovenian Ministry of Education, Science, and Sport (Contract No. S2-796-010/21301/2000), and by the European funding in the 5th Framework project entitled NEW (Contract No. QLK6-2000-00139). References 1. Merletti R: Surface electromyography: possibilities and limitations, J. of Rehab. Sci., Vol. 7, No. 3, 1994, pp. 25-34. 2. Knaflitz M, Balestra G: Computer analyisis of the myoelectric signal, IEEE Micro, No. 10, 1991, pp. 12-58. 3. Farina D, Colombo R, Merletti R, Olsen HB: Evaluation of intra-muscular EMG signal decomposition algorithms, Journal of Electromyography and Kinesiology, No. 11, 2001, pp. 175-187. Informatica Medica Slovenica 2003; 8(1) 13 4. Farina D, Fortunato E, Merletti R: Noninvasive estimation of motor unit conduction velocity distribution using lnear electrode arrays, IEEE Trans. on biomed. eng., Vol 47, No. 3, 2000, pp. 380-388. 5. Farina D,Rainoldi A: Compensation of the effect of sub-cutaneous tissue layers on surface EMG: a simulation study, Med. eng. & Physics, No. 21, 1999, pp. 487-496. 6. Farina D, Merletti R: Effect of electrode shape on spectral features of surface detected motor unit action potentials, Acta Physiol. Pharmacol. Bulg., Vol. 26, pp. 63-66,2001. 7. Zazula D, Korze D, Sostaric A, Korosec D: Study of Methods for Decomposition of Superimposed Signals with Application to Electromyograms, Pedotti et al. (Eds.), Neuroprosthetics, Berlin: Springer Verlag, 1996. 8. Gerber A, Studer RM, de Figueiredo RJP, Moschytz GS: A new framework and computer program for quantitative EMG signal analysis, IEEE Trans. on Biomedical Engineering, Vol. 31, No.12, Dec. 1984, pp.857-863. 9. Zazula D, Sostaric A: Possible Approaches to Surface EMG Decomposition, SENIAM, Vol. 7, 1999, pp. 169-176. 10. Sostaric A, Zazula D, Doncarli C: Time-Scale Decomposition of Compound (EMG) Signal, Electrotechnical Review, Vol. 67, 2000, pp. 69-75. 11. Zazula D, Plevin E: An Approach to Decomposition of Muscle and Nerve Signals, Int. Conf. on Signal, Speech, and Image Processing ICOSSIP '02, Skiathos, Greece, Sept. 2002, CD-ROM. 12. Mansour E, Barros AK, Ohnishi N: Blind separation of sources: methods, assumptions and applications, IEICE trans. Fundamentals, Vol. E83-A, No. 8, 2000, pp. 1498-1512. 13. Cardoso JF: Blind signal separation: statistical principles, Proc. of the IEEE. Special issue on blind identification and estimation, Vol. 9, No. 10, 1998,pp. 2009-2025. 14. Tong L, Liu R, Soon YHVC: Indeterminacy and identifiability of blind identification, IEEE Trans. Circuits Syst., Vol. 38, May 1991, pp. 499-509. 15. Tong L, Liu R: Blind estimation of correlated source signals, Proc. Asilomar Conf., 1990, pp. 161-164. 16. Belouchrani A, Meraim KA, Cardoso JF, Moulines E: A Blind Source Separation Technique Using Second-Order Statistics, IEEE trans. On Signal Porcessing, Vol. 45, No. 2, Feb. 1997, pp. 434-444. 17. Pham DT, Cardoso JF: Blind separation of instantaneous mixtures of non stationary sources, IEEE Transactions on Signal Processing, Vol. 49, No. 9, 2001, pp. 1837-1848. 18. Belouchrani A, Meraim KA: Blind Source Separation based on time-frequency signal representation, IEEE trans. On Signal Porcessing, Vol. 46, No. 11, Nov. 1998, pp. 2888-2898. 19. Boashash B, »Time-Frequency Signal Analysis and Processing«, Prentice Hall PTR, Englewood Cliffs, New Jersey, 2001. 20. Farina D, Merletti R, Indino B, Nazzaro M, Bottin A, Pozzo M, Caruso I: Surface EMG Cross-talk evaluated from experimental recordings and simulated signals, reflections on cross-talk interpretation, quantification and reduction, Proc. Of 4th International Workshop on Biosignal Interpretation BSI 2002, 2002, Como, Italy, pp. 163-166. 21. Cardoso F, Souloumiac A: Jacobi angles for simultaneous diagonalization, SIAM J. Mat. Anal. Appl., Vol. 17, No. 1, 1996, pp. 161-164. 22. Cohen L: Time-Frequency Analysis, Prentice Hall PTR, Englewood Cliffs, New Jersey, 1995 23. Belouchrani A, Cichocki A: A robust whitening procedure in blind source separation context, Electronics Letters, Vol. 36, 2000, pp. 2050-2051. 24. Belouchrani A, Meraim KA, Hua Y: Jacobi-like algorithms for joint block diagonalization: Application to Source Localization, Proc. ISPACS, 1998, Melbourne, Australia. 25. Cardoso JF: A least-squares approach to joint diagonalization, IEEE Signal Processing Lett., Vol. 4, Feb. 1997, pp. 52-53. 26. Cardoso JF, Souloumiac A: Jacobi angles for simultaneous diagonalization, SIAM J. Mat. Anal. Appl., Vol. 17, No. 1, 1996, pp. 161-164. 27. Cardoso JF: Perturbation of joint diagonalizers, technical report Ref\# 94D027, ftp://sig.enst.fr/pub/jfc/ Papers/joint_diag_pert_an.ps, 1994. 28. Holobar A, Fevotte C, Doncarli C, Zazula D: Single autoterm selection for blind source separation in time-frequency plane, Proc. of EUSIPCO, 2002, Toulouse, France, CD-ROM. 29. Nguyen LT, Belouchrani A, Meraim KA, Boashash B: Separating More Sources Than Sensors Using Time-Frequency Distributions, Proc. of 6th ISSPA, Vol. 2 , 2001, Kuala-Lumpur, Malaysia, pp. 583 -586. 30. Merletti R, Lo Conte L, Avignone E, Guglielminotti P: Modelling of surface myoelectric 14 Holobar A, Zazula D: Surface EMG Decomposition signals; Part I: Model implementation, IEEE Trans. Biomed. Engng., Vol. 46, No. 7, pp. 810-820, 1999. 31. Farina D, Merletti R: A novel approach for precise simulation of the EMG signal detected by surface electrodes, IEEE Transactions on Biomedical Engineering, Vol. 48, 2001, pp. 637-646. 32. Weytjens JLF, van Steenberghe D: The effects of motor unit synchronization on the power spectrum of the electromyogram, Biol. Cybern., No. 51, 1984 pp. 71-77. 33. Gazzoni M, Farina D, Merletti R: Motor unit recruitment during constant low force and long duration muscle contractions investigated with surface EMG, IX International Symposium on Motor Control, 2000, Varna, Bulgaria