EMG_processing.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "EMG_processing.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyOq6fwtYNf6355YgNq9yypW",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/luuleitner/EMG_tutorial/blob/main/EMG_processing.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "<h1>EMG Tutorial</h1>\n",
        "\n",
        "This short tutorial is a hands-on introduction to electromyography (EMG) signal processing and analyses. During this exercise, we will explore a single differential EMG's amplitude and spectral features recorded on two biceps brachii contractions. This practical example is part of the Graz University of Technology lecture series *Applied Electrophysiology and Sensors*.\n",
        "\n",
        "I wish everyone a great dive into the topic, and please do not hesitate to <a href=\"mailto:christoph.leitner@tugraz.at\">contact</a> me in case of **any** questions!\n",
        "\n",
        "yours,<br>\n",
        "Christoph<br><br>\n",
        "\n",
        "<h4>Free EMG Ressources:</h4>\n",
        "\n",
        "*   <a href=\"https://delsys.com/support/documentation/#tutorial\">Practical Introduction into EMG</a> from Carlo de Luca (provided open-source by Delsys)\n",
        "*   <a href=\"https://www.robertomerletti.it/en/emg/material/teaching/\">Lecture slides</a> and other relevant EMG resources provided by Roberto Merletti\n",
        "*   Publication on a <a href=\"https://ieeexplore.ieee.org/document/9630034\">tattoo-based EMG system</a> [1] - one recording of this study is used in this tutorial.<br><br>\n",
        "\n",
        "\n",
        "<h4>Contact</h4>\n",
        "Christoph Leitner<br>\n",
        "e <a href=\"mailto:christoph.leitner@tugraz.at\">christoph.leitner@tugraz.at</a><br>\n",
        "g <a href=\"https://github.com/luuleitner/EMG_tutorial\">github.com/luuleitner/EMG_tutorial</a><br><br>\n",
        "\n",
        "---\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "YoR5UhYG_Puy"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<h2>Getting Started</h2>\n",
        "\n",
        "\n",
        "```\n",
        "[#] Areas shown like this are executable code. Use the mousover play button to run these cells.\n",
        "```\n",
        "\n",
        "\n",
        "First we need to install the EMG tutorial package from the the <a href=\"https://github.com/luuleitner/EMG_tutorial\">GitHub repository</a>."
      ],
      "metadata": {
        "id": "pXK6Gx5UOos3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install git+https://github.com/luuleitner/EMG_tutorial"
      ],
      "metadata": {
        "id": "-fUXe5OLU-1Z"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "We load all necessary libraries into our Colab notebook in this next step."
      ],
      "metadata": {
        "id": "vi9cZ8cDPj_T"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np # numpy scientific computing package\n",
        "import pandas as pd # pandas data analysis library\n",
        "import matplotlib.pyplot as plt # matplotlib library to plot graphs and data\n",
        "from scipy import signal # package providing fundamental algorithms for scientific computing\n",
        "\n",
        "from emg_processing.filter import Butter # Butterworth filter function from the EMG-tutorial package\n",
        "from emg_processing.analyses import rms # Root Mean Square function from the EMG-tutorial package"
      ],
      "metadata": {
        "id": "iF3BqddePiIi"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "We download an EMG sample file from the repository and plot the time trace. The EMG data were recorded using two differential electrodes on the biceps brachii. Please have a look at our publication of the <a href=\"https://ieeexplore.ieee.org/document/9630034\">tattoo-based EMG system</a> to get an impression of the applied test setup and electrode placement.\n"
      ],
      "metadata": {
        "id": "LAALa1ExLl9l"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "url = 'https://raw.githubusercontent.com/luuleitner/EMG_tutorial/master/data/bbcontractions.csv'\n",
        "EMG_data = pd.read_csv(url, sep=',').to_numpy()\n",
        "EMG_data = np.ravel(EMG_data[:, 0])"
      ],
      "metadata": {
        "id": "McNf4dNMqPD5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "fig = plt.figure(figsize=(6,4), dpi = 100)\n",
        "ax_1 = fig.add_subplot(111)\n",
        "ax_1.plot(EMG_data,\n",
        "          markersize = 1,\n",
        "          color='blue',\n",
        "          label = 'x')\n",
        "ax_1.set_ylabel('EMG raw [uV]')\n",
        "ax_1.grid()\n",
        "ax_1.set_xlabel('Samples')\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "GyDwRJBHRzxe"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "<h2>Assignment</h2>"
      ],
      "metadata": {
        "id": "ovNxjPiDWqjv"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<h4>1. Filtering</h4>\n",
        "\n",
        "We recorded our data with a sampling frequency of 1kHz. To exclude unwanted frequency components in our raw EMG signal we want to apply a fourth-order (n=4) bandpass filter (BPF) with cutoff frequencies at 20Hz (LOW) and 320Hz (HIGH). For the filtering, you will use the provided function `Butter()`, which requires inputs in the following order: `Butter(raw data, cutoff low, cutoff high, sampling rate, order= order number)`"
      ],
      "metadata": {
        "id": "GMiLm80SVZ7-"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Enter the necessary filter parameters."
      ],
      "metadata": {
        "id": "SjLynWVaL7S5"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "ARG_SAMPLING_RATE = ## --> enter data values <--\n",
        "\n",
        "# Bandpass Filter (BPF)\n",
        "ARG_BPF_N = ## --> enter data values <--\n",
        "ARG_BPF_LOW = ## --> enter data values <--\n",
        "ARG_BPF_HIGH = ## --> enter data values <--"
      ],
      "metadata": {
        "id": "HAUe4TyoDxCe"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Apply the provided Butterworth bandpass filter to process the EMG raw data."
      ],
      "metadata": {
        "id": "SKIoseXiMZVa"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "EMG_data_filtered = #--> enter function to filter the raw EMG signal <--"
      ],
      "metadata": {
        "id": "G6qocwmnD5mh"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Plot the raw and the filtered EMG data. Attention: the `Butter()` class returns a nested structure - use `EMG_data_filtered.data_filtered` to access the filtered signal."
      ],
      "metadata": {
        "id": "yJZuE4SlMcq3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "fig = plt.figure(figsize=(6,4), dpi = 100)\n",
        "ax_1 = fig.add_subplot(211)\n",
        "ax_1.plot(#--> enter raw data array <--,\n",
        "          markersize = 1,\n",
        "          color='blue',\n",
        "          label = 'x')\n",
        "ax_1.set_ylabel('EMG raw [uV]')\n",
        "ax_1.grid()\n",
        "\n",
        "ax_2 = fig.add_subplot(212)\n",
        "ax_2.plot(#--> enter filtered data array<--,\n",
        "          markersize = 1,\n",
        "          color='red',\n",
        "          label = 'x')\n",
        "ax_2.set_ylabel('EMG filtered [uV]')\n",
        "ax_2.set_xlabel('Samples')\n",
        "ax_2.grid()\n",
        "\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "aPHaRSJwD_z5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "<h4>2. Amplitude analyses</h4>\n",
        "\n",
        "We first rectify the filtered EMG signal for the amplitude analyses using the numpy function `np.abs()` on the filtered EMG signal. Be aware that the output of `Butter()` has a nested structure. To access the filtered EMG signal you need to enter `name_of_your_output_variable.data_filtered`.\n",
        "\n",
        "For ease of use, I have already implemented a function called `rms()` to calculate the root mean square (RMS) of the EMG signal. However, this is not magic and could easily be implemented by yourself using the formula for the RMS:\n",
        "\n",
        "\\begin{align}\n",
        "        RMS = \\sqrt (1/n \\sum_{i} x_{i}^{2})\n",
        "\\end{align}\n",
        "\n",
        "RMS...root mean square<br>\n",
        "n    ...number of samples<br>\n",
        "$x_i$...EMG value<br>"
      ],
      "metadata": {
        "id": "ObMgVHslcgz3"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Rectify the EMG data."
      ],
      "metadata": {
        "id": "FW1U41MvNEhe"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "EMG_data_filtered.data_filtered_rect = #--> enter function to rectify the filtered EMG signal <--"
      ],
      "metadata": {
        "id": "WTZG1m0uNVqh"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Plot the rectified EMG signal."
      ],
      "metadata": {
        "id": "yQVjFht_NjpK"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "fig = plt.figure(figsize=(6,4), dpi = 100)\n",
        "ax_1 = fig.add_subplot(211)\n",
        "ax_1.plot(#--> enter filtered and rectified data array <--,\n",
        "          markersize = 1,\n",
        "          color='red',\n",
        "          label = 'x')\n",
        "ax_1.set_ylabel('EMG raw [uV]')\n",
        "ax_1.set_xlabel('Samples')\n",
        "ax_1.grid()\n",
        "\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "TOZSQ8vbNgVB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Analyse the EMG signal amplitudes using RMS."
      ],
      "metadata": {
        "id": "y8wmlb6nODOO"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "RMS = #--> enter function to calculate the rms <--\n",
        "print(f'RMS = {RMS} [uV]')"
      ],
      "metadata": {
        "id": "MdDPOWZOOKAB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "<h4>3. Power spectral analyses</h4>\n",
        "\n",
        "For the spectral analyses of the EMG signal, we will use Welch's method. For more details on the computation of the power spectrum using Welch's method have a look at some of the resources here: <a href=\"https://youtu.be/YK1F0-3VvQI\">Video tutorial by Mike X. Cohan</a>, Formalism provided by <a href=\"https://ccrma.stanford.edu/~jos/sasp/Welch_s_Method.html\">Stanford Center for Computer Research in Music and Acoustics</a>.\n",
        "\n",
        "This computation is already implemented in python's `scipy` library. You can calculate the power spectrum from the filtered EMG data using the command `signal.welch()`. The function requires the following inputs: `signal.welch(filtered EMG data, sampling rate, nperseg=1024)`. The `nperseg=1024` defines the length of the used window (1024 samples). \n"
      ],
      "metadata": {
        "id": "1fu9g48kgwRi"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Evaluate the power spectral density using Welch's method."
      ],
      "metadata": {
        "id": "Qmy0b_8FUX5m"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "f, P_den = #--> enter function to calculate the power spectrum using Welch's method <--"
      ],
      "metadata": {
        "id": "tKHyIawYUV9Z"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Plot the rectified power spectrum."
      ],
      "metadata": {
        "id": "oA-reQXJlpiP"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "fig1 = plt.figure(figsize=(6,3), dpi = 300)\n",
        "ax_1 = fig1.add_subplot(111)\n",
        "ax_1.plot(f, P_den,\n",
        "          markersize = 1,\n",
        "          color='blue',\n",
        "          label = 'x')\n",
        "ax_1.set_xlabel('Frequency [Hz]')\n",
        "ax_1.set_ylabel('Power [W/Hz]')\n",
        "plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))\n",
        "ax_1.yaxis.major.formatter._useMathText = True\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "_h8vYjR7UrwP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Frequency domain features such as mean and median frequencies are frequently used to assess muscular fatigue [2]. In this next exercise we calculate these characteristics of the power spectrum and plot the results using the following formalism:\n",
        "\n",
        "\\begin{align}\n",
        "        MNF = \\sum_{i=1}^{n} f_{i}*P_{i} / \\sum_{i=1}^{n} P_{i}\n",
        "\\end{align}\n",
        "\n",
        "\\begin{align}\n",
        "      MDF = \\sum_{i=1}^{n_{MDF}}P_{i} = \\frac{1}{2} \\sum_{i=1}^{n}P_{i}\n",
        "\\end{align}\n",
        "\n",
        "MNF...mean frequency<br>\n",
        "MDF...median frequency<br>\n",
        "$P_{i}$...power spectrum at frequency of of $i^{th} sample$<br>\n",
        "$f_{i}$...frequency of $i^{th}$ sample<br>\n",
        "n...number of samples<br>\n",
        "$n_{MDF}$...sample number at median frequency"
      ],
      "metadata": {
        "id": "nisps4ful3LV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "mean_frequency = np.sum(f.T * P_den / np.sum(P_den))\n",
        "median_frequency = np.squeeze(f[np.argwhere(np.diff(np.cumsum(P_den)<=(0.5 * np.sum(P_den))))])"
      ],
      "metadata": {
        "id": "vx6Ua46iVJdZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "fig1 = plt.figure(figsize=(6,3), dpi = 300)\n",
        "ax_1 = fig1.add_subplot(111)\n",
        "ax_1.plot(f, P_den,\n",
        "          markersize = 1,\n",
        "          color='blue',\n",
        "          label = 'x')\n",
        "ax_1.set_xlabel('Frequency [Hz]')\n",
        "ax_1.set_ylabel('Power [W/Hz]')\n",
        "plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))\n",
        "ax_1.yaxis.major.formatter._useMathText = True\n",
        "\n",
        "ax_1.vlines(mean_frequency, 0, np.max(P_den), linestyles =\"dotted\", colors =\"r\")\n",
        "ax_1.annotate(f'Mean frequency = {np.round(mean_frequency,2)} [Hz]', xy=(mean_frequency, np.max(P_den)),  xycoords='data',\n",
        "            xytext=(0.8, 0.95), textcoords='axes fraction',\n",
        "            arrowprops=dict(facecolor='black', shrink=0.01, width=1, headwidth=5),\n",
        "            horizontalalignment='right', verticalalignment='top',\n",
        "            )\n",
        "\n",
        "ax_1.vlines(median_frequency, 0, np.max(P_den), linestyles =\"dashed\", colors =\"r\")\n",
        "ax_1.annotate(f'Median frequency = {np.round(median_frequency,2)} [Hz]', xy=(median_frequency, np.max(P_den)/2),  xycoords='data',\n",
        "            xytext=(0.8, 0.5), textcoords='axes fraction',\n",
        "            arrowprops=dict(facecolor='black', shrink=0.01, width=1, headwidth=5),\n",
        "            horizontalalignment='right', verticalalignment='top',\n",
        "            )\n",
        "\n",
        "\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "5FLcQDy1amK9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Great job! \n",
        "## You have completed the EMG tutorial\n",
        "\n",
        "For more open-source tools visit my github account at: <a href=\"https://github.com/luuleitner\">github.com/luuleitner</a><br><br>\n"
      ],
      "metadata": {
        "id": "E-gAXBu3mW2Y"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<h2>References</h2>\n",
        "\n",
        "\n",
        "[1] Leitner, C. et al. (2021). UStEMG: an Ultrasound Transparent Tattoo-based sEMG System for Unobtrusive Parallel Acquisitions of Muscle Electro-mechanics. in Proc. 43nd Conferences of the IEEE Engineering in Medicine and Biology Society, Guadalajara, Mexico. DOI: <a href=\"https://ieeexplore.ieee.org/document/9630034\">10.1109/EMBC46164.2021.9630034</a>\n",
        "\n",
        "[2] Phinyomark, A. et al. (2012). The Usefulness of Mean and Median Frequencies in Electromyography Analysis. In Computational Intelligence in Electromyography Analysis – A Perspective on Current Applications and Future Challenges, edited by Ganesh R. Naik. InTech. DOI: <a href=\"https://doi.org/10.5772/50639\">10.5772/50639</a>"
      ],
      "metadata": {
        "id": "PtVdCjBhhf3M"
      }
    }
  ]
}