Python: Analysing EMG signals – Part 4
We have seen how Python can be used to process and analyse EMG signals in lessons 1, 2 and 3. When EMG signals are filtered, how does changing filter settings change the appearance of the filtered EMG signal?
A low pass filter allows frequencies below the cut-off frequency to pass through (ie. higher frequencies are removed). The filtered EMG signal can be used to interpret different physiological properties. For example, scientists investigating muscle force and muscle activity often use a low pass filter to capture the shape or “envelope” of the EMG signal as this is thought to better reflect force generated by a muscle. Let’s see how changing the low pass cut-off frequency changes the rectified signal.
We want to write a bit of code to filter and rectify an EMG signal and plot graphs comparing the unfiltered and filtered signal, for different low pass cut-off frequencies. Since this involves doing things repeatedly (eg. EMG processing, plotting) for different values (low pass cut-offs), we will write a function to do the EMG processing etc. and send the function different values in a
for loop. For quick refreshers, see this series on functions and this post on loops in Python.
Here is our function and the
for loop at the end. Key bits of code are discussed below:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
def filteremg(time, emg, low_pass=10, sfreq=1000, high_band=20, low_band=450): """ time: Time data emg: EMG data high: high-pass cut off frequency low: low-pass cut off frequency sfreq: sampling frequency """ # normalise cut-off frequencies to sampling frequency high_band = high_band/(sfreq/2) low_band = low_band/(sfreq/2) # create bandpass filter for EMG b1, a1 = sp.signal.butter(4, [high_band,low_band], btype='bandpass') # process EMG signal: filter EMG emg_filtered = sp.signal.filtfilt(b1, a1, emg) # process EMG signal: rectify emg_rectified = abs(emg_filtered) # create lowpass filter and apply to rectified signal to get EMG envelope low_pass = low_pass/(sfreq/2) b2, a2 = sp.signal.butter(4, low_pass, btype='lowpass') emg_envelope = sp.signal.filtfilt(b2, a2, emg_rectified) # plot graphs fig = plt.figure() plt.subplot(1, 4, 1) plt.subplot(1, 4, 1).set_title('Unfiltered,' + '\n' + 'unrectified EMG') plt.plot(time, emg) plt.locator_params(axis='x', nbins=4) plt.locator_params(axis='y', nbins=4) plt.ylim(-1.5, 1.5) plt.xlabel('Time (sec)') plt.ylabel('EMG (a.u.)') plt.subplot(1, 4, 2) plt.subplot(1, 4, 2).set_title('Filtered,' + '\n' + 'rectified EMG: ' + str(int(high_band*sfreq)) + '-' + str(int(low_band*sfreq)) + 'Hz') plt.plot(time, emg_rectified) plt.locator_params(axis='x', nbins=4) plt.locator_params(axis='y', nbins=4) plt.ylim(-1.5, 1.5) plt.plot([0.9, 1.0], [1.0, 1.0], 'r-', lw=5) plt.xlabel('Time (sec)') plt.subplot(1, 4, 3) plt.subplot(1, 4, 3).set_title('Filtered, rectified ' + '\n' + 'EMG envelope: ' + str(int(low_pass*sfreq)) + ' Hz') plt.plot(time, emg_envelope) plt.locator_params(axis='x', nbins=4) plt.locator_params(axis='y', nbins=4) plt.ylim(-1.5, 1.5) plt.plot([0.9, 1.0], [1.0, 1.0], 'r-', lw=5) plt.xlabel('Time (sec)') plt.subplot(1, 4, 4) plt.subplot(1, 4, 4).set_title('Focussed region') plt.plot(time[int(0.9*1000):int(1.0*1000)], emg_envelope[int(0.9*1000):int(1.0*1000)]) plt.locator_params(axis='x', nbins=4) plt.locator_params(axis='y', nbins=4) plt.xlim(0.9, 1.0) plt.ylim(-1.5, 1.5) plt.xlabel('Time (sec)') fig_name = 'fig_' + str(int(low_pass*sfreq)) + '.png' fig.set_size_inches(w=11,h=7) fig.savefig(fig_name) # show what different low pass filter cut-offs do for i in [3, 10, 40]: filteremg(time, emg_correctmean, low_pass=i)
Line 1. Define a function called
filteremg to accept
emg values (for plotting on the x- and y-axes) with default values for a low pass filter (for the EMG envelope), sampling frequency, and high and low pass filter frequencies. The default values mean that when we call the function, if we don’t specify different values, the function uses the defaults during code execution.
Line 28 onwards. Plot 3 subplots to see (1) the unfiltered, unrectified EMG signal, (2) the filtered, rectified signal, (3) the rectified signal with a low pass filter to get the EMG envelope and (4) a zoomed-in section of the signal from (3) over the time period indicated by the red line to see the underlying shape of the final signal. This line says to index values from 0.9 to 1.0 seconds, convert these numerical float values to integers, and plot time on the x- and EMG signal on the y-axis.
Line 70. Write a
for loop that says: for each low pass cut-off value in the list
[3, 10, 40], when
i is 3, call the function
filteremg and assign the value of
i to the local variable
low. When i is 10, call the function and do the same. Etc.
For each of the cut-off frequencies supplied, we see in Figure 1 below that low cut-off frequencies produce a “smoothed” EMG signal which has a more noticable shape (envelope) compared to when higher cut-off frequencies are used.
There is still discussion around which lowpass cut-off frequencies are appropriate to get the EMG envelope under different conditions (eg. sustained contractions, gait). The biomechanist David Winter argues that choosing a lowpass cut-off value should be based on how the biomechanics of production of muscle force is modelled; see here for a discussion on the linear EMG envelope.
Changing the cut-off frequencies changes the shape and appearance of the filtered EMG signal. We now have a function that can be used to see how changing these frequencies changes the appearance of the signal. As a learning activity, stick the function into an IPython notebook (or similar) and play around with the filter settings to see just how much the shape can change.
Thank you very much for your excellent knowledgable series. I search and search two days finally find your site. Very helpful to do my project work. Once again Thanks a lot.
Since this is an EMG series, for reference, Martin Héroux wrote a good Python package to import EMG and transducer signals from Spike2 files. The package performs some basic signal processing, including EMG filtering and rectification. See:
Scientifically Sound blogpost: https://scientificallysound.org/2020/12/07/spike2py/
Stable package on PyPI: https://pypi.org/project/spike2py/
Development package on Github: https://github.com/MartinHeroux/spike2py
Thank you for creating this code. This is really helpful. I have one question, what is the difference between low_pass and low_band?
Hi there, thanks for your question.
I wanted to create two filters in that function: a low pass filter, and a band pass filter. The low pass filter should only let through frequencies 10 Hz and lower (from the default cutoff) while the bandpass filter should let through frequencies from 20 to 450 Hz. I.e. the low pass cutoff for the bandpass filter is 450 Hz. So I created the function to receive two arguments for the different low pass frequency cutoffs.
low_passis the cutoff for the low pass filter while
low_bandis the low pass cutoff for the bandpass filter.
In retrospect, it would be better (more Pythonic!) to create separate functions for each of the low pass and band pass filters.
Thank you so much for clarification. I have another question, I want to do RMS signal averaging. Can you let me know what function I can use for this. Do I need to do RMS of unfiltered signal and follow rest of code similarly as you described?
Hi Joanna, thank you for the clarification. Adding to my previous question, I want to do RMS processing of the signal. My data is in pandas dataframe and not sure how to best do it.
You would need to write a function for RMS EMG. Briefly: define a window length (e.g. 50 ms) and sample EMG values. Within the window, square the values, calculate the mean, then square root the mean. Move the window along the signal and do the same. At the start and end of the signal, shrink the window length because you will run out of samples as the window gets closer to the edges.
Dear Joanna, you don’t know how glad i’m that i found this web page. The work that you’ve done (as a team) is incredible. Thanks you so much for this tutorial. I’m making my bachelor project and it has really helped me a lot!!!
Thanks for your feedback, glad the info helps
Excited to see your research report! There is a little question on your code in Line 24, just wonder why you did not divide two? because you did it during bandpass procedure :p. Looking forward to your update.
Thanks for picking that up. Code has been updated.