Autoregressive and cepstral analyses of motor unit action potentials

Medical Engineering & Physics 21 (1999) 405–419 www.elsevier.com/locate/medengphy

Autoregressive and cepstral analyses of motor unit action potentials Constantinos S. Pattichis a

a,*

, Andreas G. Elia

b

Department of Computer Science, University of Cyprus, Kallipoleos 75, P.O. Box 20537, CY1678, Nicosia, Cyprus b The Cyprus Institute of Neurology and Genetics, P.O. Box 23462, CY1683, Nicosia, Cyprus Received 10 March 1999; accepted 20 July 1999

Abstract Quantitative electromyographic signal analysis in the time domain for motor unit action potential (MUAP) classification and disease identification has been well documented over recent years. Considerable work has also been carried out in the frequency domain using classical power spectrum analysis techniques. Although MUAP autoregressive (AR) spectral analysis has been suggested as a diagnostic tool by a number of studies, it has not been thoroughly investigated yet. This work investigates the application of AR modeling and cepstral analysis for the diagnostic assessment of MUAPs recorded from normal (NOR) subjects and subjects suffering with motor neuron disease (MND) and myopathy (MYO). The following feature sets were extracted from the MUAP signal: (i) time domain measures, (ii) AR spectral measures, (iii) AR coefficients, and (iv) cepstral coefficients. Discriminative analysis of the individual features was carried out using the univariate and multiple covariance analysis methods. Both methods showed that: (i) the duration measure is the best discriminative feature among the time domain parameters, and (ii) the median frequency is the best discriminative feature among the AR spectral measures. Furthermore, the classification performance of the above four feature sets was investigated for three classes (NOR, MND and MYO) using artificial neural networks. Results showed that the highest diagnostic yield was obtained with the time domain measures followed by the cepstral coefficients, the AR spectral parameters, and the AR coefficients. In conclusion, MUAP autoregressive and cepstral analyses combined with time domain analysis provide useful information in the assessment of myopathology.  1999 IPEM. Published by Elsevier Science Ltd. All rights reserved. Keywords: Electromyography; Motor Unit Action Potential; Autoregressive analysis; Cepstral analysis; Neural networks

1. Introduction Electromyography (EMG) is the study of the electrical activity of muscle and forms a valuable aid in the assessment of neuromuscular disorders. EMG findings are used to detect and describe different disease processes affecting the motor unit, the smallest functional unit of the muscle. At slight voluntary muscle contraction the motor unit action potential is recorded (MUAP) that reflects the electrical activity of a single anatomical motor unit. In the early days of EMG, neurophysiologists used to assess MUAPs from their shape by watching the oscilloscope and hearing their audio characteristics. With these simple means, it was possible for an experienced electrophysiologist to detect abnormalities with reasonable accuracy. However, MUAP assessment although satis-

* Corresponding author. Tel.: +357-2-892240; fax: +357-2-339062. E-mail address: [email protected] (C.S. Pattichis)

factory for the detection of obvious abnormalities was not adequate for less apparent deviations or mixed patterns of abnormalities. These ambiguous cases called for quantitative MUAP analysis. Buchthal and coworkers [1–3] introduced manual methods for quantitative MUAP analysis where MUAPs were recorded photographically and then selected for analysis. Manual measurement of time domain parameters like duration, amplitude and number of phases was carried out. Twenty potentials were collected from each muscle, and their data was then compared to values obtained from the corresponding muscle in age-matched normal subjects. Manual methods although important at the time were time consuming and liable to variable sources of error due to the subjective measurement of the MUAP parameters of interest. The advent of computer technology in the last twenty years allowed scientists to renew their efforts in the analysis of EMG so as to improve the assessment of neuromuscular disorders. Computer aided EMG pro-

1350-4533/99/$ - see front matter  1999 IPEM. Published by Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 0 - 4 5 3 3 ( 9 9 ) 0 0 0 7 2 - 7

406

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

cessing saves time, standardizes the measurements and enables the extraction of features which could not be calculated manually. Computer aided EMG analysis has been carried out both in the time domain and in the frequency domain. Time domain techniques that rely on extracting features for classification directly from the MUAP waveform, such as duration, are somewhat difficult to measure automatically, whereas their manual measurement depends on one’s intuition [4]. On the other hand, signal analysis in the frequency domain reveals the frequency characteristics of the MUAP signals [5]. However, this analysis based on the Fourier Transform (FT) power spectrum estimates, suffers from reduced frequency resolution and spectral leakage effects. The limitations inherent in the classical power spectrum density (PSD) estimation methods can be overcome (or improved) by parametric modeling techniques such as autoregressive (AR) or autoregressive moving average (ARMA). These techniques give better PSD estimators, higher spectral resolution, and avoid spectral leakage effects compared to nonparametric or classical spectral estimation techniques [6]. In clinical EMG, AR MUAP analysis was suggested as a diagnostic tool by Coatrieux [7], Maranzana et al. [8] and Basmajian et al. [9] but has not as yet been investigated in detail. Furthermore, cepstral analysis that has been applied to speech recognition for a long time and more recently to surface EMG movement patterns [11], has not yet been applied in MUAP analysis. The objective of this work was to investigate the usefulness of AR and cepstral analysis in the diagnostic assessment of MUAPs recorded from normal (NOR) subjects as well as from subjects suffering with motor neuron disease (MND), and myopathy (MYO). The following MUAP feature sets were studied: (i) time domain measures, (ii) AR spectral measures, (iii) AR coefficients, and (iv) cepstral coefficients. The discriminative power of various selected features from these feature sets was evaluated using the univariate, and the multivariate feature selection analysis methods. Furthermore, classification performance of the four feature sets was investigated using artificial neural networks. The paper is organized as follows. Section 2 describes the material. Section 3 describes the EMG recording method, and time domain analysis. Section 4 covers the AR and cepstral MUAP modeling. Section 5 describes the methodology followed for MUAP feature analysis, and Section 6 describes neural network classification. Section 7 covers the results and discussion, and Section 8 the concluding remarks.

2. Material In this study, a total of 800 MUAPs were collected from the biceps brachii muscle from 12 NOR, 15 MND,

and 13 MYO subjects. Diagnostic criteria for the subjects selected were based on clinical opinion, biochemical data and muscle biopsy. Only subjects with no history or signs of neuromuscular disorders were included in the normal group. Furthermore, the biceps brachii muscle was examined because it is a proximal muscle of the shoulder girdle that is usually affected at an early stage in both MND and MYO. Also, its easy accessibility has made it attractive to study and widely reported in the literature. The time domain, frequency domain, autoregressive, and cepstral findings for the material investigated are given in Section 7.

3. Data collection and time domain analysis The EMG signal was acquired from the biceps brachii muscle using a concentric needle electrode. The EMG signal was freely triggered, bandpass filtered at 3–10 kHz, and sampled at 20 kHz for 5 s with 12 bits resolution. The signal was then lowpass filtered at 8 kHz and downsampled by two. An automatic parametric pattern recognition (PPR) algorithm [12,13] was used to identify MUAPs recorded from the same motor unit. MUAPs belonging to the same group were aligned for averaging at the maximum positive peak. This method yielded similar results except in the case of MUAPs with high shape variability. In this case, MUAP alignment was based on the maximum cross correlation coefficient of MUAPs shifted around the main positive peak. MUAPs were then averaged to obtain the averaged MUAP envelope, which is referred to as MUAP hereafter. A MUAP waveform consisted of 256 points (25.6 ms) adjusted to have the maximum peak in the region of 100 points (10 ms). The following MUAP time domain features were measured automatically as shown in Fig. 1. 1. Duration (Dur): MUAP beginning and ending were identified by sliding a measuring window of length 3 ms and width ±10 µV. 2. Spike duration (SpDur): measured from the first to the last positive peak. 3. Amplitude (Amp): difference between the minimum positive peak and the maximum negative peaks. 4. Area: sum of the rectified MUAP integrated over the duration measure. 5. Spike area (SpArea): sum of the rectified MUAP integrated over the spike duration. 6. Phases (Ph): number of baseline crossings that exceed 25 µV plus one. 7. Turns (T): number of positive and negative peaks separated from the proceeding and following peak by 25 µV. Twenty MUAPs were recorded from each subject.

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

Fig. 1.

407

MUAP time domain parameters.

4. Autoregressive and cepstral modeling For each MUAP waveform, the mean was computed and subtracted from the MUAP signal so as to avoid spectral bias or distortion [6]. Therefore, each subject was represented with a vector of size 20×256 samples.

time series of residual values (prediction errors), and N is equal to 256. The model can be interpreted as a linear system with e(n) as its input (white noise) and x(n) as its output (MUAP time series). The transfer function H(z) for the AR process is:

4.1. Stationarity test

H(z)⫽

AR modeling requires that the signal fitted be stationary over the given interval. The EMG due to its random nature, is considered to be nonstationary. However, Makhoul [14] has shown that an AR model can be applied on a nonstationary process provided that the process can be classified as locally stationary or wide sense stationary. Two tests were applied to investigate MUAP data stationarity, the run test and the reverse arrangement test as given by Bendatt and Piersol [15]. These tests were applied on the mean and standard deviation values of the MUAP signal analysed in 16 segments of 16 samples length. 4.2. AR modeling and model order selection The autoregressive model of the current sample of the signal x(n) is described as a linear combination of previous samples plus an error term e(n) which is independent of past samples:

冘

1

X(z) ⫽ E(z)

冘

(2)

p

1+

akz−k

k⫹1

The spectrum of the model can be estimated from Eq. (2) as: 1

PAR(w)⫽

冘

(3)

p

|1+

jwk 2

ak e |

k⫽1

It was assumed that the spectrum of e(n) satisfies |E(w)| i.e. for the appropriate model order p, it approaches a white noise sequence. The AR coefficients, ak, were calculated using the modified covariance algorithm as described by Marple [6]. The Akaike information criterion (AIC) [16], was used for estimating the optimum AR model order. The model order which minimizes the criterion function was selected. The formula for this criterion is as follows:

p

x(n)⫽⫺

akx(n⫺k)⫹e(n) n⫽0,1,...,N⫺1

(1)

k⫽1

where x(n) are the samples of the modeled signal, ak are the AR coefficients, p is the model’s order, e(n) is the

AIC(p)⫽NLn(ρˆ p)⫹2p

(4)

where rˆ p is the estimated linear prediction error variance for the model with order p. After a model order was selected, the model fitness was checked using the

408

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

residuals e(n) periodogram, the normalized residual energy periodogram, the normalized cumulative residuals periodogram, and the Portmanteau test statistic as described by Box and Jenkins [17]. 4.3. Spectral estimation

1. Bandwidth (BW) is the difference of frequencies at the upper (F2) and lower (F1)—3 dB points of the power spectrum and is given as: BW⫽F2⫺F1

(5)

2. Quality factor (Q) is the ratio of the dominant peak frequency F0 divided by BW and is expressed as: F0 BW

(6)

3. Moments of order 0, 1 and 2: A moment Mj of order j is defined as given by Lindstrom and Petersen [5]:

冘

(7)

4. Moments M0, M1 and M2 were computed for j=0,1 and 2. 5. Centre frequency (CF) or mean power frequency is defined as: N⫺1

f(n)PAR(f(n))

CF⫽

⫽

n⫽0

N⫺1

PAR(f(n))

M1 M0

(8)

n⫽0

6. Median frequency (FMED) is the frequency at which the power spectrum is divided into two regions with equal power defined as:

冘

FMED⫺1

n⫽0

PAR(f(n))⫽

cn⫽⫺an⫺

冘

(1⫺k/n)akcn−k for 1⬍nⱕp and

(10)

k⫽1

冘

n⫺1

cn⫽⫺

(1⫺k/n)akcn−k for n⬎p

k⫽1

where an and cn denote the nth AR and cepstral coefficients respectively. 5. Feature analysis Several features have been proposed in the literature by various researchers for MUAP EMG signature discrimination. Although the classification ability of each feature or of the different feature sets has often been reported, quantitative analysis for evaluating their effectiveness has not been studied thoroughly yet. In this study two feature selection strategies were implemented using the TOOLDIAG software package [18] : (i) the univariate, and (ii) the multivariate covariance methods. 5.1. Univariate analysis

N⫺1

2 Mj ⫽ f(n)j PAR(f(n)) (2p)j+1 n⫽0

冘 冘

c1⫽⫺a1 n⫺1

The AR spectrum was estimated as given by Eq. (3), with the AR coefficients zero padded from p+1 to 4N (1024 samples). The AR power spectrum estimate of each MUAP was normalized with its maximum power value, and the following frequency and power measures were derived:

Q⫽

nal. Alternatively, the cepstral coefficients can be derived directly from the AR coefficients using the formulae [10]:

冘

The univariate analysis (Univar) method ignores the interdependency among the features and analyses them one at a time. A class separability measure ⌬ between two classes is first calculated:

冉 冊

⌬(i,j)⫽1⫺

si+sj m ¯ i−m ¯j

2

(11)

where mi and mj are the means, and si and sj are the standard deviations of classes i and j respectively. This distance measure indicates how well a certain feature can distinguish between two classes. The optimum value of the above separability measure is one and implies that the two classes are well separated. Negative values indicate class overlapping. The Univar criterion is finally obtained by averaging all mutual distances between all classes and is expressed as:

冘冘

k⫺1

N⫺1

PAR(f(n))

(9)

FMED

4.4. Cepstral analysis The cepstrum of a signal is defined as the inverse Fourier transform of the log power spectrum of the sig-

k

2 Criterion (Univar)⫽ ⌬(i,j) k(k−1) i⫽1 j⫽i⫹1

(12)

where k is the number of classes. Once the criterion is calculated, features are ranked according to this criterion. Finally, the best features are selected. This feature selection criterion is limited that it cannot detect highly correlated features. This suggests that correlation analysis should be carried out first before applying the univariate analysis method.

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

409

5.2. Multiple covariance analysis

7. Results and discussion

The multiple covariance (M-Covar) analysis takes into account the dependence among the features. The MCovar algorithm chooses one arbitrary feature from the whole group of features as the first control variable, and for each of the unselected features it calculates the univariate F-ratio. It then orders all features following the M-Covar criterion. The feature with the highest criterion value is passed to the controls and substitutes the first which was initially selected arbitrarily. The above process is repeated for all unselected features. At each iteration the selected feature with the highest F-ratio joins the already selected features until the desired number of selected features is reached. Details of the multiple covariance analysis method can be found in Michalski and Larson [19]. The computational time required for the M-Covar algorithm is higher than that of the Univar algorithm. The advantage of the method is that it represents highly correlated features by only one discriminative feature. The rest of the features are classified as features with low univariate F-ratio.

A total of 800 MUAPs were analysed, recorded from 12 NOR, 15 MND and 13 MYO subjects. Time domain, autoregressive and cepstral feature sets are presented and discussed in this section followed by statistical analysis and neural network classification.

6. Neural network classification Artificial neural network (ANN) models were developed for the classification of MUAP features into three groups: NOR, MND, and MYO. The motivation behind the use of ANN lies in their capacity for making no assumptions about the underlying probability density functions of the input data, finding near-optimum solutions from incomplete data sets, and the fact that learning is accomplished through training [20]. For each subject, the average vector of 20 MUAPs per subject for each feature set was computed and used as input to the classifiers. The following four different feature sets extracted from the MUAP waveforms as described in Sections 3 and 4 were used: 1. Time domain parameters (n=7): Dur, SpDur, Amp, Area, SpArea, Ph, T 2. Frequency domain parameters (n=6): M0, M1, M2, FMED, BW, Q 3. Autoregressive coefficients (n=12): a1 to a12 4. Cepstral coefficients (n=12): c1 to c12 where n is the number of features. Three different neural network classifiers were investigated, the back-propagation (BP) [21], the radial-basis function network (RBF) [22], and the self-organizing feature map (SOFM) [23]. These algorithms were implemented as given in the MATLAB neural network toolbox [24].

7.1. Time domain analysis Time domain parameters are the most widely used parameters in clinical neurophysiology for the interpretation of EMG findings [4]. The following time domain parameters were computed from the MUAP waveforms as presented in Section 3, and illustrated in Fig. 1: duration, spike duration, amplitude, area, spike area, number of phases, and number of turns. The statistics for the NOR, MND, and MYO groups are tabulated in Table 1. For the NOR group, the mean and standard deviation values computed for duration, amplitude, and number of phases were 9.60±2.75 ms, 0.376±0.306 mV, and 2.6±0.8, respectively. For the MND group, the mean and standard deviation values for duration, amplitude, and number of phases were 13.42±3.86 ms, 0.614±0.426 mV, and 4.0±1.8, respectively. Typically the mean durations of the MUAPs are longer and the mean amplitudes, and the number of phases are increased in patients suffering from motor neuron disease. In this disease there is an increase in the number or density of the fibers in the motor unit due to denervation of the motor units. There is also an increase in the temporal dispersion of the activity picked up by the recording electrode. This effect is the result of slowed conduction along the terminal branches of individual nerve fibres, increase in the end-plate zone, or both. For the MYO group, the mean and standard deviation values for duration, amplitude, and number of phases were 7.15±2.34 ms, 0.314±0.250 mV, and 2.7±1.0, respectively. MUAPs with short duration and reduced amplitude are typical findings in patients suffering from myopathy. These findings are attributed to fiber loss within the motor unit, with the degree of reduction of these parameters reflecting the amount of fiber loss. 7.2. Autoregressive and cepstral analyses 7.2.1. Stationarity Autoregressive modeling can be applied on a nonstationary process provided that the process is locally stationary or wide sense stationary. Two tests were used to investigate whether the MUAP data satisfy the stationarity criterion in the wide sense, the run test and the reverse arrangements test, applied on the mean and standard deviation values, of the MUAP signals of 256 points, of 25.6 ms epoch. Results showed that for the NOR, MND, and MYO groups, 85% of the MUAP sig-

410

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

Table 1 Statistics of time domain MUAP parameters for the NOR, MND, and MYO groups Class

NOR MND MYO

Duration (ms) Mean

SD

Spike duration (ms) Mean SD

9.60 13.42 7.15

2.75 3.86 2.34

5.40 6.86 4.21

2.46 4.09 1.83

Amplitude (mV)

Area(m Vms)

Mean

SD

Mean

0.376 0.614 0.314

0.306 0.426 0.250

0.370 0.832 0.234

nals in average passed the stationarity criterion for both tests at a significance level of a=0.05. Wide sense stationarity of EMG signals was also shown for surface EMG of intervals shorter than 1.0 s using the run test by Inbar and Noujaim [25], Hefftner et al. [26], and Popivanov and Todorov [27]. 7.2.2. AR modeling and model order selection The modified covariance method was used to compute the AR coefficients, and the model order was selected using the AIC criterion. The Fast Fourier Transform (FFT), and the AR spectra for two MUAPs are plotted in Fig. 2. AR spectra are shown in dotted lines, and FFT spectra in solid lines. Generally, the AR spectra either followed closely the FFT spectra, or formed the spectral envelope in the case of spiky waveforms. The mean and standard deviation values of model order p for the AIC criterion, and for the AR coefficients

Fig. 2.

Phases

SD

Spike area (mVms) Mean SD

Mean

SD

Mean

SD

0.198 0.590 0.193

0.232 0.520 0.163

2.6 4.0 2.7

0.8 1.8 1.0

3.0 4.7 3.1

2.0 2.5 1.2

0.118 0.369 0.132

Turns

a1 to a12 are given in Table 2. The mean and standard deviation values of model order p were 6.08±3.01 for the NOR group, 7.67±4.19 for the MND group, and 5.98±3.28 for the MYO group. Results showed that higher model orders were needed for the MND subjects compared to the NOR subjects due to the higher complexity in shape. On the other hand, slightly lower model orders were obtained in the MYO group compared to the NOR group due to the rather simpler form of MUAPs in the MYO group. The above observations were also confirmed by histograms of model order p for the NOR, MND, and MYO groups as illustrated in Fig. 3. Furthermore, these results clearly indicated that twelve AR coefficients were adequate for MUAP modeling. Once a model order has been identified and the parameters estimated, diagnostic checks must be applied to the model to determine its adequacy. These checks are customarily applied to the residuals of the model, i.e the

Sample MUAPs, and their corresponding FFT (———) and AR spectra (· · ·). MUAP model order p is also given.

0.081842 0.15883 0.097425 0.093753 ⫺0.002984 0.188069 ⫺0.000086 0.105679 0.010371 0.136583 0.220224 0.117944 ⫺0.000872 ⫺0.005655 0.013743 0.017705 ⫺0.045432 ⫺0.020974 0.192322 0.261214 0.194694 0.113423 0.293977 ⫺0.091136 0.114510 0.345164 ⫺0.096369 0.105285 0.271477 ⫺0.060686 NOR MND MYO

0.234101 0.296878 0.219653

0.018550 0.079257 0.046554

SD a9 Mean a7 Mean SD a6 Mean Class

SD

a8 Mean

SD

0.182741 0.237154 0.158845

SD a10 Mean

0.000612 0.019720 ⫺0.008572

a11 Mean

SD

0.371162 0.417175 0.365019 ⫺0.130027 ⫺0.165833 ⫺0.192223 0.486064 0.522950 0.463620 0.210669 0.339362 0.330669 0.641230 0.671467 0.571509 ⫺0.417521 ⫺0.635927 ⫺0.631906 0.691752 0.718538 0.589850 0.420822 0.439856 0.415501 3.01 4.19 3.28 6.08 7.67 5.98 NOR MND MYO

⫺1.668010 ⫺1.723160 ⫺1.579977

1.092608 1.214877 1.204518

SD a3 Mean SD a2 Mean SD a1 Mean SD p Mean Class

Table 2 Statistics of model order p, and AR coefficients a1 to a12 for NOR, MND, and MYO groups

a4 Mean

SD

a5 Mean

SD

a12 Mean

SD

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

411

error between the actual signal and the signal predicted by the model. The order of the model must so be chosen to guarantee that the residual is nearly white noise. AR model validity was initially checked by a visual inspection of the residuals periodograms, and the normalized residuals periodogram of each MUAP. It was observed that the residuals periodogram generally resembled the spectrum of a white noise sequence with slight variations. This indicated that the AR model represented adequately the MUAP data. The model fitness was also checked by the normalized cumulative residuals periodograms and the portmanteau test statistic Q [17]. Results showed that for all MUAPs, the estimated model orders produced white noise residuals, and the Q statistic hypothesis was satisfied, thus accepting the hypothesis of model adequacy. The modified covariance method was implemented in this study for computing the AR coefficients. This method was preferred over the Yule–Walker, the Burg and the covariance methods, as it has been documented that it provides higher frequency resolution, and avoids frequency estimation biases and spectral line splitting problems that are associated with one or more of the above power spectral estimators [6]. The AIC criterion for needle EMG AR model order selection was also used by Coatrieux [7], and Maranzana et al. [8]. Coatrieux [7] suggested model orders in the range of 11–13 depending on muscle force, whereas Maranzana et al. [8] proposed model orders up to nine. Basmajian et al. [9] suggested an AR model of six in a needle MUAP study. Various model orders in the range of 3–20 for surface EMG modeling were reported in the literature using different selection criteria [25,28–31]. 7.2.3. AR spectral estimation The following frequency and power measures were determined from the AR model power spectrum of each MUAP: CF, FMED, F0, F1, F2, BW, Q, M0, M1, and M2. The mean and standard deviation values for these measures are given in Table 3. Mean and standard deviation values of CF for the NOR, MND and MYO groups were: 522±332, 446±283 and 755±355 Hz. Similarly, FMED for the NOR, MND, and MYO groups were: 413±285, 339±241 and 629±332 Hz. It is clearly shown that in the case of the MYO group, CF and FMED were higher compared to the NOR group, illustrating that the power spectrum in the case of myopathies was shifted towards higher frequencies. This was also reported in a number of studies by Richardson [32], Watlon [33], Fex and Krakau [34], Gersten et al. [35], Kopec and Hausmanowa-Petrusewicz [36], Carlsson et al. [37], Sandstedt [38], and Lindstrom and Petersen [5]. On the other hand, CF and FMED in the case of the MND group were lower compared to the NOR group, illustrating that the power spectrum was shifted in lower frequencies. This finding was also reported by Kaiser and Petersen [39],

Fig. 3. Histograms of median frequency, quality factor, and model order for the NOR, MND, and MYO groups. Median, semi interquartile range (SIQR) and percentage of polyphasic MUAPs are given for each group. Crossed bars show polyphasic MUAPs.

412 C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

413

Table 3 Statistics of AR spectral measures for the NOR, MND, and MYO groups Class Centre frequency [CF] (Hz) Mean SD NOR 522 MND 446 MYO 755

322 283 355

Median frequency [FMED] (Hz) Mean SD

Peak frequency [F0] (Hz)

Lower-3 db Upper-3 db frequency frequency [F1] (Hz) [F2] (Hz)

Bandwidth [BW] (Hz)

Quality factor [Q]

Moment M0 Moment M1 Moment M2 (µV2) (µV2/s×103) (µV2/s2×106)

Mean SD

Mean SD

Mean SD

Mean SD

Mean SD

Mean SD

413 339 629

215 213 411

35 61 104

560 442 871

525 381 767

0.47 0.71 0.63

20.71 15.02 15.34 23.06 24.52 51.95 15.30 11.80 9.54 17.59 14.12 39.70 29.70 16.88 27.68 27.18 45.41 64.95

285 241 344

278 240 389

104 124 200

Gersten et al. [35], Kopec and Hausmanowa-Petrusewicz [36], Larsson [40], Shigiya et al. [41], Sandstedt et al. [38], and Lindstrom and Petersen [5]. Furthermore, BW in the MYO group is higher than in the NOR group, with the reverse being valid for the MND group. The histograms of the AR spectral measures FMED and Q are given in Fig. 3. Spectral moments of order 0, 1 and 2 are also tabulated in Table 3. Spectral moment M0 is actually the average power of the spectrum. Spectral moment M0 was found to be higher in the MYO group than in the NOR group. M0 in the MYO group was 29.70 µV2, and in the NOR group 20.71 µV2. Conversely, M0 in the MND group was smaller than M0 in the NOR group. Actual mean values of M0 in the MND and NOR groups were 15.30 and 20.71 µV2, respectively. Spectral moments M1 and M2 represent the area of the power spectrum after being weighted by frequency raised to the power of one and two respectively. Thus, since all the frequency parameters in the MYO group were found to be of greater value and of higher power level, it follows that M1 and M2 in the MYO group would be greater than M1 and M2 in the NOR group. This is true since M1 and M2 were 27.68 µV2/s×103 and 45.41 µV2/s2×106 in the MYO group respectively and 15.34 µV2/s×103 and 24.52 µV2/s2×106 in the NOR group respectively. A similar justification is valid for M1 and M2 of the MND group when compared to the NOR group. M1 and M2 in the MND group were 9.54 µV2/s×103 and 14.12 µV2/s2×106 respectively, whereas M1 and M2 in the NOR group were 15.34 µV2/s×103 and 24.52 µV2/s2×106 respectively. Histograms of moments M0, M1, and M2 are given in Fig. 4. It is noted that although these measures were proposed by Lindstrom and Petersen [5], they have not been further investigated. 7.2.4. Cepstral analysis Twelve cepstral coefficients, c1 to c12 were determined directly from the estimated AR coefficients a1 to a12. Their group statistics for the NOR, MND, and MYO groups are given in Table 4. The classification performance of the cepstral coefficients c1 to c12 was also examined by neural network analysis as given in the following section. The motivation of using cepstral coefficients as

428 361 549

432 349 539

0.57 0.98 0.78

Mean SD

Mean SD

diagnostic features in this study was due to their documented higher identification accuracy in speech recognition compared to the AR coefficients [10]. Furthermore, in a recent study for classifying patterns of movement via surface EMG signals, it was reported that cepstral coefficients showed significantly better separability compared to the AR coefficients [11]. 7.3. Feature analysis One way analysis of variance (ANOVA) of the following selected time domain and AR parameters was carried out as tabulated in Table 5: Dur, SpDur, Amp, T, M0, M1, M2, CF, FMED, F0, BW, Q, and AR model order p. t-Test analysis was also carried out between the NOR–MND, NOR–MYO, and MND–MYO groups. In the NOR–MND group, there was significant difference for all parameters except in peak power frequency, F0. In the NOR–MYO group, there was no significant difference in the parameters Amp, and p. Furthermore, in the MND–MYO group there was significant difference in all parameters except from quality factor, Q. Assumptions for ANOVA analysis to be valid are given in Table 5. Assumptions number 3 and 4 require that samples are normally distributed within each group, and that homogeneity exists between groups. These assumptions may be validated for the MUAPs under investigation. However, the assumption of normality for each parameter in each group is violated as demonstrated by histograms given in Figs. 3 and 4. Correlation analysis was also carried out to establish the interrelationship between the following features: Dur, SpDur, Amp, T, M0, M1, M2, CF, FMED, F0, BW, Q, and p. Correlation analysis results are tabulated in Table 6. Time parameters were not strongly correlated with the rest of the features. There was a high correlation between spectral moments M0, M1, M2, and frequency parameters CF, FMED, and BW. The remaining parameters were not highly correlated between each other or with the aforementioned parameters. The discriminative power of a small number of features was investigated using the univariate and the multiple covariance analysis methods (see Table 7). The following features were analysed: Dur, SpDur, Amp, CF,

Fig. 4. Histograms of moments M0, M1 and M2 for the NOR, MND, and MYO groups. Median, semi interquartile range (SIQR) and percentage of polyphasic MUAPs are given for each group. Crossed bars show polyphasic MUAPs.

414 C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

415

Table 4 Statistics of cepstral coefficients c1 to c12 for the NOR, MND, and MYO groups Class

c1 Mean

SD

NOR MND MYO

1.668010 1.723160 1.579977

0.420822 0.386697 0.439856 0.366177 0.415501 0.129635

Class

c7 Mean

SD

0.040310 0.057157 ⫺0.003685

0.137804 0.037209 0.177345 0.022531 0.119866 ⫺0.012868

NOR MND MYO

c2 Mean

c8 Mean

SD

c3 Mean

0.431402 0.200811 0.421948 0.30473 0.384323 0.080628

SD

c9 Mean

0.144269 0.020981 0.199724 0.007635 0.109165 ⫺0.005609

SD

c4 Mean

SD

c5 Mean

c6 Mean

SD

SD

0.215739 0.123910 0.151817 0.101349 0.118718 0.049142 0.116364 0.243998 0.89913 0.180572 0.131919 0.151253 0.074525 0.146572 0.180367 0.014515 0.128453 0.021473 0.132818 0.006215 0.124227

SD

c10 Mean

SD

c11 Mean

c12 Mean

SD

SD

0.179839 0.023168 0.212216 0.024834 0.291446 0.040019 0.366013 0.225184 0.025021 0.261658 0.017820 0.308422 0.027882 0.386145 0.109275 0.005091 0.134606 0.017915 0.138857 0.019449 0.131812

Table 5 One way analysis of variance (ANOVA) of various selected time and AR spectral parametersa

Duration Spike Duration Amplitude Turns Moment M0 Moment M1 Moment M2 Centre frequency Median frequency Peak frequency Bandwidth Quality factor Model Order p a b

H0: mnNOR=mnMND=mnMYO Significant difference Pair comparisons (t-tests) NOR–MYOb NOR–MNDb

MND–MYOb

Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes

Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes No Yes

Yes Yes Yes Yes Yes Yes Yes Yes Yes No Yes Yes Yes

Yes Yes No Yes Yes Yes Yes Yes Yes Yes Yes Yes No

Assumptions: (1) Random samples; (2) independent groups; (3) normal distributions [N (mn, s2)]; (4) homogeneity (s2NOR=s2MND=s2MYO). Significant difference at level a=0.05.

FMED, M0, Q, and p. The number of turns T, though not highly correlated with any of the above features was removed so as to have a small number of time parameters in the final feature set. Additionally, M2 though is highly correlated with FMED, was included in the feature selection analysis aiming at diversity. Results given in Table 7 show that the univariate analysis ordered features in the following series: Dur, FMED, SpDur, M2, Amp, Q, and p; whereas the order of features for the modified covariance analysis was: FMED, Dur, p, Amp, M2, Q, and SpDur. It is clearly shown that duration, Dur, is the best discriminative feature among time domain parameters, and FMED is the best feature among AR spectral measures. Model order p was ranked among the poorest classification features in the univariate method, whereas in the modified covariance method, it was classified third. Pfeiffer and Kunze [42] investigated the diagnostic usefulness of parameters Dur, SpDur, Amp, T, Ph, and

CF for the classification of NOR, neuropathic, and MYO subjects. They showed that CF and SpDur were the most useful parameters. In this study, FMED and Dur were found to be the best discriminative features. Comparing the findings of the two studies, it can be inferred that there is agreement regarding the frequency features, as CF and FMED are highly correlated, whereas for the time features there is no agreement, as Dur and SpDur are moderately correlated. 7.4. Neural network classification For training the ANN classifiers, a total of 480 MUAPs obtained from 24 subjects, 8 NOR, 8 MYO and 8 MND, were used, whereas for evaluation, a total of 320 MUAPs, obtained from 16 subjects, 4 NOR, 5 MYO and 7 MND were used. The system was trained and evaluated using five different bootstrap sets where in each set 24 different subjects were selected at random

Duration Spike duration Amplitude Turns Moment M0 Moment M1 Moment M2 Centre frequency Median frequency Peak frequency Bandwidth Quality factor Model order p

1.000 0.172 0.352 ⫺0.076 ⫺0.067 ⫺0.054 ⫺0.051

⫺0.087

⫺0.119

⫺0.092 ⫺0.002 0.103

⫺0.229

⫺0.165

⫺0.210 0.032 0.157

Spike duration (SpDur)

1.000 0.474 0.402 0.308 ⫺0.228 ⫺0.182 ⫺0.142 ⫺0.224

Duration (Dur)

0.110 0.026 0.091

0.122

0.172

1.000 0.200 0.124 0.163 0.170 0.182

Ampiture (Amp)

⫺0.020 0.084 0.146

0.028

0.090

1.000 0.020 0.026 0.018 0.128

Turns (T)

0.956 ⫺0.140 ⫺0.339

0.607

0.905

1.000 0.952 0.877 0.888

0.901 ⫺0.063 ⫺0.264

0.648

0.916

1.000 0.978 0.899

0.827 ⫺0.054 ⫺0.241

0.601

0.853

1.000 0.844

0.760 ⫺0.015 ⫺0.183

0.577

0.959

1.000

Moment M0 Moment M1 Moment M2 Centre frequency (CF)

0.807 0.086 ⫺0.157

0.692

1.000

Median frequency (FMED)

Table 6 Correlation matrix between various selected time domain, and AR spectral parameters for the NOR, MND, and MYO groups

0.633 0.431 ⫺0.010

1.000

Peak frequency (F0)

1.000 ⫺0.157 ⫺0.366

Bandwidth (BW)

1.000 0.439

Quality factor (Q)

1.000

Model order (p)

416 C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

417

Table 7 Feature selection using the univariate and the multiple covariance analysis methods Feature no.

Feature

Univar criterion

Feature

M-Covar criterion

1 2 3 4 5 6 7

Duration Median frequency Spike duration Moment M2 Amplitude Quality factor Model order p

0.325 ⫺1.749 ⫺3.284 ⫺3.520 ⫺8.363 ⫺9.497 ⫺114.172

Median frequency Duration Model order p Amplitude Moment M2 Quality factor Spike duration

25.706 12.381 5.424 2.188 1.669 1.198 0.096

for training, and 16 different subjects for evaluation as described above. For each subject, the average vector of 20 MUAPs per subject for each feature set was computed and used as input to the classifiers. Table 8 tabulates the mean and standard deviation of the percentage of correct classifications, i.e. diagnostic yield (DY) for the five bootstrap sets for each classifier for the evaluation set. As shown in the average diagnostic yield column of Table 8, the highest diagnostic yield was obtained for the time domain parameters, being 78.3%, followed by the cepstral coefficients, the frequency domain parameters, and the autoregressive coefficients with diagnostic yields of 63.7, 62.5 and 51.2%, respectively. It is important to note that the time domain parameters which are empirically defined and are the only parameters which are extensively used in the clinical neurophysiology laboratory for the assessment of neuromuscular disorders gave the highest diagnostic yield. A similar diagnostic yield in the region of 78–80% for the time domain parameters was also obtained using neural networks [12,13], and genetics-based machine learning [43] based on the same material as used in this study. Furthermore, the better performance of the cepstral coefficients over the autoregressive coefficients for surface EMG was also reported by Kang et al. [11]. In a very recent paper, we have shown that MUAP

time domain measures, AR spectral measures, AR coefficients, cepstral coefficients, as well as wavelet coefficients can be combined into a modular neural network decision support system in order to improve the classification performance of individual features or classifiers [44]. This method was based on the same material as used in this study. More specifically, a decision making system was developed where BP, RBF, and SOFM neural network classifiers were trained as documented in this study (as shown in Table 8). In addition, the BP, RBF, and SOFM classifiers were trained with wavelet coefficients computed using the following techniques [45] Dubechies with 4 and 20 coefficients, Chui, and Battle– Lemarie. The corresponding average diagnostic yields for these algorithms were 66.2, 59.6, 63.3 and 65.8%, respectively. (The motivation behind the use of wavelet analysis is that it provides a linear time-scale frequency representation of localized statistical measures (the scalogram) for nonstationary signal analysis like MUAP’s. See Appendix A on wavelet analysis). The 24 output results (of the eight feature sets applied to the three classifiers) were combined using majority voting, i.e a subject was assigned to the class with the maximum number of votes. The diagnostic yield achieved with the combination of the 24 classification results was 82.5%. This diagnostic yield was higher than the diagnostic yield of the best feature set from the best classifier,

Table 8 Mean and standard deviation of the diagnostic yield, DY%, for the evaluation set of the neural network models after bootstrapping the available data for five different sets of subjects Feature set

na

Back-propagation BP (DY%)

Radial basis function network RBF (DY%)

Self-organizing Average DY% feature map SOFM (DY%)

(i) Time domain (duration, spike duration, amplitude, area, spike area, phases, and turns) (ii) Frequency domain (moments M0, M1, and M2, median frequency, bandwidth, and quality factor) (iii) Autoregressive (a1 to a12) (iv) Cepstral (c1 to c12)

7

81.2±11.8

72.5±8.5

81.2±5.5

78.3±8.6

6

63.7±2.5

60.0±6.3

63.7±10.1

62.5±6.3

12 12

56.2±7.9 62.5±15.8

50.0±6.8 60.0±5.0

47.5±10.1 68.7±3.9

51.2±8.3 63.7±8.2

a

Feature set vector size. For each subject, the average vector of 20 MUAPs per subject for each feature was computed.

418

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

which was 81.2% for the time domain parameters and the SOFM and BP classifiers as given in Table 8. On the other hand, it should be emphasized that although, the overall diagnostic yield of the system proposed in [44] is close to the best feature set and the best classifier diagnostic yield of 81.2%, the overall diagnostic yield score is more reliable since its computation is based on different representations of the MUAP data and different generalizations of the neural network classifiers. Furthermore, it is noted that selecting the best classifier or the best feature set is not necessarily the ideal choice, since potentially valuable information contained in the less successful feature sets or classifiers may be lost. The combination of the results of the different features and the different classifiers increases the probability that the errors of the individual features or classifiers may be compensated by the results of the rest [46]. 8. Concluding remarks Autoregressive and cepstral analyses were investigated in this study for the assessment of MUAPs recorded from NOR, MND, and MYO subjects. Results showed that the cepstral coefficients provided a better diagnostic yield than the AR spectral measures, and the AR coefficients. However, the diagnostic yield obtained with the time domain parameters was better than the cepstral coefficients. On the other hand, the diagnostic yield of the cepstral coefficients and the AR spectral measures was similar to the diagnostic yield obtained with wavelet analysis. Thus, AR and cepstral analysis, combined with time domain and wavelet analyses, through a modular neural network decision support system can provide useful information for the assessment of neuromuscular disorders. Furthermore, it was shown that the MUAP duration and median frequency carry similar discriminative power which is higher than the rest of the time and AR spectral measures investigated. Future work is needed to investigate the interpretation of AR and cepstral coefficient changes in relation to muscle pathophysiology: MUAP propagation velocity, velocity dispersion, size of innervation zone, number of fibers in the motor unit, and electrode recording site. Acknowledgements The authors would like to express their sincere thanks to the Cyprus Institute of Neurology and Genetics for their support in this work. Appendix A. Wavelet analysis The following four different wavelets were investigated for the analysis of the same MUAP dataset that

was used in this study [45]: Daubechies with 4 (DAU4) and 20 (DAU20) coefficients, Chui (CH) and Battle– Lemarie (BL). The main characteristics of these wavelets which motivated their application on MUAP analysis are outlined below: 1. The orthogonal wavelet transform (WT) decomposes a signal into a set of orthogonal basis functions. Unlike the short time fourier transform (STFT) which uses orthogonal sinusoids over the window, while the basis functions from overlapping windows are not orthogonal, the WT basis functions are always orthogonal. This means that every WT coefficient represents an entirely different signal component. 2. Assuming that high frequencies change rapidly, while low frequencies change slowly, the WT uses long duration windows for capturing the low frequency components and short duration windows for capturing the high frequency components. This means that for MUAP signals we look at the high frequency coefficients for locating abrupt signal changes, while we look at the low frequency coefficients for capturing the average behaviour of the signal over long durations. In contrast, the STFT uses a fixed window and does not provide a varying time-frequency resolution. (For the STFT, the use of a short window results in low spectral resolution, whereas the use of a longer window improves the frequency resolution but results in a loss of time-resolution). 3. For the case of the DAU4, the WT coefficients at the highest frequency scales provide high time-resolution of only four signal samples. This allows the DAU4 wavelet to effectively track the MUAP main spike transient signal at a time resolution that the STFT simply cannot match. 4. The energy in the MUAP signal is distributed among a small number of well localized in time WT coefficients. In addition to the time-frequency features, the WT supports multiresolution analysis of the MUAP signals. Briefly, at each stage (scale) the WT provides two sets of coefficients: the scaling function coefficients, and the wavelet function coefficients. Each successive stage replaces the scaling function coefficients by another two sets of scaling function and wavelet coefficients. In terms of analysing MUAP signals, we compare the difference between the two sets of scaling function coefficients, coming from successive scales, to recognize the signal features that the set of wavelet coefficients represents.

References [1] Buchthal F, Pinelli P, Rosenfalck P. Action potential parameters in normal human muscle and their physiological determinations. Acta Physiol Scand 1954;32:219–29.

C.S. Pattichis, A.G. Elia / Medical Engineering & Physics 21 (1999) 405–419

[2] Buchthal F, Guld C, Rosenfalck P. Action potential parameters in normal human muscle and their dependence on physical variables. Acta Physiol Scand 1954;32:200–15. [3] Buchthal F. An introduction to electromyography. Copenhagen: Gyldendal, 1957. [4] Stalberg E, Andreassen S, Falck B, Lang H, Rosenfalck A, Trojaborg W. Quantitative analysis of individual motor unit potentials: a proposition for standardised terminology and criteria for measurement. J Clin Neurophysiol 1986;3:313–48. [5] Lindstrom L, Petersen I. Power spectrum analysis of EMG signals and its applications. In: Desmedt JE, editor. Computer-Aided Electromyography Prog Clin Neurophysiol, vol. 10. Basel: Karger, 1983:1–51. [6] Marple SL Jr. Digital spectral analysis with applications. Englewood Cliffs: Prentice Hall, 1987. [7] Coatrieux JL. Interference electromyogram processing—Part II: experimental and simulated EMG AR modeling. Electromyogr Clin Neurophysiol 1983;23:481–90. [8] Maranzana MF, Molinari R, Somma-Riva G. The parametrization of the electromyographic signal: an approach based on simulated EMG signals. Elect Clin Neurophysiol 1984;24:47–65. [9] Basmajian JV, Gopal DN, Ghista DN. Electrodiagnostic model for motor unit action potential (MUAP) generation. American J Phys Medic 1985;64:279–94. [10] Atal BS. Automatic recognition of speakers from their voices. IEEE Proc 1976;64:460–75. [11] Kang W, Shiu JR, Chen CK, Lai JS, Tsao HW, Kuo TS. The application of cepstral coefficients and maximum likelihood method in EMG pattern recognition. IEEE Trans Biomed Eng 1997;42(8):777–85. [12] Pattichis CS. Artificial neural networks in clinical electromyography. PhD dissertation, QMW College, University of London, 1992. [13] Pattichis CS, Schizas CN, Middleton LT. Neural network models in EMG diagnosis. IEEE Trans Biomed Eng 1995;42:486–96. [14] Makhoul J. Linear prediction: a tutorial review. Proc IEEE 1975;63:561–80. [15] Bendat JS, Piersol AG. Random data: analysis and measurement procedures. New York: Wiley Interscience, 1986. [16] Akaike H. A new look at the statistical model identification. IEEE Trans Autom Contr 19, 716–23, December 1974. [17] Box GEP, Jenkins GM. Time series analysis: forecasting and control. 2nd ed. San Francisco: Holden-Day, 1976. [18] Rauber TW, Barata MM, Steiger-Garcao AS. A toolbox for analysis and visualization of sensor data in supervision. Proceedings of the International Conference on Fault diagnosis, Toulouse, France, 1993. [19] Michalski RS, Larson JB. Selection of most representative training examples and incremental hypotheses: the underlying methodology and the description of programs AE SEL and AQ11V. Research report UIUCD CS-R 78-867, University of UrbanaChampaign, Illinois, USA, 1978. [20] Haykin S. Neural networks—a comprehensive foundation. Macmillan College Publishing Company, 1994. [21] Rumelhart DE, Hinton GE, Williams RJ. Learning internal representations by error propogation. In: Rumelhart DE, McClelland JL, editors. Parallel distributed processing, vol. 1. Cambridge, MA: MIT Press, 1986:318–62. [22] Broomhead DS, Lowe D. Multivariable functional interpretation and adaptive networks. Complex Syst 1988;2:321–55. [23] Kohonen T. Self-organizing maps. Berlin: Springer, 1995. [24] Demuth H, Beale M. Neural network toolbox user’s guide. Natick, MA: The Maths Work Inc, 1994. [25] Inbar GF, Noujaim AE. On surface EMG spectral characterization and its application to diagnostic classification. IEEE Trans Biomed Eng 1984;31(9):597–604.

419

[26] Hefftner G, Zucchini W, Jaros GG. The electromyogram (EMG) as a control signal for functional neuromuscular stimulation— Part I: autoregressive modeling as a means of EMG signature discrimination. IEEE Trans Biomed Eng 1988;35:597–604. [27] Popivanov D, Todorov A. Statistical procedures for interference EMG power spectra estimation. Med Biol Eng Comput 1986;24:344–50. [28] Graupe D, Cline W. Functional separation of EMG signals via ARMA identification methods for prosthesis control purposes. IEEE Trans Syst Man Cybern 1975;5:252–9. [29] Doerschuk PC, Gustafson SE, Willsky AS. Upper extremity limb function discrimination using EMG signal analysis. IEEE Trans Biomed Eng 1983;30:18. [30] Paiss O, Inbar GF. Autoregressive modeling of surface EMG and its spectrum with applications to fatigue. IEEE Trans Biomed Eng 1987;34:761–70. [31] Triolo RJ, Nash DH, Moskowitz GD. The identification of times series models of lower extremity EMG for the control of prostheses using Box-Jenkins criteria. IEEE Trans Biomed Eng 1988;35:584–93. [32] Richardson AT. The analysis of muscle action potentials in the differential diagnosis of neuromuscular diseases. Arch Phys Med 1951;32:199–206. [33] Walton JN. The electromyograph in myopathy: analysis with the audio frequency spectrometer. J Neurol Neurosurg Psychiat 1952;15:219–26. [34] Fex J, Krakau CET. Some experiences with Walton’s frequency analysis of the electromyogram. J Neurol Neurosurg Psychiat 1957;20:178–84. [35] Gersten JW, Cenkovich FS, Jones GD. Harmonic analysis of normal and abnormal electromyograms. Am J Physic Med 1965;44:235–40. [36] Kopec J, Hausmanowa-Petrusewicz I. Application of harmonic analysis to the electromyogram’s evaluation. Acta Physiol Polonica 1966;17:597–608. [37] Carlsson C, Dencker J, Henriksson KG, Magnusson I, Petersen I, Riman E. Clinical, histological, and electromyographical studies in chronic alcoholics. Acta Neurol Scand 1972;48(51, Suppl.):425–7. [38] Sandstedt P. Quantitative examination in neuromuscular disorders studied by muscle biopsy and electromyography. Linkoping University Medical Dissertations No 121, Linkoping, 1981. [39] Kaiser E, Petersen I. Muscle action potentials studied by frequency analysis and duration measurement. Acta Neurol Scand 1965;41:213–36. [40] Larsson LE. On the relation between the EMG frequency spectrum and the duration of symptoms in lesions of the peripheral motor neuron. Electroenceph Clin Neurophysiol 1975;38:69–78. [41] Shigiya R, Itoh K, Suhara K, Sameshima M, Suzuki H. Spectral analysis of surface electromyogram. Electroenceph Clin Neurophysiol 1973;34:799–800. [42] Pfeiffer G, Kunze K. Frequency analysis and duration of motor unit potentials: Reliability and diagnostic usefulness. Electroenceph Clin Neurophysiol 1993;89:365–74. [43] Pattichis CS, Schizas CN. Genetics based machine learning for assessment of certain neuromuscular disorders. IEEE Trans Neur Net 1996;7:427–39. [44] Christodoulou CI, Pattichis CS, Fincham WF. A modular neural network decision support system in EMG diagnosis. Special Issue on Computational Intelligent Diagnostic Systems in Medicine. J Intell Syst 1998;8(12):99–143. [45] Pattichis CS, Pattichis MS. Time scale analysis of motor unit action potentials. IEEE Trans Biomed Eng 1999;46(11): in press. [46] Tumer K, Ghosh J. Theoretical foundations of linear and order statistics combiners for neural pattern classifiers. Technical report, Dept. of Electrical and Computer Engineering, University of Texas, Austin, 1996.