Neuromuscular disease detection based on feature extraction from time–frequency images of EMG signals employing robust hyperbolic Stockwell transform

In this paper, a novel technique for detection of healthy (H), myopathy, (M) and amyotrophic lateral sclerosis (ALS) electromyography (EMG) signals is proposed employing robust hyperbolic Stockwell transform (HST). HST is an efficient signal processing technique to analyze any nonstationary signal in joint time–frequency (T–F) plane. However, a major issue with HST is the optimum selection of Gaussian window parameters since the resolution in the T–F plane depends on the shape of the window. Considering the aforesaid fact, in this article, a genetic algorithm (GA) based optimized HST is proposed for improved EMG signal analysis in T–F plane. Several novel features were extracted from HST spectrum and features with high statistical significance were selected for classification using several benchmark classifiers. It was observed that optimized HST resulted in better classification accuracy of EMG signals which indicates its potential for clinical applications.


| INTRODUCTION
Electromyography (EMG) is a graphical recording of complex and nonstationary electrical signals generated by the motor units innervating the skeletal muscles and associated motor neurons that control voluntary and reflexive muscle actions. EMG signals are widely used in pathology labs to diagnose several neuromuscular and neurological disorders, among which myopathy (M) and amyotrophic lateral (ALS) are very common. Myopathy is a disease associated with muscles causing symptoms such as muscular spasm, muscle cramps. 1 Amyotrophic lateral sclerosis (ALS) is a highly progressive neurological disorder that is responsible for the death of motor neurons. 2 Therefore, the ability of the human brain to control voluntary muscle actions is lost, resulting in fatigue, numbness, loss of reflex action, and so forth. If left untreated, ALS can be sometimes fatal. A significant number of unnatural deaths resulted from ALS have been reported. Therefore, it is necessary to accurately detect these diseases early to avoid further progression. It is not easy task to differentiate between M and ALS from each other based on medical symptoms only. In clinics and hospitals, detection of M and ALS disorders is usually done by visual screening of EMG data by expert neurologists. However, this conventional method is unsophisticated and unreliable since it relies solely on human intervention and is subjected to human error, especially for longer EMG recordings. To address this issue, researchers all over the world are trying to develop a computer-aided fast and accurate disease detection system using modern signal processing and machine learning techniques so that the diagnosis of neuromuscular and neurological disorders can be done quickly and accurately by extracting meaningful features from EMG signals.
In existing literature, several methods have been proposed by researchers for detection and classification of neuromuscular disorders based on EMG signals. Time domain analysis of EMG data using cross correlation has been implemented for detection and classification of neuromuscular diseases in Reference 3. Several time domain features like zero crossing, root mean square, autoregressive cepstral analysis for detection of EMG signals have been reported in References 4-6. Detection of abnormal EMG signals using time and frequency domain features has been reported in References 7,8. Classification of EMG signals using fast Fourier transform (FFT) has been reported in Reference 9. In References 10-12, Bajaj et al. proposed empirical mode decomposition (EMD) based musculoskeletal disease detection and classification framework. Classification of EMG signals employing Stockwell transform (ST) has been reported in Reference 13. In Reference 14, convolutional neural network (CNN) based approach was proposed for the classification of neuromuscular disorders from EMG signals. Classification of EMG signals corresponding to different healthy and diseased conditions employing nonlinear dynamic method has been reported in Reference 15. Optimally parameterized weighted visibility graph (WVG) has been proved to be an effective tool for automatic diagnosis of neuromuscular diseases. 16 Since EMG signals represent firing patterns of motor neurons, they are recurrent and nonstationary in nature. Therefore, EMG signal analysis in joint time-frequency (T-F) domain provides a better understanding of dynamics of motor neurons. Short-time Fourier transform (STFT) based T-F domain analysis for detection of myopathy from normal EMG signal has been reported in Reference 17. Detection and classification of EMG signals using continuous wavelet transform (CWT) have been reported in References 18,19. In References 20-22, researchers proposed a discrete wavelet transform (DWT) based automated diagnosis system for neuromuscular disorders. Cross wavelet transform (XWT) based disease detection framework from EMG signals has been reported in Reference 23. In References 24,25, a novel tunable Q-factor wavelet transform (TQWT) based EMG signal classification framework had been proposed. Detection of M and ALS employing modified window Stockwell transform (MWST) has been reported in Reference 26. Thus, it is evident that several techniques for analysis of nonstationary EMG signals in time-frequency frame have been reported for classification of neuromuscular disorders.
In joint time-frequency (T-F) analysis, ST is a popular technique that has been used to analyze nonstationary biomedical signals in exiting literature. 26,27 However, in the generalized version of ST, a fixed Gaussian window is used, which is not adaptive in nature. To overcome this limitation, HST is proposed in Reference 28 as an extension of the ST. HST is proposed by replacing the standard Gaussian window with a hyperbolic Gaussian window, where the front and the back taper of the Gaussian window can be varied depending on the type of signal to be analyzed. Considering its advantages, HST has been widely used to analyze many nonstationary time series in T-F plane. 28,29 However, one major issue associated with HST is the choice of appropriate hyperbolic window (HW) parameters. In T-F analysis, proper selection of Gaussian window parameters is extremely important since the resolution in T-F plane depends on the shape of the Gaussian window. 26,30 Improved resolution in T-F plane can lead to better discrimination of spectral images of nonstationary signals. This in turn leads to better feature extraction and eventually classification of EMG signals. In previous works, window parameters of the hyperbolic function have been chosen empirically, depending on the type of signal. 29 The process is not automated and adaptive in nature. Considering this issue, a robust HW parameter selection based on evolutionary computation technique is proposed in this article using maximum energy concentration measure (ECM) based constrained optimization. The parameters should be chosen in such a way such that the highest ECM in the T-F plane can be achieved. The main contributions of the proposed work are as follows: i. HST based analysis of EMG signals in T-F plane for classification of neuromuscular disorders is reported for the first time. ii. Besides, instead of arbitrary selection of HW parameters using trial-and-error method, a method based on GA approach is proposed for robust selection of HW parameters such that resolution of EMG signals in the T-F plane is maximum.
iii. From the transformed EMG signals in the T-F plane, several new features were extracted and Student's ttest of the features was conducted between each class to examine their statistical significance.
iv. Finally, using the highly statistically significant features, the performance is evaluated using three classifiers.
v. Four classification tasks are performed by using different combinations of EMG signals of each class.
A brief outline of our work is illustrated using a flowchart in Figure 1.

| EMG DATABASE
In the present study, EMG signal are acquired from the publicly available online EMGLAB database. 31 The database comprises of six male and four female healthy (H) subjects between 21 and 37 years of age. None of them have any history of neurological disorders. In addition, the database also contains EMG recordings from five male and two female myopathy (M) patients between 19 and 63 years of age and four male and four female ALS patients aged between 35 and 67 years, respectively. Recording of EMG data was done from five different places in the muscle by inserting a concentric needle electrode. The electrode was inserted at three different insertion depths and at constant and low voluntary levels of muscle contraction employing audio-video feedback. The EMG data were digitized at a sampling frequency of 23437.5 Hz using a 16-bit resolution digitizer. After acquisition of EMG signals from H, M, and ALS patients, they were filtered using a digital high-pass filter (firstorder Butterworth, 50 Hz, 3 dB limit, with zero-phase distortion) to suppress the low-frequency baseline movements. In addition, the signals were further filtered using a high and low pass filter with cutoff frequencies set at 2 Hz and 10 kHz, respectively prior to analysis. The detailed description of the database has been reported in Reference 31. Figure 2 shows the snapshot of sample H, M, and ALS EMG signals, respectively.

| Description of tasks
In this work, four classification tasks have been performed. In Table 1, detail of those classification tasks is shown.
The main significance of selecting those four different classification tasks is to deeply investigate the performance of proposed paradigm in detection of neuromuscular diseases. The purpose of performing T-I and T-II is to distinguish both M and ALS from H EMG signal. The purpose of T-III is to discriminate between M and ALS signals. T-IV is a general classification task where it is possible to distinguish between H and disease (either M or ALS) signals, without identifying the exact type of neuromuscular disease. In the next section, we will discuss about the methodology in detail, which is used in this paper for EMG classification.

| Hyperbolic Stockwell transform
Generalized version of ST of any function g(t) can be mathematically expressed as Equation (1) 28 : Here, f is the Gaussian window (GW), whose standard deviation σ(f ) varies inversely as frequency f, and τ is the translation parameter. Due to fixed shape of the GW, generalized version of ST offers high time resolution for low-frequency components and higher frequency resolution for high-frequency components, which imposes problem in dealing with nonstationary signals. To address this issue, hyperbolic Stockwell transform (HST) was developed as an extended version of generalized ST, which replaces the generalized GW with a modified HW whose shape can be varied by varying the window parameters. The HST of any function g(t) can be expressed as Equation (2) 29 : In the above equation, γ f and γ b are forward and backward taper parameters, respectively, λ is the positive curvature parameter, which controls the shape of the window, and V is the hyperbola, which is expressed as Equation (3): In Equation (3), ξ is the translation factor which is used to make certain that the peak of HW occurs at (τ À t) = 0 and is expressed as Equation (4):

| Selection of HW parameters
The HST described in Equation (2) has three distinct parameters (γ f , γ b , and λ) that control the shape of the HW. Interestingly, the resolution in T-F plane can be altered by varying these three parameters. In existing studies, the parameters of the HW are usually selected by trial-and-error method depending on the nature of the signal to be analyzed. In this contribution, a method for robust window parameter selection based on GA approach is proposed in such a way that the chosen parameters of the HW can provide maximum energy concentration measure (ECM) in the T-F plane. The objective function of the present optimization problem is given by Equation (5): Here, s.t Equation (6) is the normalized HST. In the earlier works by the authors, similar type of objective function was formulated to optimize window parameters of modified ST in Reference 26. Inspired by Reference 26, similar optimization problem is formulated here to optimize window parameters of HST. Also, in the earlier works, this is the standard fitness function used to determine the resolution of any nonstationary signal in the T-F plane. 26,32,33 Once the optimum values of γ f , γ b , and λ are determined, the optimum parameter ξ is computed using Equation (4). Now, the constraints of the optimization used in this work are based on the criterion given in Reference, 28 and f p is the peak frequency of the signal to be analyzed. As mentioned earlier, GA-based optimization technique has been utilized in this work to optimize the HST spectrum. GA is a very popular optimization algorithm, proposed in. 34 GA has proved to be a very efficient and stable searching algorithm for global optimization problems with three major parts: generation, genetic operation, and replacement. 35 In Figure 3, the optimization process utilizing this algorithm has been reported. GA has been successfully used in different domain such as mechanical engineering, 36 flowshop scheduling, 37 finding multiple shortest routes in street networks, 38 cloud computing, 39 signal processing, 35 and biomedical engineering 30 to solve complex optimization problems. In this work, GA has been used to obtain the optimum HW parameters that correspond to maximum ECM of the HST spectrum in the T-F plane. A flowchart showing the use of GA for optimum selection of window parameters is given below.

| Advantages of HST over existing methods
In existing literature, several techniques have been proposed by researchers for analysis of EMG signals in the T-F frame. Compared with existing methods, HST has certain advantages. In comparison with CWT, HST can retain the absolute phase information of the signal and is less sensitive to noise. 27 Also, unlike CWT, signal analysis using HST is dependent on the type of mother wavelet used. Moreover, in CWT, the shape of the mother wavelet is fixed irrespective of the type of signal to be analyzed. On the other hand, HST offers the provision to adaptively vary the front and back taper of the HW depending on the nature of the input signal. Thus, HST is more signal adaptive compared with CWT. In addition, compared with EMD, signal analysis using HST is free from modemixing and end effect problems. Thus, considering its advantages, HST with tunable window parameters is proposed in this work to analyze H, M, and ALS signals in the T-F plane. To determine the optimum window parameters, GA is used in this work such that it results in maximum ECM in the T-F plane. From the transformed signals in the T-F plane, several features have been extracted, which are described below.

| Feature extraction
In this present work, 14 novel features have been selected for classification purpose. Application of HST transforms any time series to a matrix consisting of complex elements. In this study, the features were extracted from the magnitude of the HST matrix. In Table 2, selected 14 features have been reported.

| Statistical test
In any classification task, statistical tests play a vital role in feature ranking, selection, and elimination in higherdimensional features. This not only reduces the processing time but also increases the classification performance, which may get compromised because of redundant features. Among different statistical tests, ttest has been proved to be one of the most popular and efficient tests that can be used to compare the means of two groups. 40 In this present work, t-test based statistical analysis of features has been done to examine the statistical significance of the extracted features between different groups of EMG signals prior to classification.

| Classifiers
To classify EMG signals based on HST based features, three different machine learning (ML) classifiers were used in this study. A brief theory of the classifiers used in the present work is described below.
3.6.1 | Support vector machines (SVM) SVM is a well-known ML algorithm, which has been used by researchers in different fields, mainly to solve binary classification tasks. Since, four binary classification tasks are chosen in this study, SVM is used in this paper. In an SVM, data points from different classes are classified by finding an optimum separating hyper-plane that has maximum margin. The maximum margin is determined using structural risk minimization principle. For nonlinear SVMs, the training data are mapped to feature space of high dimension employing several kernel functions and satisfying Mercer's theorem. Some of the most frequently used kernel functions used in SVM are linear, polynomial, and radial basis functions (RBF). In this study, classification of EMG data is done using different kernel functions.

| k-nearest neighbor (kNN)
kNN classifier is a widely used supervised ML algorithm with low computational burden and easy implementation technique. This algorithm assigns class to an object depending upon the class label of its nearest neighbor. The choice of neighboring distance k makes this algorithm very efficient and easy to tune according to the data. The detailed algorithm of kNN was reported in Reference 41. Here, we used Euclidean distance as the distance parameter of kNN for classification of EMG signals.

| Naïve Bayesian (NB)
NB classifier is a very well-known ML classifier with a wide range of application. This algorithm works on Bayes' theorem and posteriori hypothesis. Owing to its probabilistic framework, NB classifier is robust against the outliers and can quickly update itself with the addition of new data points. Hence, NB classifier has been proved to be very efficient in classification of clinical data as noisy trials. Detail about the NB algorithm can be found in Reference 42. In previous literature, 16 NB classifier had been proved to be an efficient classification algorithm for EMG signal classification, and hence NB classifier is used in this present contribution.

| EMG signal analysis using HST
In this study, a total of 491 EMG signals (300 H, 106 M, and 85 ALS) have been used from this database to validate our proposed framework. In the pre-processing stage, the EMG signals were segmented into 64 equal time frame, and 25 frames (from 31st to 55th segment) from the middle steady portion of the signal have been HST operation was performed on each class of EMG signal, which transforms the time domain EMG signals into T-F frame. Now, as mentioned earlier, for the HST operation, the parameters γ f , γ b , and λ for each class of EMG signals were determined using the GA-based optimization, such that the resolution in the T-F plane is maximum. These parameters were obtained by solving Equation (5). In Figure 4, sample T-F spectrum images obtained after HST operation on three categories of EMG signals are shown. The corresponding optimum values of γ f , γ b , and λ used to compute the T-F images are also reported. From Figure 4, significant differences can be observed among H, M, and ALS signals. In Figure 5, T-F images of three classes of EMG signals obtained using the standard GW of ST are also shown for the purpose of comparison. It is evident from Figures 4 and 5 that the proposed HST with optimum selection of HW parameters has resulted in better resolution of EMG signals in the T-F frame compared with standard ST. To prove the aforesaid fact, the ECM values are obtained using the proposed optimum hyperbolic GW parameters and are reported in Table 3. Moreover, the ECM values obtained by using the standard GW of ST are also reported for the purpose of comparison. It can be seen from Table 3 that the ECM values computed using the optimum window parameters of HST spectrum is significantly higher than ST resulting better spectrum resolution in the T-F plane.
In this work, we selected the GA parameters empirically for optimum selection of the hyperbolic GW parameters. To avoid local optima, we empirically set low mutation rate (0.01) and high crossover probability (0.8) along with a population size of 50, and the initial population range was varied from 10 to 10. We also set the number of iterations at 500, so that the algorithm reaches near desired global optima instead of local optima. It was observed that the fitness function saturates near the global optima, yielding the same value with no information exchange. Now, since the goal is to achieve maximum resolution in the T-F plane, therefore the objective function as described in Equation (5) Figure 6 shows the minimization of the cost function with the number of iterations during the optimization process for sample H, M, and ALS signals, respectively. From Figure 6 we can see that the loss function decreases with each iteration and becomes almost zero at the final iteration. This optimization of hyperbolic GW yields maximum resolution of the EMG signals in the T-F plane.

| Classification performance
In this study, the time domain EMG signals were transferred to the T-F plane using HST employing optimum selection of HW parameters. After extraction of the features from the magnitude of HST spectrums, t-test has been performed, and p-values of the selected 14 features are reported in Table 4. It can be observed from Table 4 that for all classification tasks, all 14 selected features showed high statistical significance (p < 0.001).
The performance of the proposed model is evaluated in terms of three statistical measures, namely accuracy (ACC), sensitivity (SEN), and specificity (SPE). The mathematical expressions for the above three statistical measures are given in Equations (8)-(10): where TP, TN, FP, and FN symbolize true positive, true negative, false positive, and false negative, respectively. Since the number of EMG signals per class was different, prior to classification, the dataset has been augmented so that each class contains the same number of EMG signals. For this purpose, the minimum number of EMG signals present in each class is chosen. As mentioned earlier, in the present work, there are 2125 ALS signals, 2650 M signals, and 7500 H signals. So, to avoid the problem of class unbalance, 2125 signals per class is chosen in this study to evaluate the performance of the classifiers.
Here, we performed 10-fold cross validation, and finally the mean performance parameters are reported in terms of percentage (%) with standard deviation in Table 5.   Table 5, it can be clearly seen that the performance of SVM is better than the kNN and NB classifier for all classification tasks. In the case of SVM classifier, SVM with RBF kernel has delivered the best performance for T-I and T-II, whereas for T-III and T-IV, linear and polynomial kernels delivered better performance. The optimum value of RBF kernel width and the polynomial order delivering best performance is also mentioned in Table 5. Also, in the case of kNN classifier, the value of k has been varied from 2 to 15, and the k value delivering best performance is indicated in Table 5 for all four tasks. The highest accuracy has been obtained for T-III, whereas the lowest accuracy has been reported for T-I. One of the major reasons behind this observation is that the nature of the firing patterns of the motor neurons in M and ALS is well reflected by the HST (as shown in Figure 3). For M signals, EMG signals show rapid and irregular firing of motor neurons, where ALS disease causes loss of motor neurons, resulting in less firing. As M does not cause death of any motor neurons, and firing patterns of M and H appear to be quite similar (except the irregularity in M), resulting in lower classification accuracy for T-I. Nevertheless, the proposed model delivered reasonably satisfactory for all different classification tasks.

| Comparison of classification performance with ST
Classification performance of the proposed framework employing HST has also been compared with classification results, obtained by employing ST. In Figure 7, a comparative study of the classification performance between the standard ST and proposed HST is presented in terms of classification accuracy. From Figure 7, we can observe that the proposed method employing HST with optimum selection of window parameters has delivered better classification performance for all four classification tasks. This may be explained by the fact that due to the optimization of window parameters of HST, it can extract more information from the signals compared with standard ST, which employs conventional GW. Hence, the classification accuracy is high for all four classification tasks.

| Performance comparison with existing literature
The classification performance of the proposed model is also compared with some of the existing literature studied   Table 6. From Table 6, it can clearly be observed that the proposed framework outperforms most of the existing methods, and the overall performance of proposed technique is quite satisfactory. Thus, the proposed method can be potentially implemented for computer-aided diagnosis of neuromuscular disorders from EMG signals.

| CONCLUSION
In the present work, HST based feature extraction technique is proposed for automated detection and classification of neuromuscular and neurological disorders.
Instead of selecting HW parameters using trial-and-error method, a robust HW parameter selection using GAbased optimization is proposed in this article, such that the resolution in the T-F plane is maximum. From the transformed signals in the T-F plane, 14 new statistically significant features have been extracted from the HST matrix of EMG signals to distinguish between H, M, and ALS patients. Four classification tasks are addressed in this paper, and for each problem t-test of the extracted features is conducted to examine their statistical significance. Finally, the highly discriminative feature parameters are served as inputs to several benchmark classifiers for classification of EMG signals. It is observed that the robust HST proposed in this work can differentiate between H, M, and ALS EMG signals, successfully achieving a high degree of classification accuracy, which can be implemented in future for analysis of other nonstationary signal analysis. The pros and cons of the present study are discussed briefly. This is the first study where HST has been applied to analyze nonstationary EMG data in the T-F plane. Besides, optimization of HW parameters for improved resolution of EMG signal in the T-F plane is also proposed for the first time in this type of study. However, the present study has certain limitations. This study has been validated on a small number of EMG signals. Moreover, the number of EMG sample points considered for this study was 4096 only. Before clinical application, it is necessary to test the method on a large dataset of EMG signals, containing large number of sample points. This will be done as a part of the future work to develop a real-time computer-aided disease detection system. Also, the parameters of the HW were determined using GA-based optimization. In future, other optimization techniques like particle swarm optimization, simulated annealing technique, and so forth, can be implemented to determine the optimum window parameters. Besides, the proposed HST based approach has been implemented on the entire EMG signal instead of constituent MUAPs. It would be an interesting piece of work to investigate the feasibility of HST based feature extraction from MUAPs for detection of neuromuscular diseases. Further, the EMG data used in this study was acquired from biceps brachii muscle fiber. Future work will focus on classification of EMG data recorded from some other muscle fibers located in other parts of the human body. In addition, the proposed HST based feature extraction method will also be potentially implemented in future for analysis of different other nonstationary bio-potentials like ECG, EEG, and so forth, for diagnosis of cardiac as well as neurological diseases.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are openly available in 'EMGLAB' at http://www.emglab.