Denoising electromyographic signals via stationary wavelet decomposition and filtering
Rehman Ali, Gareth Guvanasen & Stephen Deweerth, PhD
Abstract
The measurement of electromyographic (EMG) activity may be used to diagnose neuromuscular disorders, to control electrical devices through muscle-computer interfaces, or to return motor function via neuroprosthetic devices. Unfortunately, EMG measurements will often contain noise, introduced by mechanical vibrations, electromagnetic radiation, or the electrical stimulation of tissue. We compare wavelet-based filtering to traditional Fourier-based filtering methods and demonstrate that wavelet-based approaches better conserve EMG signals in the frequency domain while removing noise measurable within the time domain.
Introduction
Electromyographic (EMG) signals are composed of spatially and temporally superimposed motor unit action potentials (MUAP) from various muscle motor units. The contractile force of an isometrically held muscle can be determined by the frequency of these MUAPs, however, signal noise can obscure the identification of the firing frequencies of these motor units (Ding, Wexler, & Binder-Macleod, 2002). Furthermore, chronic low back pain (CHLBP) and various pathologies such as motor dysfunction and neuromuscular junction disorders can be identified based on altered shape and frequency of these MUAPs (Huppertz, Disselhorst-Klug, Silny, Rau, & Heimann, 1997). Each MUAP can be mathematically interpreted as a finite energy signal localized in time and encoded in frequency. Based on this interpretation, if the time or the frequency domain of an EMG signal are considered separately, information would be lost.
Traditional Fourier transforms limit our interpretation of the EMG to the frequency domain of the input signal (McClellan, Schafer, & Yoder, 2003). Thus, if a narrow band in the frequency domain were isolated, resolution about the time interval over which the narrow band of frequencies was active would be lost. Conversely, if we analyze a signal using narrow windows of time, we lose frequency resolution over those time intervals. Therefore, it is impossible to isolate a narrow range of time and a narrow band of frequencies at which the signal occurs. Fourier analysis, by itself, cannot simultaneously render information about the localized energies in the time domain and the frequencies expressed by the signal. An inability to process EMG regarding both of these aspects of MUAPs makes noise removal, and, consequently, EMG analysis difficult. Therefore, a we need a processing algorithm that perform both time-domain and frequency-domain analysis.
A promising approach is wavelet analysis, which provides information in the form of scale and time shift, as opposed to only frequency in the case of Fourier analysis. Wavelet analysis processes a signal by decomposing its constituent frequencies into logarithmically spaced bins in the frequency domain. Therefore as the resolution of frequency is increased, the time resolution is reduced, and vice-versa. Despite this trade off in resolution, wavelets can be optimized for isolating time and frequency bands to identify MUAPs of specific energies and dilations localized in time. The wavelet transform has been used to remove noise from EMG signals. Wavelets such as the Mexican Hat and the Morlet, in particular, have been used to analyze EMG signals because of their resemblance to MUAPs (Laterza & Olmo, 1997).
We will demonstrate the results of noise removal by wavelet and Butterworth (a Fourier-based filter) filtering on EMG data in order to determine the tradeoffs of each method on EMG. This will be accomplished in the time domain by correlating the derivatives of the filtered EMGs (wavelet and Butterworth) and in the frequency domain by quantifying the deviation in the derivative of the decibel spectra for each signal. The filtering method that demonstrates the greatest correlation of time-domain derivatives and minimizes the deviation in the frequency domain would be shown to best conserve the original EMG signal.
Methods
To explain our application of wavelets, we present the mathematical theory behind wavelets and their application under the discrete wavelet transform (DWT). Equation (1) is the basic definition of a child wavelet. The mother wavelet is , and the children wavelets are time-shifted (b is time-shift) and normalized dilations (a is dilation) of . The defining property of is that its integral over all space is 0, the integral of its square over all space is 1. Because the child wavelets inherit these properties, they are the ideal finite energy signals on which a signal may be projected. Discretization (2) of a and b into j and k enables the wavelet decomposition of x(t) to be expressed in series form.
The series coefficients are calculated by orthogonal projection in equation (3) and x(t) can be reconstructed by using the series in equation (4). In addition to the mother wavelet, the scaling function (or father wavelet) was devised to decompose the residual function, because it would take an infinite number of terms to describe series (4). The key properties of are that it integrates and square integrates to 1 over all space. Decomposition and reconstruction formulas similar to (3) and (4) are applied to the residual function in terms of 's. These derivations give the foundation of the continuous wavelet transform (CWT), the continuous time form of the DWT.
Wavelet Filter Design Theory
Wavelet decomposition in discrete time occurs in the DWT tree in Figure 1. The low-pass filter g[n] and high-pass filter h[n] come from the decompositions of a mother wavelet and a father wavelet into father wavelets that are narrowed by a factor of 2. This is why outputs of the filters are then downsampled by a factor of 2 on the nodes of the decomposition tree in Figure 1. The g[n] and h[n] satisfy equations (5) and (6):
According to the definitions in equations (5) and (6), the mother wavelet definition actually comes from the father wavelet. Equations (5) and (6) also imply relations (7) and (8):
Furthermore, the DWT tree in Figure 1 can be reversed by upsampling at the nodes and passing the decomposition through the reconstruction filters whose coefficients are simply in reverse order of the deconstruction filters. The wavelet-based filtering method applied in this work involves scaling each decomposition level to remove wavelets of certain dilations while maintaining others, perhaps with enhancement, and then reconstructing the signal with the modified coefficients by traversing through the aforementioned reconstruction tree.
Stationary Wavelet Transform Decomposition of a Signal
A wavelet decomposition using the Haar wavelet with one approximation level and five detail levels illustrates DWT decomposition (Figure 2C). The mother Haar wavelet is shown in Figure 2A and the corresponding father wavelet is shown in Figure 2B. The d_{1} level is composed of the most binary-narrowed (narrowed by factor of a power of 2) and time-indexed Haar wavelet. Going up the detail levels from d_{1} to d_{5}, the Haar Wavelet dilates by a factor of two at each stage, making the blocks from the Haar wavelet become more apparent. The graphs in d_{1} to d_{5} integrate to zero because integrates to 0, but the approximation level integrates to some positive number because integrates to 1. This observation will explain why the wavelet decomposition produces local offsets in a filtered signal
One disadvantage of a pure DWT is that time shifting is discrete and dependent on the degree of dilation of the wavelet at that level, so the decompositions of a periodic signal such as the ECG in Figure 3 are not necessarily periodic. The first six ECG signals in Figure 3 are almost perfectly periodic, but their corresponding decompositions are not periodic. This means that discrete time-shifting produces phase compensation issues at each level. A modification of the DWT known as the stationary wavelet transform (SWT) overcomes this lack of translational invariance by upsampling the filters rather than downsampling the signal at each level. The SWT tree (Figure 4) demonstrates the method of decomposition and reconstruction used in this work.
Filtering Techniques and Applied Parameters
DC components were removed from the original EMG signal obtained from the MEA in order to specifically address EMG morphology rather than voltage offset. To serve as an example of traditional Fourier-based filtering, a 20^{th} order non-causal (centered) low-pass Butterworth filter with normalized cutoff frequency of 0.2 was applied to the EMG signal. To demonstrate wavelet-based filtering, the Stationary Wavelet Transform (SWT) of the EMG signal was performed by using the coif3 wavelet to decompose the EMG signal down to ten levels (or scales) in the wavelet domain. The coif3 wavelet is often used to decompose EMG because it resembles the MUAP (Liu, Chen, & Chen, 2005). Each scale was given a certain weight: one weight for the residual approximation level at the 10^{th} scale (a multiplication factor of 1.2) and ten weights for the detail coefficients at each scale (in ascending order from the 1^{st} to the 10^{th} scale: 0.135, 0.675, 0.945, 1.215, 1.35, 1.35, 1.35, 1.35, 1.35, and 1.35). The weights were given the aforementioned values because they most effectively removed noise. A wavelet-filtered EMG signal was formed by applying the inverse SWT on each of the weighted scales.
Time-Domain Analysis
The filtering method that better conserves the time-domain characteristics of the signal must best conserve the shape of the original EMG and consequently its derivatives after filtering. Thus, to determine which filtering method better conserved the original EMG signal after increasing orders of differentiation, a Spearman's correlation was calculated between the points of the original EMG and each of the filtered EMGs, and between the points of each of their higher order derivatives as well. Furthermore, p-values for each correlation values were obtained in order to evaluate statistical significance. The filtering method that demonstrates a higher statistically significant correlation at each derivative would be deemed to best conserve the time-domain content of the original signal.
Frequency-Domain Analysis
Time-domain analysis provides an incomplete evaluation of how well a particular filtering method removes noise, so it is necessary to examine the frequency domain for a more comprehensive assessment. Single-sided spectrums were generated using the Fast Fourier Transform (FFT) for the original EMG, Butterworth-filtered EMG, and wavelet-filtered EMGs. The amplitudes were represented on a decibel scale to accentuate details in the frequency domain. The relative amplitude content within the frequency domain is best described by the shape of the decibel FFT spectrum. Thus, the derivative of the amplitude in decibels with respect to frequency in hertz was used to quantify shape. Root mean square (RMS) error values were calculated between the derivatives of the FFT spectra for the original EMG and the wavelet-filtered EMG as well as between derivatives of the FFT spectra for the original EMG and the Butterworth-filtered EMG. The lower error value would indicate which method best conserves frequency-domain information.
Results
The SWT-filtering and Butterworth-filtering methods were compared in the time and frequency domain by correlation analysis and RMS error respectively.
Impact of the Filters in the Time Domain
The selected segments of the original and filtered EMGs are shown in Figure 5. We observe that at localized regions within the EMG signal, the stationary wavelet transform (SWT) produces a signal that best matches the shape of the input signal. Although the low-pass Butterworth filter effectively removes noise from the original EMG, it produces a smoothed version of the original EMG and significant shape characteristics of the original signal are lost. Although the SWT reconstructed signal is slightly offset elsewhere, it retains shape characteristics better than the Butterworth filter.
Figure 6 illustrates that for the first derivative and higher, the correlation of the SWT-filtered EMG to the original EMG was consistently higher than the correlation of the Butterworth-filtered EMG to the original EMG. Thus, SWT-filtered EMG signal better conserves the original EMG when taking progressively higher order derivatives than the Butterworth-filtered EMG. Similar results were obtained with regards to the qualitative shape characteristics and the correlations over the derivatives when these filtering methods were applied over other intervals of EMG.
Impact of the Filters in the Frequency Domain
The FFT spectra in decibels and the derivative of the FFT spectra with respect to frequency are shown in Figure 7A. Although both the Butterworth filter and the SWT decomposition attenuate frequencies outside the baseband of the EMG, the SWT decomposition better conserves the content of the entire frequency spectrum (within and outside the baseband). The preservation of the original EMG spectrum become more noticeable in Figure 7B, which show the derivative of FFT spectra of the Butterworth-filtered, SWT-filtered, and original EMG signals. A comparison of the spectra demonstrates that the stationary wavelet transform better conserves the shape characteristics of the original signal without perturbing relative amplitude content over the entire frequency domain.
Despite the frequency attenuation, the SWT derivative in the frequency domain best matches the frequency-domain derivative of the original EMG (RMS error of 0.33); however, the Butterworth-filtered EMG derivative does not overlap with the original beyond the cutoff frequency for the Butterworth filter (RMS error of 4.48). As a result, the Butterworth filter produces an RMS error that is 13.6 fold greater than through the SWT. Thus, greater frequency detail is lost via the Butterworth filter than through the SWT. Similar results are obtained when these filtering techniques applied to other intervals of EMG.
Discussion
SWT filtering was shown to better conserve the content of the EMG by producing a greater correlation at each derivative and a lower RMS error on the frequency-domain derivative, but the SWT filtered signal also appears to have introduced baseline wander to the original EMG signal. Regardless, the properties of SWT filtering should enable improved characterization of the MUAPs contained within an EMG.
Source of Baseline Wander
The localized offsets produced by the wavelet-filtered EMG in the time domain is called baseline wander, which the stationary wavelet transform appears to modify in the original signal. The reasons for this departure originate from the convergence properties of the wavelet space: as the wavelet decomposition converges, less baseline modification is observed. This space for wavelet decomposition L^{2}( ) is defined as the sequence of finite energy signal subspaces (Misiti, 2007). As j decreases, the subspace is inclusive of all the 's that came before it (higher j's) and is used to make progressively better approximation of x(t). The element that adds more detail to is the child wavelet narrowed by factor of two: . Thus, is composed of and the incremental level of detail , which is the subspace of wavelets of a certain dilation.
Convergence, thus, depends on diminishing levels of detail coming from as j increases. The recursive set definition from (9) means the following after iterative dilations in discrete time:
Because of sampling from continuous time to discrete time, the wavelet cannot be infinitely narrow; therefore, denotes the level of the most narrowed wavelet in discrete time. The set notation (10) describes the Wavelet decomposition performed up to n levels. In Figures 2 and 3, the approximation level a _{5} correspond to , and the detail levels d_{1} to d_{5} correspond to to The nth level residual or the approximation level can be decomposed into the scaling functions . If , would be the DC component of a periodic signal.
These decomposition properties indicates that if wavelet decomposition were not carried out to the sufficient number of levels, would present a significant amount of baseline wander. Because the windows of analysis were 1024 sample long, the largest n could have been was 10 (2^{10} = 1024); thus still represented a significant amount of baseline wander. The only way to ameliorate this baseline wander is to expand the analysis window exponentially; however, this window expansion would lead to significant increases computation time. Thus, our choice of n = 10 was based on striking the balance between the convergence of and minimizing computation time.
Characterizing Shape by Derivatives
Software packages that employ wavelet decomposition on EMG have been shown to identify and extract individual MUAPs based on their shape characteristics (Zennaro, Wellig, Koch, Moschytz, & Laubli, 2003). Because the shape of the original EMG is best defined by its derivatives, a Spearman's correlation was conducted between the points of the original EMG and the filtered EMGs (Figure 6). It was discovered that for the first derivative and higher, the correlation of the SWT-filtered EMG to the original EMG was consistently higher than the correlation of the Butterworth-filtered EMG to the original EMG. This preservation of time-domain derivatives may be of great use in analyzing muscle mechanics through EMG because in the literature we find that neural activation of muscle is best defined by the derivatives of the EMG signal rather than the EMG signal itself (Buchanan, Lloyd, Manal, & Besier, 2004).
Our method of determining the efficacy of wavelets deviates from the convention, which involves minimizing the square of the error and having the wavelet decomposition provide better convergence to the original signal (Daubechies, 1990). This convention, however, does not provides information regarding the conservation of the derivatives of the EMG. Furthermore, it assumes that convergence to the original signal determines efficacy of the method. It is because wavelets conserve the key features within a signal that methods such as the wavelet functional ANOVA have been developed to identify statistically significant differences between EMGs (McKay, Welch, Vidakovic, & Ting, 2013).
Characterizing Shape in Frequency Domain
The frequency domain was examined because time-domain analysis alone remains an incomplete evaluation of how well a particular filtering method removes noise. As such, frequency-domain analysis was used to supplement the aforementioned time-domain analysis. It was discovered that the SWT better conserved the frequency domain of the original EMG while attenuating the frequencies outside the baseband. The weights given to the levels from d_{1} to d _{10} compress the signal at certain scales (or dilations) rather than frequencies. As a result, wavelets enable us to remove noise at certain frequencies without removing EMG information at those frequencies. Furthermore, since scale is related to logarithmically spaced bins in the frequency domain, wavelets can also be exploited to attenuate a signal at certain frequencies just as in this work.
Future Work
The work presented herein does not provide a complete analysis of Fourier-based and wavelet-based filtering techniques for analyzing EMG signals. Further research should attempt to demonstrate which methods among Fourier-based and wavelet-based filters are the most effective with regard to time-domain derivative correlation and frequency-domain error minimization. Many wavelet-based filtering methods simply involve giving weights to the coefficients at each scale just as presented in this work, but other wavelet-based filtering techniques, common in the field of data compression, are based on applying separate thresholds to the coefficients at each scale in order to reduce low energy wavelets in the original signal (Santoso, Powers, & Grady, 1997). Avenues such as these need to be thoroughly explored to compare the effects of various wavelet-based filters on the time and frequency domain. For example, "soft thresholding" and "hard thresholding" are still subcategories under the aforementioned threshold-based method and should also be investigated. The methods of analysis presented in this work should be applied to different kinds of wavelet-based filtering approaches in the future to obtain a more complete understanding of their impact on the signal.
References
- Buchanan, T. S., Lloyd, D. G., Manal, K., & Besier, T. F. (2004). Neuromusculoskeletal modeling: Estimation of muscle forces and joint moments and movements from measurements of neural command. Journal of Applied Biomechanics, 20(4), 367-395. Daubechies, I. (1990). The Wavelet Transform, Time-Frequency Localization and Signal Analysis. IEEE Transactions on Information Theory, 36(5), 961-1005. doi: Doi 10.1109/18.57199
- Ding, J., Wexler, A. S., & Binder-Macleod, S. A. (2002). A predictive fatigue model--I: Predicting the effect of stimulation frequency and pattern on fatigue. IEEE Trans Neural Syst Rehabil Eng, 10(1), 48-58. doi: 10.1109/TNSRE.2002.1021586
- Huppertz, H. J., Disselhorst-Klug, C., Silny, J., Rau, G., & Heimann, G. (1997). Diagnostic yield of noninvasive high spatial resolution electromyography in neuromuscular diseases. Muscle Nerve, 20(11), 1360-1370.
- Laterza, F., & Olmo, G. (1997). Analysis of EMG signals by means of the matched wavelet transform. Electronics Letters, 33(5), 357-359. doi: Doi 10.1049/El:19970250
- Liu, H. H., Chen, X. H., & Chen, Y. G. (2005). Wavelet transform and real-time learning method for myoelectric signal in motion discrimination. 7th International Symposium on Measurement Technology and Intelligent Instruments, 13, 250-253. doi: Doi 10.1088/1742-6596/13/1/058
- McClellan, J. H., Schafer, R. W., & Yoder, M. A. (2003). Signal processing first. Upper Saddle River, NJ: Pearson/Prentice Hall. McKay, J. L., Welch, T. D. J., Vidakovic, B., & Ting, L. H. (2013). Statistically significant contrasts between EMG waveforms revealed using wavelet-based functional ANOVA. Journal of Neurophysiology, 109(2), 591-602. doi: DOI 10.1152/jn.00447.2012
- Misiti, M. (2007). Wavelets and their applications. London ; Newport Beach, CA: ISTE. Santoso, S., Powers, E. J., & Grady, W. M. (1997). Power quality disturbance data compression using wavelet transform methods. IEEE Transactions on Power Delivery, 12(3), 1250-1257. doi: Doi 10.1109/61.637001
- Zennaro, D., Wellig, P., Koch, V. M., Moschytz, G. S., & Laubli, T. (2003). A software package for the decomposition of long-term multichannel EMG signals using wavelet coefficients. IEEE Transactions on Biomedical Engineering, 50 (1), 58-69. doi: Doi 10.1109/Tbme.2002.807321
Acknowledgements
The author would like to acknowledge the mentorship of graduate student Gareth S. Guvanasen working under Dr. Stephen P. Deweerth.