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 time and 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.
Figure 1: Unfiltered vs filtered EMG signals with different low pass cut-off frequencies
Summary
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.
Scientifically Sound 
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.
LikeLike
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
LikeLike
Dear Joanna,
Thank you for creating this code. This is really helpful. I have one question, what is the difference between low_pass and low_band?
LikeLike
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 whilelow_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.
LikeLike
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?
LikeLike
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.
LikeLike
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.
LikeLike
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!!!
LikeLike
Thanks for your feedback, glad the info helps
LikeLike
Dear Joanna:
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.
LikeLike
Hi there,
Thanks for picking that up. Code has been updated.
LikeLike