{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Heart rate estimation, signal de-noising and mobile health\n", "The goal of this lab is to study an algorithm for correctly retrieving the timing instants of an ECG (EKG) time series. We will implement a simplified variant of the [Pan-Tompkins algorithm](https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4122029) and we will test it on real ECG data recordings. Towards the end of the lab you will be asked to apply your algorithms to the problem of heart rate estimation using a mobile device.\n", "\n", "### Objectives\n", "By the end of this session you should be able to understand the main challenges in the estimation of EGC signals. How to work from modeling assupmtions in order to exploit your prior knowledge in the estimation task.\n", "\n", "* _Note: Everything below this line should be truned in as part of your lab report._\n", "\n", "***\n", "\n", "### Signal generation and template matching\n", "In this task we will implement a simple version of a QRS-complex detector by using template matching, thresholding and non-maximum suppresion on the\n", "observed ECG signal. We will then use those detections to provide an estimate of the heart rate.\n", "\n", "* **Signal Generation -** We are going to create a synthetic EGC signal that we will be using to test our detector algorithm. For that purpose generate a stream of $10$ equally spaced pulses over a time span of $10\\ s$ and with a sampling rate of $f_s = 256\\ Hz$. The goal of this task is to create a basic model for a ECG signal. You can think of the signal as the convolution of a canonical pulse shape $\\varphi(t)$ with a stream of Dirac delta functions: \n", "\n", " $$x(t) = \\sum_{k=1}^K a_k\\, \\varphi(t-t_k) = \\varphi(t) \\ast \\sum_{k=1}^K a_k\\, \\delta(t-t_k),$$ \n", "\n", " where $K$ is the number of pulses observed, $a_k$ represent the amplitudes of the different pulses and $t_k$ correspond to the locations of the pulses in time. For the purpose of this task use the canonical pulse shape generated by the `ecg_wave()` function defined in the notebook and constant unit amplitudes $a_k=1$. Plot the generated signal over time adding the appropriate labels for the axis.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ============================================================================\n", "# import modules\n", "# ============================================================================\n", "# Note that this part of the code needs to be run prior to any other code cell\n", "\n", "import numpy as np\n", "import json\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "from scipy import signal\n", "\n", "# inline plots\n", "%matplotlib inline\n", "sns.set()\n", "\n", "# ============================================================================\n", "# functions\n", "# ============================================================================\n", "def ecg_wave(x):\n", " \"\"\"\n", " This function generates a synthetic ECG template of unit length (support).\n", "\n", " \"\"\"\n", " \n", " # compute signal (superposition of splines)\n", " return 0.3 * signal.bspline(3*3*x-7.5,2)\\\n", " + 0.15* signal.bspline(3*4*x-2,3)\\\n", " + signal.bspline(3*2*2*x-6,1) - 0.2 * signal.bspline(3*4*x-5,1) - 0.4 * signal.bspline(3*4*x-7,1)\n", "\n", "\n", "# ============================================================================\n", "# Start your code below this line\n", "# ============================================================================\n", "# ...\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Simulate Noise -** Generate a noisy version of the synthetic ECG signal generated before by adding Gaussian noise with standard deviation $\\sigma=0.5$. Plot the noisy observations. Can you identify the locations of the QRS-complex?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ...\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "* **Find Peaks -** Implement a QRS-complex detector on the noisy signal you just generated. For that purpose we will use template matching, thresholding, and non-maximum suppression. We recommend that you normalize the signal before thresholding. For instance, you can normalize the signal to take values in the range $[0,1]$ or you could standardize the signal by making it zero mean and unit variance. After normalization define a threshold value and keep only those values of the signal that are above the given threshold. Once you have done that, perform a non-maximum suppression (i.e., keep a value if it is greater than the previous and following values). To implement the thresholding and non-maximum suppression operations you can use `signal.find_peaks` function." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def normalize_range(x):\n", " \"\"\"\n", " Normalizes the signal such that its minimum value is 0 and its maximum value is 1,\n", " provided the signal is not identically zero.\n", " \"\"\"\n", " # check that there are non-zero elements\n", " if np.any(x):\n", " \n", " # subtract minimum\n", " minx = np.min(x)\n", " z = x - minx\n", " \n", " # divide by max value\n", " maxz = np.max(z)\n", " return z/maxz\n", " \n", " else:\n", " \n", " return x\n", "\n", "# ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **RR Intervals -** Plot the original (synthetic) ECG signal as well as the locations of the peaks to verify your method. Keep in mind to compensate for any delay you might have introduced by filtering. From the peak binary signal, estimate the $RR$ interval sequence $r_n$ and its average value as:\n", "\n", "$$ \\bar R = \\frac{1}{N}\\sum_{n=0}^N r_n, $$\n", "\n", "where $N$ is the number of peaks detected. Provide an estimate of the average heart beat rate in beats per minute $[bpm]$.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ...\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### QRS complex detection via Pan-Tompkins algorithm\n", "Now let us consider a more realistic scenario where we have a noisy ECG signal where we don't know a priori the shape of the QRS-complex waeform. There are different\n", "sources of noise that can be present during an ECG acquisition. In addition to high-frequency thermal\n", "noise, it is also common to observe the presence of low-frequency interference coming from breathing. \n", "When it comes to pathologies, different non-additive distortions might be present on the ECG signal \n", "that alter the shape of the QRS-complex itself but for the purpose of this task we just assume a healthy \n", "individual where signal distortion comes solely from the acquisition process. The procedure that we\n", "will employ in order to estimate the locations of the QRS-complex is based on the Pan-Tompkins\n", "algorithm [Pan-Tompkins algorithm](https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4122029).\n", "\n", "* **Effect of noise** - Load the signal **ecg_mitnst.json** and plot the signal over time. The signal corresponds to a sample ECG from the [MIT noise stress dataset](https://physionet.org/physiobank/database/nstdb/)\n", "that has been downloaded from [physionet.org](https://physionet.org/physiobank/database/).\n", "Plot the signal and observe the presence of a strong low-frequency component." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# access files stored on Google Drive from Colab\n", "from google.colab import drive\n", "drive.mount('gdrive/')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# load data from Google Drive\n", "with open('gdrive/My Drive/bmdslab/lab-01/ecg_mitnst.json','r') as infile:\n", " data = json.load(infile)\n", "\n", " # print data to see dictionary fields\n", " print(data)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Pan-Tompkins Algorithm -** In order to deal with noise we will implement a simplified version of the Pan-Tompkins algorithm. For that purpose proceed as follows:\n", "\n", " 1. _Bandpass filtering:_ Follow the steps described in the manual in order to perform a bandpass filtering operation as the concatenation of lowpass and highpass filtering operations. Display the obtained signals at the intermediate steps.\n", " \n", " 2. _Differentiation:_ Use the discrete filter $d_n$ to approximate the derivative of the underlying signal:\n", " \n", " $$ d_n = \\frac{1}{8}\\big(\\delta_{n-2} + 2\\delta_{n-1} -2\\delta_{n+1} - \\delta_{n+2}\\big) $$\n", " \n", " The filter is intended to localize the steepest region in the QRS-complex.\n", " \n", " 3. _Signal squaring:_ Square the obtained signal after differentiation and plot the obtained waveform. Why was squaring helpful in revealing the peaks of the QRS complex?\n", " \n", " 3. _Integration:_ Integrate the resulting signal from the squaring operation with a box window of length $L=50$. Display the resulting signal.\n", " \n", " 4. _Peak detection:_ Use a peak detector to reveal the locations of the peaks. Estimate the average $RR$ interval and plot the sequence of estimated peaks.\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mobile health\n", "In this last task we will consider the estimation of the heart rate in a less controlled scenario. In this case,\n", "the signal comes from a mobile device and it is not a one-dimensional time-series but a video recording. Those with a\n", "smartphone can record and use their own signals for the experiment. Alternatively, you can use the video\n", "provided for that purpose. The recorded signal consists of a video of the tip of the finger placed\n", "right in front of the camera while the camera flash is on. Play this video on your notebook and try to observe\n", "the intensity variations due to the heart beat.\n", "\n", "* Using the video signal and the methods studied in previous sections estimate the average heart rate from it.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import io\n", "import base64\n", "from IPython.display import HTML\n", "\n", "# this code displays the video\n", "video = io.open('gdrive/My Drive/bmdslab/lab-01/ppg_resized.mp4', 'r+b').read()\n", "encoded = base64.b64encode(video)\n", "HTML(data=''''''.format(encoded.decode('ascii')))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import cv2\n", "\n", "# create video reader object\n", "vread = cv2.VideoCapture(\"gdrive/My Drive/bmdslab/lab-01/ppg_resized.mp4\")\n", "\n", "# get number of frames\n", "nframes = int(vread.get(cv2.CAP_PROP_FRAME_COUNT))\n", "\n", "# get frame rate\n", "fps = vread.get(cv2.CAP_PROP_FPS)\n", "\n", "# loop over frames\n", "for ff in range(nframes):\n", " # read frames\n", " ret, img = vread.read()\n", " \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ..." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }