Sometimes when we plot data we want to indicate a relevant feature or have a user select a section of the data. For example, if we record electromyography (EMG) and force signals, we might want to have the user mark the baseline EMG level so that it can be removed, and determine the amount of EMG produced during a maximum voluntary contraction (MVC).

This tutorial shows how to mark regions of interest on plots using Python’s Matplotlib library. We begin by reading in real EMG and force data (collected during two MVCs) saved in a text file, then use a number of functions to process the data. Many of these functions have been covered in previous tutorials, so working through this tutorial will help you revise and consolidate material from before. Links to previous lessons are provided where relevant. If you understand the commands, follow the way they are structured and can begin to write them from memory, you are doing very well.

The full set of commands is given below. The commands (saved in the Python file mvc.py) and the dataset (mvc.lvm) are available for download here.

 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

import numpy as np
import scipy as sp
from scipy import signal
import matplotlib.pyplot as plt

# Read in force and EMG data
infile = open('mvc.lvm', 'r')
line = infile.readlines()[23:]
infile.close()

data = [row.strip().split(',') for row in line]

time, force, emg = [], [], []
time = np.array([float(row[0]) for row in data])
force = np.array([float(row[1]) for row in data])
emg = np.array([float(row[2]) for row in data])

# Write a function to process EMG, and process the EMG signal
sampling_rate = 2000
def getEMG(emg, sampling_rate, highpass=20, lowpass=450):
    # remove mean
    emg = emg - np.mean(emg)
    # create filter parameters and digitally filter signal
    high, low = highpass/sampling_rate, lowpass/sampling_rate
    b, a = sp.signal.butter(4, [high, low], btype='bandpass')
    emg = sp.signal.filtfilt(b, a, emg)
    # take absolute of emg
    emg = abs(emg)
    return emg

emg = getEMG(emg, 2000)

# Find MVC force, mark corresponding EMG and 50 ms window across it on graph
force_mvc = np.max(force)
index_mvc = np.where(force == force_mvc)[0]
index_mvc = int(index_mvc)
# find no. of samples over 25 ms
width = int(0.025 * sampling_rate)
# find mean EMG and highest possible EMG voltage
emg_mvc = np.mean(emg[index_mvc-width : index_mvc+width])
emg_max = np.max(emg)
print('MVC EMG over 50 ms (a.u.): {:.2f}'.format(emg_mvc))

Line 1-4. Import necessary libraries.

Line 6-16. Read in force and EMG data from the text file mvc.lvm using the readlines() function. The file comes with a header followed by time, force and EMG data separated by commas; so, only data from line 23 onwards are read in. Data are stored as string values. For each row, using a list comprehension, strip the new line characters from the end and split the values by commas. Create lists to contain the data. Use more list comprehensions to read in string values from each row, convert them to floating point numbers, and finally convert the lists into Numpy arrays. The series on working with text files explains how to read data stored as text into Python.

Line 18-31. Write a function getEMG() that accepts arguments specifying the dataset, sampling rate and filter cut offs to process the EMG signal. Assign default values to highpass and lowpass arguments. Process the EMG signal. The series on functions and how to process EMG signals explains how to write functions in Python, and the background and techniques to process EMG signals.

Line 33-42. Find the force at MVC and use this to find the corresponding index (ie. sample number) value. We want to use force at MVC to find the corresponding EMG and calculate the mean EMG over a 50 ms window there. The width of half the window will be the number of samples over 25 ms. Use these values to calculate the mean EMG over 50 ms, and also find the highest possible EMG value. Print out the EMG at MVC to two decimal places. The series on Python’s print() statement explains the print functionality.

 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

# Have user select baseline EMG region
plt.plot(time, emg)
x = plt.ginput(2)
print('Start and stop x-values: {:.2f}, {:.2f}'.format(x[0][0], x[1][0]))
plt.show()
plt.close()

# Find index values for start and stop points
start_time = x[0][0]
stop_time = x[1][0]
for i,j in enumerate(time):
    if j == np.around(start_time, decimals=2): # round to 2 dec pla
        start_index = i
    if j == np.around(stop_time, decimals=2):
        stop_index = i
print('Start and stop index values: {:d}, {:d}'.format(start_index, stop_index))

# Plot force and EMG on time
fig = plt.figure()
plt.clf()
plt.subplot(2,1,1)
plt.plot(time, force)
plt.xlabel('Time')
plt.ylabel('Force (V)')
plt.subplot(2, 1, 2)
plt.plot(time, emg)
plt.xlabel('Time')
plt.ylabel('EMG (a.u.)')

# Mark baseline and MVC regions on plot
x1, x2 = time[index_mvc-width], time[index_mvc+width]
y1, y2 = emg_max+0.01, emg_max+0.01
plt.plot([x1, x2], [y1, y2], color='r', linestyle='-', linewidth=5,)
plt.plot(time[start_index:stop_index], emg[start_index:stop_index], color='g')

# Format figure spacing and save
fig.tight_layout()
fig.savefig('mvc.png')

Line 1-6. Plot the EMG data. When the graph pops up, the user is to click two points on the graph to select the region between them. This is accomplished using the ginput() function. Eg. the user is to select start and stop points of the EMG signal when the subject was relaxed. Print out the coordinates of these points (in time and EMG values) to two decimal places.

Line 8-16. Use the time and EMG coordinates to find the corresponding index values using a for loop. The tutorial on for loops explains how these work. Print out the index values as integers.

Line 18-34. Plot subplots of force and EMG on time. We would like to mark with a red line the region over which MVC EMG was calculated. The x and y coordinates of the start and end of the line are assigned to x1, x2, y1, and y2. These are passed to plt.plot() and the features of the line are specified. Also, index the time and EMG data that were selected by the user and plot the indexed values in green.

Line 36-38. Finally, fix the relative spacing of the subplots and save the figure.

After all this work, we produce the following output and Figure 1.

1
2
3

MVC EMG over 50 ms (a.u.): 0.90
Start and stop x-values: 0.13, 0.54
Start and stop index values: 260, 1080

Of course, remember to version control your code with Git.

Figure 1:

 

Summary

We learned how to mark regions of interest on plots using Matplotlib, and have a user select regions of interest using the ginput() function. Along the way, we revised material on reading in data from text files, functions, for loops, print statements, plotting, and analysing EMG data.