Robust Classification of Intramuscular EMG Signals to Aid the Diagnosis of Neuromuscular Disorders

Goal: This article presents the design and validation of an accurate automatic diagnostic system to classify intramuscular EMG (iEMG) signals into healthy, myopathy, or neuropathy categories to aid the diagnosis of neuromuscular diseases. Methods: First, an iEMG signal is decimated to produce a set of “disjoint” downsampled signals, which are decomposed by the lifting wavelet transform (LWT). The Higuchi's fractal dimensions (FDs) of LWT coefficients in the subbands are computed. The FDs of LWT subband coefficients are fused with one-dimensional local binary pattern derived from each downsampled signal. Next, a multilayer perceptron neural network (MLPNN) determines the class labels of downsampled signals. Finally, the sequence of class labels is fed to the Boyer-Moore majority vote (BMMV) algorithm, which assigns a class to every iEMG signal. Results: The MLPNN-BMMV classifier was experimented with 250 iEMG signals belonging to three categories. The performance of the classifier was validated in comparison with state-of-the-art approaches. The MLPNN-BMMV has resulted in impressive performance measures (%) using a 10-fold cross-validation—accuracy = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$99.87\pm 0.25$\end{document}, sensitivity (normal) = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$99.97\pm 0.13$\end{document}, sensitivity (myopathy) = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$99.68\pm 0.95$\end{document}, sensitivity (neuropathy) = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$99.76\pm 0.66$\end{document}, specificity (normal) = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$99.72\pm 0.61$\end{document}, specificity (myopathy) = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$99.98\pm 0.10$\end{document}, and specificity (neuropathy) = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$99.96\pm 0.14$\end{document}—surpassing the existing approaches. Conclusions: A future research direction is to validate the classifier performance with diverse iEMG datasets, which would lead to the design of an affordable real-time expert system for neuromuscular disorder diagnosis.


I. INTRODUCTION
A. iEMG Signal Changes in Neuromuscular Disorder Middle: In myopathy patients, the number of functional muscle fibers per MU decreases, leading to short-duration, small-amplitude, and polyphasic MUAPs. Right: In neuropathy patients, the number of muscle fibers in the MU increases, which causes long-duration, highamplitude, and polyphasic MUAPs [1].

A. iEMG Dataset
The iEMG recording was monitored visually and audibly by means of a computer monitor and a speaker. The amplitude and duration of an MUAP captured by the needle electrode were measured on an oscilloscope panel. An MUAP was selected if it was more than a threshold and generated by a single MU. No more than three different MUAPs were recorded at a single insertion site within an epoch of 50 ms. The signals, which were not too complex or noisy and thus suitable for evaluating MUAP parameters, were recorded for a period of 11.2 s. Care was taken to record only once from an MU and to ensure that the entire muscle was explored. The recorded iEMG signals were amplified 4000 times with a differential amplifier and filtered by an analog band-pass filter with a lower and upper cut-off frequency of 2 Hz and 10 kHz, respectively. The signals were then sampled at a rate of 23 438 Hz and digitized with a 16-bit resolution.

B. Lifting Wavelet Transform
4) Normalization or scaling. The outputs of the lifting scheme are weighted by the factors k e and k o , where k e is the inverse of k o . These values are used to normalize the energy of the scaling and wavelet functions, respectively.
1. Construct k new times series x k o defined as: where o = 1, 2, . . . , k denotes the initial time value, k represents the discrete time interval between points, and · returns the integer part of the argument. 2. Compute the average length of each time series x k o as given by: 3. Calculate the average length of all x k o 's having the same delay k as: 4. Estimate FD (denoted as D) to be the slope of the least square linear fit of the curve, ln(L(k)) versus ln(1/k), which implies that L(k) ∝ k −D . Therefore, when L(k) is plotted against 1/k, with k = 1, 2, . . . , k max , on a double logarithmic scale, the data must remain on a straight line with a slope equal to D. Scale parameter (k max ) Fractal dimension (D) For the LWT coefficients in each subband from level-five decomposition, D was computed with Higuchi's method by selecting an integer value for k max in the interval [2,50]. Around k max = 10, D plateaus meaning that there is no benefit from further calculations, i.e., computing D with k max > 10.
In order to select a suitable value for the scaling parameter k max , D values computed for the wavelet subband coefficients were plotted against k max ∈ {2, 3, . . . , 50}. The value of k max for which D reaches a plateau is regarded as the saturation point as suggested in [4], [5]. In our case, k max was selected as 10 as shown in Fig. 2.

D. One-Dimensional Local Binary Pattern
For a given data point s[n] in a signal, the 1-D LBP operator is defined as an ordered set of binary comparisons between s[n] and its neighboring S/2 data points occurring before and after s[n]. For our application, four neighboring data points are considered before (s [1], . . . , s [4]) and after (s [5], . . . , s [8]) each data point s[n]. The 1-D LBP code is constructed for s[n] by comparing its value with each one of its neighboring data points, S = {s [1], s [2], . . . , s [8]}, and thresholding the respective differences, Γ (s[k] − s[n]), k = 1, 2, . . . , card(S), to produce a set of binary numbers Thus, an eight-digit binary code is obtained for a neighborhood S surrounding a data point s[n]. The decimal equivalent of the binary code ranging between 0 and 255 is then obtained as given by where card(·) is the cardinality of the set in the argument. The LBP code in (3) represents the local structural information around s[n]. After sequentially performing this operation for all the data points in a signal, a histogram is constructed with the decimal equivalents of the LBP codes, which serves as the LBP feature for our iEMG classification task. BMMV inputs: 9 downsampled signal labels BMMV output: iEMG label Fig. 3. An illustration to demonstrate the class assignment by the BMMV algorithm. Four iEMG recordings are arranged along the rows, each of which is downsampled by a factor nine, i.e., M = 9. Each row of colored squares (left) represents a set of nine class labels corresponding to "disjoint" downsampled signals as determined by the MLPNN. The colored hexagon in the respective row (right) denotes the class assigned by the BMMV to an iEMG datum based on the majority vote rule. The filled-in colors imply the following categories: red ⇒ healthy, green ⇒ myopathy, blue ⇒ neuropathy, and black ⇒ indeterminate.

A. Choice of Wavelet Function for LWT
Since TF characteristics of wavelet functions vary, it is imperative to select an appropriate wavelet function to perform LWT. This means that wavelet coefficients reflect the degree of similarity between the wavelet function and the decomposed iEMG signals, i.e., larger wavelet coefficients imply a higher degree of similarity [6]. The downsampled iEMG signals are nonstationary and nondecaying as illustrated in Fig. 4. The mother wavelets in Fig. 5, namely, db4, sym5, coif2 Fig. 4. A section of a downsampled iEMG recording (first row). The frequency subbands, a 5 , d 1 , d 2 , . . . , d 5 , from level-five LWT decomposition of the downsampled iEMG signal contain scaling/wavelet coefficients (rows two to seven). From the shape of the signal in the first row, one can infer that the downsampled iEMG signals are nonstationary and nondecaying.     6. The best Ac value returned by a grid search (Ac gridBest ) for each wavelet function, with which the LWT was implemented, is shown; the associated hyperparameter values are displayed in the right panel. Note that rbio3.7 wavelet yielded the maximum Ac gridBest (90%), and hence we selected rbio3.7 for our implementation.
To further narrow down the choice, fused feature sets-FDs of LWT subband coefficients computed with a chosen mother wavelet and 1-D LBP-from the validation iEMG dataset were supplied to MLPNN. For each wavelet function, the best accuracy achieved by the classifier (Ac gridBest ) via a grid search was noted. Among the Ac gridBest values, the maximum was attained by rbio3.7 wavelet function (see Fig. 6), and therefore we chose rbio3.7 wavelet for our classifier framework. The decimation factor M was determined to be nine by the following empirical procedure. The iEMG signals from the validation dataset were downsampled by an integer factor, which was varied from two to 15 in steps of one in every instance. The fused feature sets constructed from the downsampled signals were fed to MLPNN performing 70-30 hold-out validation.

B. Determining Decimation Factor
The hyperparameter values were varied in a grid search. The Ac gridBest values were noted for M values in the interval [2,15], whose maximum was attained for M = 9.
To evaluate how much the performance of MLPNN-BMMV classifier (particularly the BMMV algorithm) is reliant on the choice of M, we conducted the following investigation. The iEMG signals in the experimental dataset were decimated by the integer values in the range [2,15]. For each M value, the features extracted from the downsampled signals were classified using the MLPNN-BMMV classifier (10-fold CV) and the respective Ac (%) was recorded. The empirical results are plotted in Fig. 7, which show that the BMMV-based classifier performance is almost insensitive to the choice of M ≥ 5.

C. Performance Comparison of Binary iEMG Classifiers
The following studies have reported the classification performance measures, i.e., Ac, Se, and Sp, obtained from a relatively simple binary classification 1 -healthy or ALS-using the same iEMG dataset as ours. Interestingly, the MLPNN-BMMV approach resulted in a higher Ac value (99.87%) than those listed in Table I.