Week 6: Jupyter Notebook – Calculating Live BPM

The next challenge for me was to find a way to calculate live BPM from my new ECG data, and then publish this to be visualised elsewhere. After feeling lost while researching for methods in C++, another summer studio facilitator Trang Nguyen recommended for me to instead look at using Python. Although I was unable to find her work, I was able to find a couple of examples using Python to calculate BPM from ECG data.

The example below used a csv to analyse a discrete amount of ECG data.

As I was still very unfamiliar with how to use Python, the process of understanding the example code and how to incorporate live data was quite challenging for me and I needed to ask for frequent additional help. Without understanding the fundamentals, it is hard to know where to even begin, but luckily I was able to access help and advice from both my studio team and our facilitator.

To still utilise this example for its peak detection and BPM calculation, I learnt that I should used a pandas data frame to store a set amount of my own ECG data to then be analysed. To do this, following some advice, I created an initial empty dataframe called “ECGDataSet”, then within the on_message function, created a new dataframe “newData” to be concatenated with the initial “ECGDataSet”. Also, I learnt that ECGDataSet needed to be a global variable as before it wasn’t concatenating and only recognising the initial empty data frame. I could tell this was occurring by adding the output “ECGDataSet.tail()” or “ECGDataSet.head()”. The data frame was also limited to 2500 inputs by adding “ECGDataSet.tail(2500)” to the concatenate line. To make the BPM calculation update automatically, our facilitator helped me set up a loop with a sleep time of 0.1s, meaning the plot will update every 0.1s and recalculate the average BPM.

# !pip install "notebook>=5.3" "ipywidgets>=7.2"
# !pip install paho-mqtt
# Get the pandas dataset ready for incoming data
import pandas as pd
ECGDataSet = pd.DataFrame([0], columns=['ECG'])
import paho.mqtt.client as mqtt

# The callback for when the client receives a CONNACK response from the server.
def on_connect(client, userdata, flags, rc):
    print("Connected with result code "+str(rc))
    # Subscribing in on_connect() means that if we lose the connection and
    # reconnect then subscriptions will be renewed.
    client.subscribe("ECG")
    
# The callback for when a PUBLISH message is received from the server.
def on_message(client, userdata, msg):
    payloadStr = str(msg.payload).split("'")[1]
    payloadInt = float(payloadStr)
    #print(payloadInt)
    newData = pd.DataFrame([payloadInt], columns=['ECG'])
    global ECGDataSet
    ECGDataSet = pd.concat([ECGDataSet.tail(2500), newData], ignore_index=True)

client = mqtt.Client()
client.on_connect = on_connect
client.on_message = on_message
client.username_pw_set("username", password="password")
client.connect("139.180.169.188", 5269, 60)
client.loop_start()
import time
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math
from IPython import display

# Setup a loop with the ability to interrupt
try:
    while True:
        # begin analyzing the dataset for peaks
        dataset = pd.DataFrame()
        dataset['hart'] = ECGDataSet['ECG']
        #Calculate moving average with 0.75s in both directions, then append do dataset
        hrw = 0.75 #One-sided window size, as proportion of the sampling frequency
        fs = 100 #The example dataset was recorded at 100Hz
        mov_avg = dataset['hart'].rolling(int(hrw*fs)).mean() #Calculate moving average
        #Impute where moving average function returns NaN, which is the beginning of the signal where x hrw
        avg_hr = (np.mean(dataset.hart))
        mov_avg = [avg_hr if math.isnan(x) else x for x in mov_avg]
        mov_avg = [x*1.2 for x in mov_avg] #For now we raise the average by 20% to prevent the secondary heart contraction from interfering, in part 2 we will do this dynamically
        dataset['hart_rollingmean'] = mov_avg #Append the moving average to the dataframe
        #Mark regions of interest
        window = []
        peaklist = []
        listpos = 0 #We use a counter to move over the different data columns
        for datapoint in dataset.hart:
            rollingmean = dataset.hart_rollingmean[listpos] #Get local mean
            if (datapoint < rollingmean) and (len(window) < 1): #If no detectable R-complex activity -> do nothing
                listpos += 1
            elif (datapoint > rollingmean): #If signal comes above local mean, mark ROI
                window.append(datapoint)
                listpos += 1
            else: #If signal drops below local mean -> determine highest point
                maximum = max(window)
                beatposition = listpos - len(window) + (window.index(max(window))) #Notate the position of the point on the X-axis
                peaklist.append(beatposition) #Add detected peak to list
                window = [] #Clear marked ROI
                listpos += 1
        ybeat = [dataset.hart[x] for x in peaklist] #Get the y-value of all peaks for plotting purposes
        plt.title("Detected peaks in signal")
        plt.xlim(0,2500)
        plt.plot(dataset.hart, alpha=0.5, color='blue') #Plot semi-transparent HR
        plt.plot(mov_avg, color ='green') #Plot moving average
        plt.scatter(peaklist, ybeat, color='red') #Plot detected peaks
        display.clear_output(wait=True)
        plt.show()
        
        
        # Calculate the average bpm
        RR_list = []
        cnt = 0
        while (cnt < (len(peaklist)-1)):
            RR_interval = (peaklist[cnt+1] - peaklist[cnt]) #Calculate distance between beats in # of samples
            ms_dist = ((RR_interval / fs) * 1000.0) #Convert sample distances to ms distances
            RR_list.append(ms_dist) #Append to list
            cnt += 1
        bpm = 60000 / np.mean(RR_list) #60000 ms (1 minute) / average R-R interval of signal
        print ("Average Heart Beat is: %.01f" %bpm) #Round off to 1 decimal and print
        
        client.publish("BPM", (bpm))

        # sleep
        time.sleep(0.1)
except KeyboardInterrupt:
    pass
Output plot showing detected peaks and average BPM value
New BPM data stored in Influxdb

To then publish my BPM value being calculated every 0.1 seconds, it only required the simple addition of “client.publish(“BPM”, (bpm))”, with the topic being BPM and the message bpm. The Node-RED function node phrasing didn’t need any alterations and I could begin sending BPM data into our Influx database to then be visualised on Grafana once again.

Grafana dashboard with new BPM data
laser cutting the lid
Autodesk Fusion 630 sketch of lid

The final step to complete by project was to laser cut a clear acrylic lid for my enclosure. Although we went able to operate the laser cutter ourselves as it requires an induction, we were able to watch our facilitator set up each type of operation for the machine. Firstly, the outer 2mm perimeter was ‘vector rastered’ down by 1mm to create an indented lip to fit perfectly into the enclosure. The text was then rastered as well, however the text first needed to be flipped as the rastering was being executed onto the underneath side of the lid. The final stage was to ‘vector cut’ the outer perimeter and the hole for the switch. Upon putting the whole device and enclosure together, I slightly mis-measured where the hole needed to be as it was a bit off center and quite hard to turn the switch on and off without either long nails or a pen. If I were to make an enclosure again I would definitely do a test laser cut on plywood before jumping straight into the acrylic to correct any errors.

S U M M A R Y

Accomplished

– Calculated BPM using Python and visualised this data in Grafana
– Produced a clear acrylic lid for enclosure


Learnt

– Python: data frames, functions, global variables and loops
– Watched how to operate laser cutter machine


What next

– Continue to learn programming in C++ and Python
– Reflect on everything I have learnt in the past 6 weeks
– Begin to think about the next IoT project I can start to work on


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s