د.مهند كاظم ونور Al-Khwarizmi Engineering Journal Al-Khwarizmi Engineering Journal, Vol. 9, No. 1, P.P. 39-46 (2013) An Autocorrelative Approach for EMG Time-Frequency Analysis Mohannad K. Sabir* Noor K. Muhsin** *,**Department of Biomedical Engineering/ Al-Khwarizmi College of Engineering/ University of Baghdad mohannad3151962@gmail.com: lEmai* **Email: Noorbmemsc81@yahoo.com (Received 17 March 2011; accepted 2 December 2012) Abstract As they are the smallest functional parts of the muscle, motor units (MUs) are considered as the basic building blocks of the neuromuscular system. Monitoring MU recruitment, de-recruitment, and firing rate (by either invasive or surface techniques) leads to the understanding of motor control strategies and of their pathological alterations. EMG signal decomposition is the process of identification and classification of individual motor unit action potentials (MUAPs) in the interference pattern detected with either intramuscular or surface electrodes. Signal processing techniques were used in EMG signal decomposition to understand fundamental and physiological issues. Many techniques have been developed to decompose intramuscularly detected signals with various degrees of automation. This paper investigates the application of autocorrelation function (ACF) method to decompose EMG signals to their frequency components. It was found that using the proposed method gives a quite good frequency resolution as compared to that resulting from using short time fast Fourier transform (STFFT); thus more MU’s can be distinguished. Keywords: Digital Signal processing; Nonstationary Processing of Biomedical Signals; Time-Frequency Analysis of Biosignals. 1. Introduction The electromyogram (EMG) is an electrical manifestation of the contracting muscle. This can be either a voluntary or involuntary muscle contraction. The EMG signal is a complicated signal which is affected by anatomical and physiological properties of the muscle, the control scheme of peripheral nervous system, as well as the instrumentation used for detection and the process used to record the EMG signals. The basic functional unit of the muscle contraction is a motor unit. This muscle fiber contracts when the action potentials of the motor nerve which supplies it reach a depolarization threshold. The depolarization generates an electromagnetic field measured as a very small voltage called EMG [1]. EMG signal is a complicated signal in the body that cannot be differentiated easily by the physician because of its shape. 2. EMG Measurements Although action potentials from individual muscle fibers can be recorded under special conditions, it is the electrical activity of the entire muscle that is of primary interest. In this case, the signal is a summation of all the action potentials within the range of the electrodes, each weighted by its distance from the electrodes as shown in Figure 1. Fig. 1. Surface EMG Summed Activities of Many MU’s are Recorded. mailto:mohannad3151962@gmail.com mailto:Noorbmemsc81@yahoo.com Mohannad K. Sabir Al-Khwarizmi Engineering Journal, Vol. 9, No.3, P.P. 39-46 (2013) 40 Since the overall strength of muscular contraction depends on the number of fibers energized and the time of contraction, there is a correlation between the overall amount of EMG activity for the whole muscle and the strength of muscular contraction. The EMG potentials from a muscle or group of muscles produce a complicated waveform that varies in amplitude with the amount of muscular activity. Peak amplitudes vary from 25 μV to about 5 mV, depending on the location of the measuring electrodes with respect to the muscle and the activity of the muscle. A frequency response from about 5 Hz to well over 5000 Hz is required for faithful reproduction [3]. The EMG signals of the tested muscles (right pectoralis, right biceps and right triceps) are shown in Figure 2. Fig. 2. The Tested EMGs Signals. According to Figure 2, time domain cannot give us a lot of information about muscle components, hence it is necessary to transform the signal into a frequency domain to obtain more information regarding the muscle fiber. 3. EMG Signal Processing Raw EMG offers us valuable information in a particularly useless form. This information is useful only if it can be quantified. Various signal processing methods are applied to raw EMG to reach the accurate and actual EMG signal. Wavelet analysis is chosen so as to match the shape of the MUAP; the resulting WT yields the best possible energy localization in the time-scale plane. It was found out that WT is an alternative to other time- frequency representations with the advantage of being linear. Wavelet analysis have been used to match the shape of the MUAP, for a unipolar recorded signal and under certain hypotheses presented by Gabor in 1946. Ismail and Asfour came with a theory saying that, the most common methods used to determine the frequency spectrum of EMG are the fast and short term Fourier transforms (FFT and SFT); but they also concluded that a major drawback of these transformation methods is that they assume the signal to be stationary. However, EMG signals are non-stationary [3]. To overcome the above limitation of FFT method, we propose the Autocorrelation function method for EMG spectral analysis as explained in the following section. 4. Methodology Resolving EMG signal into its significant constituent MUAPTs requires the ability to detect the discharges (i.e., MUAPs) of the MUs that are significantly contributing to the composite signal, and to correctly associate each detected MUAP with the MU that generated it. EMG signal decomposition therefore involves two basic steps: detecting MUAPs and recognizing the detected MUAPs [4]. Correlation based processing is known to be comparatively robust against noise. The autocorrelation function method [6], is classified into this category, and may be one which provides the best performance in noisy environments [6]. Thus; autocorrelation spectral analysis, instead of FFT techniques, is used in this work to characterize the distributions of MUAPs of different MUs. Successful EMG signal decomposition can therefore involve the following steps: (1) EMG signal acquisition. (2) Slicing the incoming signal into sequence of time windows, where for each window the dc- off set (which represents the mean of the signal) is removed and then the resultant window is multiplied by hamming window. (3) Applying autocorrelation function (ACF) to each window, & (4) Resolution (frequency components) and classification of superimposed MUAPs by searching for the maximums of ACF. A block diagram of the proposed system is shown in Figure 3. 0 500 1000 1500 2000 2500 3000 3500 4000 -2 0 2 s am ple no. am p. (v ) EMG of pec torals musc le 0 500 1000 1500 2000 2500 3000 3500 4000 -0.5 0 0.5 s am ple no. am p. (v ) EM G of biceps muscle 0 500 1000 1500 2000 2500 3000 3500 4000 -1 0 1 s am ple no. am p. (v ) EM G of triceps m us cle Mohannad K. Sabir Al-Khwarizmi Engineering Journal, Vol. 9, No.3, P.P. 39-46 (2013) 41 Fig. 3. The Proposed System’s Block Diagram A. EMG Signal Acquisition The recorded EMG signal coming from the electrodes is amplified by that little EMG box using a gain of 5,000 according to the manufacturer [7]. Data files that we used are downloaded from [8] with a sampling frequency of 2000 Hz and 4000 sample length (2 sec. time duration). The EMG amplifiers inside the box have different biases for the pectorals, biceps and triceps muscles in 5 lb weight. Analysis of these data files started with converting them to (.mat) files then the small DC offset of the signal is removed by forcing the mean of the window to be equal to zero, as shown in Figure 4. x1(t) = x0(t) - dc_offset …(1) Where x1(t) is the EMG signal without a DC offset, and x0(t) is the recorded signal. Fig. 4. EMG’s after Removing DC Offset. B. Hamming Window Multiplication Multiplying the truncated input signal by a hamming window is to reduce the abruptness of the truncated ends and thereby improve the frequency response. The hamming window is given by: )/2cos(46.054.0)( Niiw π−= …(2) for i=0,1,2,…N-1, and N is the length of the window [5]. x(i)= x1(i) * w(i) …(3) for i=0,1,2 … N-1. Since the autocorrelation method required at least two periods of the lowest frequency component (20 Hz) to give appropriate estimation of the contained frequencies, then a window of 128msec is chosen. . Moreover, it is required to maintain the stationary property of the EMG signal during each window. Fig (5) illustrates the mean of each window. According to the sampling time of 500µsec (Ts=1/fs), each successive window will represent 256 samples (128msec/500µsec). Thus, there will be 15 windows for the input EMG signal of 2 sec time duration. Number of windows = 4000(samples)/256= 15.625 Fig. 5. The Mean of Each Window of 256 Samples. C. Autocorrelation Method 1) Principle Autocorrelation function ( )kR is calculated by ( ) ∑ − = += 1 0 )()( 1 N n knxnx N kR …(4) 0 500 1000 1500 2000 2500 3000 3500 4000 -2 0 2 DC offset=.0608 sample no. am p. (v ) 0 500 1000 1500 2000 2500 3000 3500 4000 -0.5 0 0.5 DC offset=0.0619 sample no. am p. (v ) 0 500 1000 1500 2000 2500 3000 3500 4000 -1 0 1 DC offset=0.0597 sample no. am p. (v ) EMG signal Acquisition Hamming Window Autocorrelation Resolution & Classification Mohannad K. Sabir Al-Khwarizmi Engineering Journal, Vol. 9, No.3, P.P. 39-46 (2013) 42 where x(n) is the EMG signal, k is the lag number, and n is the discrete time index. The autocorrelation function ( )kR has a large value when )(nx is similar to )( knx + . If )(nx has a period of P, then ( )kR has peaks at k = lP where l is an integer. Essentially, ( )0R gives the largest value among ( )kR , k= lP, for l=0, 1, 2…. The second value is given by ( )PR . Other ( )kR peaks of usually decrease as k increases [6] for signals containing a single frequency component. 2) Implementation The proposed algorithm uses the autocorrelation function to compute the EMG signal spectrum. The maximum power frequencies region of the EMG signal is to be from 20 Hz to 150 Hz, which corresponds to the region (6.67 to 50) msec of the autocorrelation function. Since, the sampling frequency of the input EMG signal is 2000Hz and each truncated window is of 256 samples, then that time region that is represented is between the samples 14 and 100 of the ACF. Thus, each peak in this region is estimation to a frequency contained in the input EMG signal. The bold region illustrated in Fig (6) represents the above interested region of the ACF of the 14th window of the pectoralis EMG signal. It is quite clear that all the major and minor peaks of the ACF represent estimation of the corresponding frequencies. The frequency is the reciprocal of multiplying the peak location by the sampling time. For example, the peak at location 98, as illustrated in Fig (6), is estimation to: f(98)= 1/(98*500µsec)= 20.4081 Hz This represents the detected MUAP of MU. So that, the autocorrelation maxima between these two locations represent the frequency components of the EMG signal. Fig. 6. The Interested Search Region. 5. Results and Discussion We applied the autocorrelation method to compute the frequency spectrum of the tested EMG signals that sampled by 2000 Hz (500 μsec.). As described previously the maximum power frequencies of EMG signals will be appear in the region from sample 14 to sample 100 of the ACF. Figures 7, 8, and 9, illustrate the estimated frequency spectrum of each window of the total 14 windows contained in the pictoralis, biceps, and triceps muscles respectively. Mohannad K. Sabir Al-Khwarizmi Engineering Journal, Vol. 9, No.3, P.P. 39-46 (2013) 43 Fig. 7. Frequency Spectrum of Each Window of the Pectoralis EMG Signal. Fig. 8. Frequency Spectrum of Each Window of the Biceps EMG Signal. Fig . 9. Frequency Spectrum of Each Window of the Triceps EMG Signal. Mohannad K. Sabir Al-Khwarizmi Engineering Journal, Vol. 9, No.3, P.P. 39-46 (2013) 44 The resolution of the frequency spectrum using ACF is not linear as illustrated in Fig. (10). It is quite clear from this figure that at low frequencies (samples 90-100) we will get a frequency accuracy of less than 2 Hz, which mean a very good resolution, as mentioned previously that almost the MUs have a low frequency values. Fig. 10. ACF Frequency Resolution Verses Sample’s Location. This resolution is about 0.2 Hz which allows us to get more details about the corresponding MUs. While the frequency resolution, for high frequency band (samples 14-20), is between 8.3 to 11 Hz, which gives less details about the corresponding MUs. However, it seems that this low frequency resolution is quite enough for the high frequency band, since the EMG signal contains less number of motors having frequencies between 100 Hz to 150 Hz. The proposed autocorrelation method of EMG spectrum analysis may introduce false peaks at the multiples of the period. These peaks will be considered as estimation to a frequency while it is either the second or the third multiple of the pre- estimated one. To overcome this drawback, Tetsuya [6] proposed a weighted autocorrelation approach to improve the extracted pitch accuracy. Where, the autocorrelation function is weighted by the inverse of the average magnitude difference function (AMDF). The AMDF is described by [6]: ψ(k)= ∑ | ( ) + )| …(5) The AMDF has the characteristic that when x(n) is similar with x(n+k), ψ(k) becomes small. Thus, if x(n) has a period of P, ψ(k) produces a deep notch at k= P. Figure (11) illustrates the AMDF of the same window shown in Figure (6). Fig. 11. The AMDF of the Same Window Illustrated in Fig.(6). As illustrated in Fig (11) that AMDF has minimum peaks at the same samples of the maximum peaks of the ACF illustrated in fig (6). The weighted ACF is given by [6]: η(k)=R(k)/ (ψ(k)+M) …(6) Where, M is a fixed number (M>0). This weighted ACF is applied to all the windows of the pictoralis, biceps, and triceps muscles. Figure (12) illustrates the weighted ACF of the same window shown in Figure (6). Fig. 12. The Weighted ACF for the Same Window Illustrated in Fig (6). Mohannad K. Sabir Al-Khwarizmi Engineering Journal, Vol. 9, No.3, P.P. 39-46 (2013) 45 Fig. 13. The Weighted ACF & ACF of the Same Window. Figure (13) illustrates that the shape of the weighted ACF does not differ from that of the non- weighted ACF. Thus, one can conclude that all the maxima of the ACF are estimation to the fundamental frequencies, and are not to any of their harmonics. This is also true for all other windows of the same muscle and the other muscles too. 6. Conclusions Since the EMG signals are non-stationary, so they need a time frequency analysis to characterize the distributions of MUAPs of different MUs. Correlation based processing is known to be comparatively robust against noise. The autocorrelation function method, is classified into this category, and may be one which provides the best performance in noisy environments. Thus; autocorrelation spectral analysis is applied to short period Hamming windowed EMG signals of different muscles. The proposed method of spectral analysis of EMG signals using autocorrelation function (ACF), gives a good estimation of frequency components for non stationary signals such as EMG. Moreover, its nonlinear frequency resolution property (having high frequency resolution at low frequency band) is quite convenient for EMG signal, since more MUs have low frequency values, thus requiring high frequency resolution at this frequency band (20Hz-50Hz), to distinguish more MUs of the EMG signal. A weighted ACF is also applied to reduce the effects of harmonic peaks that may occur at the ACF. Although the tested EMG signals are few, the results are optimistic, and more investigation is recommended. 7. References [1] Motion Lab System, Inc., "A Software User Guide for EMG Graphing and EMG Analysis", 2009. [2] Jessica Zarndt, "The Muscle Physiology of Electromyography", http://www.holycross.edu/departments/biol ogy/kprestwi/phys'02/labs/emg_lab. [3] M.B.I. Raez, M.S. Hussain, and F. Mohd- Yasin, " Techniques of EMG signal analysis: detection, processing, classification and applications", Faculty of Engineering, Multimedia University. 63100 Cyberjaya, Selangor. Malaysia. 2006. [4] Roerto Merletti and Philips A. Parker, "Electromyography Physiology, Engineering and Noninvasive Applications",2004. [5] Steven W. Smith, "The Scientist and Engineer's Guide to Digital Signal Processing", 2009. [6] Tetsuya Shimamura, Hajime Kobayashi, "Weighted Autocorrelation for Pitch Extraction of Noisy Speech", IEEE transaction on speech and audio processing, vol.9, no.7, October 2001. [7] BioResearch Inc. [8] http://www.BMES642:EMG analysis. http://www.holycross.edu/departments/biol http://www.BMES642:EMG )2013( 39-46 ، صفحة1دد، الع9 مجلة الخوارزمي الھندسیة المجلد ابر مھند كاظم ص 46 تحلیل الطیف الترددي إلشارة تخطیط العضالت باستخدام دالة التطابق الذاتي **نور كمال محسن *مھند كاظم صابر جامعة بغداد/ كلیة الھندسة الخوارزمي/ تيقسم ھندسة الطب الحیا**،* mohannad3151962@gmail.com *االلكتروني البرید: Noorbmemsc81@yahoo.com **االلكتروني البرید: الخالصة EMGتحلیإلشارة . و ذلك یؤدي إلى فھم استراتیجیات التحكم في العضالت والحاالت المرضیة. اللبنات األساسیة للنظام العصبي العضلي MUتعتبر . العضلي أو السطحيMUلھا الدور في الكشف عن نمط التداخل مع ) MUAPs(ھي عملیة تحدید وتصنیف ، كما ان وحدة العمل و اإلمكانات الحركیة وقد تم تطویر العدید من التقنیات لتتحلل اشارات الكشف . لفھم القضایا األساسیة ، والفسیولوجیة EMGت تقنیات معالجة اإلشارات في التحلل إشارة واستخدم EMGات إشارات و إیجاد الطیف الترددي إلى مكون) ACF(ھذا البحث تدارس في تطبیق دالة التطابق الذاتي . عضلیا مع درجات مختلفة من التشغیل اآللي ، وبالتالي )STFFT(فقد وجد أن استخدام الطریقة المقترحة یعطي القرار لطریقة التردد بانھا جیدة جدا بالمقارنة مع ما ھو ناتج عن استخدام . و تكرارھا .MUیمكن تمییز أكثر من mailto:mohannad3151962@gmail.com mailto:Noorbmemsc81@yahoo.com