notebooks/imu-stepcount.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysing steps in a short walk using acceleration and rotation\n",
    "\n",
    "In this noteboook we examine a recording from the Arduino IMU of a short walk with the aim to extract the step count, cadence and timings including ground time and step duration.\n",
    "\n",
    "We implement three techniques outlined in the paper `A comprehensive comparison of simple step counting techniques using wrist- and ankle-mounted accelerometer and gyroscope signals` by Matthew Rudy and Joseph Mahoney - [https://www.researchgate.net/publication/325451208_A_comprehensive_comparison_of_simple_step_counting_techniques_using_wrist-_and_ankle-mounted_accelerometer_and_gyroscope_signals](https://www.researchgate.net/publication/325451208_A_comprehensive_comparison_of_simple_step_counting_techniques_using_wrist-_and_ankle-mounted_accelerometer_and_gyroscope_signals).\n",
    "\n",
    "* Peak-finding\n",
    "* Fast Fourier Transform (FFT)\n",
    "* Autocorrelation\n",
    "\n",
    "Each of these methods allows us to count steps. The peak-finding method also identifies where the steps occur in the timeseries, so this in turn allows us to isolate steps and calculate such things as ground time and step duration.\n",
    "\n",
    "## The IMU\n",
    "\n",
    "The Arduino Nano inertial measurement unit gives us acceleration, measured in `g`s (`1g = 9.8m/s/s`), and rotation (angular velocity). In this notebook we demonstrate how to extract steps from either acceleration or rotation.\n",
    "\n",
    "## The data\n",
    "\n",
    "We expect a CSV file with columns for time, 3 axes of acceleration, and 3 axes of gyroscopic rotation.\n",
    "\n",
    "Example:\n",
    "```\n",
    "time,aX,aY,aZ,gX,gY,gZ\n",
    "363328,1.0157471,-0.0056152344,0.032714844,1.0986328,-4.8828125,-3.540039\n",
    "363338,1.0213623,0.009643555,0.017700195,1.8310547,-3.540039,-4.0893555\n",
    "363348,1.015625,0.01940918,0.020874023,1.8920898,-2.1362305,-4.3945312\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from scipy.integrate import cumtrapz\n",
    "from scipy.signal import butter, filtfilt, periodogram, spectrogram, find_peaks\n",
    "from sklearn.decomposition import PCA\n",
    "import matplotlib.pyplot as plt\n",
    "sns.set(rc={'figure.figsize':(11, 4)})\n",
    "import numpy as np\n",
    "import gpxpy\n",
    "from xml.etree import ElementTree as ET\n",
    "from datetime import timedelta\n",
    "from tqdm.notebook import tqdm\n",
    "tqdm()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_file_path = \"../data/20200914 mellor sports field/mellor sports field circles (wmf).csv\"\n",
    "output_file_path = \"../data/20200914 mellor sports field/mellor sports field circles (wmf).gpx\"\n",
    "df = pd.read_csv(data_file_path)\n",
    "df.columns = [\"time\", \"aX\", \"aY\", \"aZ\", \"gX\", \"gY\", \"gZ\"]\n",
    "df.time = df.time - df.time.min()\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Key parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`fs` is a sampling frequency of an IMU used (i.e. how many times per second accelerationa and angular velocity are measured). The default delay between readings (defined in the Arduino code) is 10 ms, which corresponds to the sampling frequency of 100 Hz. The value of sampling frequency needs to be changed accordingly if different rate was used during the recording (either with a different device, or with modified Arduino code)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 1000 / 10 # Hz, sampling frequency\n",
    "total_time = df.time.max() - df.time.min()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_time_period(a, b):\n",
    "    \"\"\"\n",
    "    a, b -- time in seconds\n",
    "    \"\"\"\n",
    "    return df.loc[(df.time >= a * 1000) & (df.time < b * 1000)]\n",
    "\n",
    "def filter_series(data, fc_low=None, fc_high=None, fs=fs):\n",
    "    if fc_low is not None and fc_high is not None:\n",
    "        return band_pass(data, fc_low, fc_high, fs)\n",
    "    elif fc_low is not None:\n",
    "        return high_pass(data, fc_low, fs)\n",
    "    elif fc_high is not None:\n",
    "        return low_pass(data, fc_high, fs)\n",
    "    else:\n",
    "        return data\n",
    "        \n",
    "def band_pass(data, fc_low, fc_high, fs):\n",
    "    w_low = fc_low / (fs / 2) # Normalize the frequency\n",
    "    w_high = fc_high / (fs / 2) # Normalize the frequency\n",
    "    b, a = butter(5, [w_low, w_high], 'bandpass')\n",
    "    return filtfilt(b, a, data)\n",
    "\n",
    "def high_pass(data, fc_low, fs):\n",
    "    w_low = fc_low / (fs / 2) # Normalize the frequency\n",
    "    b, a = butter(5, w_low, 'highpass')\n",
    "    return filtfilt(b, a, data)\n",
    "\n",
    "def low_pass(data, fc_high, fs):\n",
    "    w_high = fc_high / (fs / 2) # Normalize the frequency\n",
    "    b, a = butter(5, w_high, 'lowpass')\n",
    "    return filtfilt(b, a, data)\n",
    "\n",
    "def get_magnitude(data, fc_low=None, fc_high=None, fs=None):\n",
    "    magnitude = np.sqrt((data**2).sum(axis=1))\n",
    "    if fc_low is not None and fc_high is not None:\n",
    "        return band_pass(magnitude, fc_low, fc_high, fs)\n",
    "    elif fc_low is not None:\n",
    "        return high_pass(magnitude, fc_low, fs)\n",
    "    elif fc_high is not None:\n",
    "        return low_pass(magnitude, fc_high, fs)\n",
    "    else:\n",
    "        return magnitude\n",
    "    \n",
    "def peak_detection_steps(data, pos_kwargs=None, neg_kwargs=None, plot=False):\n",
    "    peaks, _ = find_peaks(data, **pos_kwargs)\n",
    "    neg_peaks, _ = find_peaks(-data, **neg_kwargs)\n",
    "    if plot:\n",
    "        sns.lineplot(x=range(len(data)), y=data)\n",
    "        sns.scatterplot(x=peaks, y=data[peaks])\n",
    "        sns.scatterplot(x=neg_peaks, y=data[neg_peaks])\n",
    "        plt.show()\n",
    "#     return len(peaks)\n",
    "    return max(len(peaks), len(neg_peaks))\n",
    "\n",
    "def fft_dominant_freq(data, fs, plot=False):\n",
    "    f, Pxx = periodogram(data, fs=fs)\n",
    "    if plot:\n",
    "        sns.lineplot(f, Pxx)\n",
    "        plt.xlim([0.0,10.0])\n",
    "        plt.show()\n",
    "    return f[np.argmax(Pxx)]\n",
    "\n",
    "def fft_steps(data, dt, fs, plot=False):\n",
    "    return dt * fft_dominant_freq(data, fs, plot) / 1000\n",
    "\n",
    "def autocorr(x):\n",
    "    result = np.correlate(x, x, mode='full')\n",
    "    return result[:int(len(result)/2)]\n",
    "\n",
    "def autocorr_steps(data, plot=False):\n",
    "    corr = autocorr(data)\n",
    "    peaks, _ = find_peaks(corr)\n",
    "    if plot:\n",
    "        sns.lineplot(x=range(len(corr)), y=corr)\n",
    "        sns.scatterplot(x=peaks, y=corr[peaks])\n",
    "        plt.show()\n",
    "    return len(peaks)\n",
    "\n",
    "def to_steps_per_minute(step_count, dt):\n",
    "    \"\"\"\n",
    "    dt -- time in seconds\n",
    "    \"\"\"\n",
    "    return step_count / dt * 60\n",
    "\n",
    "def get_spm_for_period(a, b, columns, fc_low=None, fc_high=None, peak_detection_kwargs={}):\n",
    "    frame = get_time_period(a, b)\n",
    "    dt = (frame.time.max() - frame.time.min()) / 1000.0\n",
    "    aMagnitude = get_magnitude(frame.loc[:,columns], fc_low=fc_low, fc_high=fc_high, fs=fs)\n",
    "    n_peak_steps = peak_detection_steps(\n",
    "        aMagnitude, \n",
    "        **peak_detection_kwargs\n",
    "    )\n",
    "    dominant_freq = fft_dominant_freq(aMagnitude, fs=fs)\n",
    "    n_autocorr_steps = autocorr_steps(aMagnitude)\n",
    "    return pd.Series({\n",
    "        \"peak_detection_spm\": to_steps_per_minute(n_peak_steps, dt), \n",
    "        \"fft_spm\": dominant_freq * 60, \n",
    "        \"autocorrelation_spm\": to_steps_per_minute(n_autocorr_steps, dt), \n",
    "    })\n",
    "\n",
    "def get_spm(data, a, b, peak_detection_fc=(None, None), fft_fc=(None, None), autocorrelation_fc=(None, None), peak_detection_kwargs={}):\n",
    "    frame = get_time_period(a, b)\n",
    "    dt = (frame.time.max() - frame.time.min()) / 1000.0\n",
    "    \n",
    "    _data = filter_series(data.loc[frame.index].values, *peak_detection_fc)\n",
    "        \n",
    "    n_peak_steps = peak_detection_steps(\n",
    "        _data, \n",
    "        **peak_detection_kwargs\n",
    "    )\n",
    "    \n",
    "    _data = filter_series(data.loc[frame.index].values, *fft_fc)\n",
    "    dominant_freq = fft_dominant_freq(_data, fs=fs)\n",
    "    \n",
    "    _data = filter_series(data.loc[frame.index].values, *autocorrelation_fc)\n",
    "    n_autocorr_steps = autocorr_steps(_data)\n",
    "    return pd.Series({\n",
    "        \"peak_detection_spm\": to_steps_per_minute(n_peak_steps, dt), \n",
    "        \"fft_spm\": dominant_freq * 60, \n",
    "        \"autocorrelation_spm\": to_steps_per_minute(n_autocorr_steps, dt), \n",
    "    })\n",
    "\n",
    "def scalar_projection(u, v):\n",
    "    \"\"\"\n",
    "    project u on v\n",
    "    \"\"\"\n",
    "    v_norm = np.sqrt(sum(v**2))\n",
    "    return (np.dot(u, v)/v_norm)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Measuring steps from acceleration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore\n",
    "We first look at what data we have. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.lineplot(x=df.time, y=df.aX, label=\"aX\")\n",
    "sns.lineplot(x=df.time, y=df.aY, label=\"aY\")\n",
    "sns.lineplot(x=df.time, y=df.aZ, label=\"aZ\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zooming in, we can clearly see how different components of acceleration contribute to a step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_df = df.iloc[4050:4090].loc[:, [\"aX\", \"aY\", \"aZ\"]]\n",
    "_df.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "For further calculations, we need to find and remove the gravity vector. We can do this by taking a mean average of the first readings (assuming the user was initially standing still)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_df = df.iloc[:9]\n",
    "sns.lineplot(x=_df.time, y=_df.aX, label=\"aX\")\n",
    "sns.lineplot(x=_df.time, y=_df.aY, label=\"aY\")\n",
    "sns.lineplot(x=_df.time, y=_df.aZ, label=\"aZ\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gravity_vector = _df.loc[:,[\"aX\", \"aY\", \"aZ\"]].mean()\n",
    "(gravity_vector, np.sqrt((gravity_vector**2).sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting vector is of the expected direction (primarily in Z+ direction) and the expected magnitude (roughly 1 `g`). Subtracting it from all the acceleration readings leaves only the unbalanced acceleration/force (i.e. the acceleration that contributes to motion). We also multiply acceleration values by 9.8 to convert it to m/s^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[:,[\"aX\", \"aY\", \"aZ\"]] = (df.loc[:,[\"aX\", \"aY\", \"aZ\"]] - gravity_vector) * 9.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we do not know the direction of movement (the device can be placed in different orientation), we perform Principal component analysis (PCA). The 0th component should ideally lie in the direction of movement, assuming there's a greater variation in acceleration in this direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[:,[\"pca0\", \"pca1\", \"pca2\"]] = PCA().fit_transform(df.loc[:,[\"aX\", \"aY\", \"aZ\"]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we expected, the 0th component shows two clearly defined peaks that correspond to a step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_df = df.iloc[4050:4090].loc[:, [\"pca0\", \"pca1\", \"pca2\"]]\n",
    "_df.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the spectrogram of the 0th component, which shows a signal at ~1.3 Hz (roughly 78 steps per minute). However the second harmonic of this signal (at ~2.6 Hz or 156 steps per minute) is more prominent, which is the expected value for cadence measured from both feet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, t, Sxx = spectrogram(df.pca0, fs)\n",
    "plt.pcolormesh(t, f, Sxx, shading='gouraud')\n",
    "plt.ylabel('Frequency [Hz]')\n",
    "plt.xlabel('Time [sec]')\n",
    "plt.ylim([0,10])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cadence calculations\n",
    "We use three methods to estimate cadence from 0th component:\n",
    "\n",
    "* Peak detection -- count all peaks (positive and negative) with minimum height of 35 m/s^2 and distance of 20 readings between (~200 ms)\n",
    "* FFT (Fast Fourier Transform) -- perform FFT (i.e. transform signal from time domain to a frequency domain) and take a frequency with the highest power\n",
    "* Autocorrelation -- perform a correlation of a signal with a delayed copy of itself (this reveals periodicity of a signal) and count the peaks.\n",
    "\n",
    "Additionally, we filter the signal differently for different methods. We apply bandpass filter with frequencies 0.5-4.0 Hz for FFT, and 0.3-4.0 Hz for Autocorrelation. No filtering is done before Peak detection, as it makes the peaks that we've seen before less prominent and makes it harder to detect them.\n",
    "\n",
    "Applying the methods to a short segment of 5 seconds gives 85.7, 192.0 and 160.7 steps per minute respectively. Note: since the second harmonic is more prominent, FFT and autocorrelation give cadence for both feet (i.e. double of peak detection method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 4500\n",
    "time = df.time[a:a+500]\n",
    "total_time = time.max() - time.min()\n",
    "display(to_steps_per_minute(peak_detection_steps(\n",
    "    df.pca0[a:a+500].values, \n",
    "    pos_kwargs={\n",
    "        \"distance\": 20,\n",
    "        \"height\": (35, None)\n",
    "    }, \n",
    "    neg_kwargs={\n",
    "        \"distance\": 20,\n",
    "        \"height\": (35, None)\n",
    "    },\n",
    "    plot=True\n",
    "), total_time / 1000))\n",
    "display(to_steps_per_minute(fft_steps(band_pass(df.pca0[a:a+500], fc_low=0.5, fc_high=4.0, fs=fs), dt=total_time, fs=fs, plot=True), total_time / 1000))\n",
    "display(to_steps_per_minute(autocorr_steps(band_pass(df.pca0[a:a+500], fc_low=0.3, fc_high=4.0, fs=fs), plot=True), total_time / 1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying these algorithms to windows of 10 seconds, we get a series of cadence measurements over time. The expected mean cadence (measured with Stryd) is ~82 spm. The results of our calculations come close to that with 71.5, 170.2 and 154.5 spm for Peak detection, FFT, and Autocorrelation respectively (as before FFT and Autocorrelation gives cadence for both feet, hence double).\n",
    "\n",
    "The cadence line plot looks choppy, since the number of steps measured is a discrete. However it can be smoothen using Exponential Moving Average (see lines labeled with EMA suffix)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "peak_detection_kwargs = {\n",
    "    \"pos_kwargs\": {\n",
    "        \"prominence\": 20,\n",
    "        \"distance\": 20,\n",
    "        \"height\": (35, None)\n",
    "    }, \n",
    "    \"neg_kwargs\": {\n",
    "        \"prominence\": 20,\n",
    "        \"distance\": 20,\n",
    "        \"height\": (35, None)\n",
    "    },\n",
    "}\n",
    "\n",
    "spms = []\n",
    "dt = 10\n",
    "step = 1.0\n",
    "for a in tqdm(np.arange(0, df.time.max() / 1000, step)):\n",
    "    b = a + dt\n",
    "    spm = get_spm(df.pca0, a, b, fft_fc=(0.5, 4.0), autocorrelation_fc=(0.3, 4.0), peak_detection_kwargs=peak_detection_kwargs)\n",
    "    spm.name = a\n",
    "    spms += [spm]\n",
    "#     break\n",
    "acceleration_spm_df = pd.DataFrame(spms)\n",
    "acceleration_spm_df[\"peak_detection_spm_EMA\"] = acceleration_spm_df.peak_detection_spm.ewm(alpha=0.3, adjust=False).mean()\n",
    "acceleration_spm_df[\"autocorrelation_spm_EMA\"] = acceleration_spm_df.autocorrelation_spm.ewm(alpha=0.3, adjust=False).mean()\n",
    "acceleration_spm_df.plot()\n",
    "display(acceleration_spm_df.describe())\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Measuring steps from rotation\n",
    "The same approach can be taken to calculate cadence from gyroscope data (angular velocity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.lineplot(x=df.time, y=df.gX, label=\"gX\")\n",
    "sns.lineplot(x=df.time, y=df.gY, label=\"gY\")\n",
    "sns.lineplot(x=df.time, y=df.gZ, label=\"gZ\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_df = df.iloc[4050:4090].loc[:, [\"gX\", \"gY\", \"gZ\"]]\n",
    "_df.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of principal components, we are using the total magnitude of angular velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_magnitude = get_magnitude(df.loc[:,[\"gX\", \"gY\", \"gZ\"]].iloc[4050:4090], fs=fs)\n",
    "sns.lineplot(x=df.time.iloc[4050:4090], y=_magnitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though the signal is weaker for angular velocity magnitude, the peaks at 1.3 and 2.6 Hz are still detectable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, t, Sxx = spectrogram(get_magnitude(df.loc[:,[\"gX\", \"gY\", \"gZ\"]], fs=fs), fs)\n",
    "plt.pcolormesh(t, f, Sxx, shading='gouraud')\n",
    "plt.ylabel('Frequency [Hz]')\n",
    "plt.xlabel('Time [sec]')\n",
    "plt.ylim([0,10])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And hence, applying the cadence measurement algorithms as before, we get the mean cadence of 65.9, 89.3, and 75.8 spm for Peak detection, FFT and Autocorrelation respectively. Which are close to the expected ~82 spm measured with Stryd."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "peak_detection_kwargs = {\n",
    "    \"pos_kwargs\": {\n",
    "        \"height\": (50, None)\n",
    "    }, \n",
    "    \"neg_kwargs\": {\n",
    "        \"height\": (50, None)\n",
    "    },\n",
    "}\n",
    "\n",
    "spms = []\n",
    "dt = 10\n",
    "step = 1\n",
    "for a in np.arange(0, df.time.max() / 1000, step):\n",
    "    b = a + dt\n",
    "    spm = get_spm_for_period(a, b, [\"gX\", \"gY\", \"gZ\"], fc_low=0.5, fc_high=2.0, peak_detection_kwargs=peak_detection_kwargs)\n",
    "    spm.name = a\n",
    "    spms += [spm]\n",
    "spm_df = pd.DataFrame(spms)\n",
    "spm_df.plot()\n",
    "spm_df.describe()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Export GPX\n",
    "The resulting can be exported as a GPX track, to use it other tools that support GPX (e.g. DC Rainmaker Analyzer). By default Peak detection cadence calculated from acceleration and smoothen with EMA is exported. \n",
    "\n",
    "Note: the timestamp of the first reading needs to be provided, since GPX assumes that all timestamps are absolute, but the Wearable My Foot device only gives relative timestamps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_cadence_extension(cadence):\n",
    "    prefix = \"gpxtrx:\"\n",
    "    element = ET.Element(f\"{prefix}TrackPointExtension\")\n",
    "    cadence_element = ET.SubElement(element, f\"{prefix}cad\")\n",
    "    # Schema only permits integers up to 254\n",
    "    cadence_element.text = str(int(cadence if cadence <= 254 else 254))\n",
    "    return element\n",
    "\n",
    "def get_point(time, cadence):\n",
    "    extensions = [get_cadence_extension(cadence)]\n",
    "    point = gpxpy.gpx.GPXTrackPoint()\n",
    "    point.extensions = extensions\n",
    "    point.time = time\n",
    "    return point\n",
    "\n",
    "def get_gpx(timestamps, cadence):\n",
    "    \"\"\"\n",
    "    data -- pandas DataFrame with time and cadence fields\n",
    "    \"\"\"\n",
    "    gpx = gpxpy.gpx.GPX()\n",
    "    gpx.nsmap[\"gpxtrx\"] = 'http://www.garmin.com/xmlschemas/GpxExtensions/v3'\n",
    "    track = gpxpy.gpx.GPXTrack()\n",
    "    gpx.tracks.append(track)\n",
    "    segment = gpxpy.gpx.GPXTrackSegment()\n",
    "    track.segments.append(segment)\n",
    "    segment.points = [get_point(t, c) for t, c in zip(timestamps, cadence)]\n",
    "    return gpx\n",
    "    \n",
    "    \n",
    "timestamp = pd.to_datetime(\"2020-09-14T17:07:14\")\n",
    "timestamps = [timestamp + timedelta(seconds=x) for x in acceleration_spm_df.index]\n",
    "cadence = acceleration_spm_df.peak_detection_spm_EMA\n",
    "\n",
    "with open(output_file_path, 'w+') as f:\n",
    "    f.write(get_gpx(timestamps, cadence).to_xml())"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}