{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Heartbeat Evoked Arpeggios (cardiac-contingent stimuli)\n\nThis tutorial illustrates how to use the ``Oximeter`` class to trigger stimuli\nat different phases of the cardiac cycle. The PPG signal is recorded for 30 seconds\nand peaks are detected online. Four notes ('C', 'E', 'G', 'Bfl') are played in synch\nwith peak detection with various delays: no delay,  1/4, 2/4 or 3/4 of the previous\ncardiac cycle length. While R-R intervals are prone to large changes over longer\ntimescales, such changes are physiologically limited from one heartbeat to the next,\nlimiting variance in the onset synchrony between the tones and the cardiac cycle.\nOn this basis, each presentation time is calibrated based on the previous RR-interval.\nThis procedure can easily be adapted to create a standard interoception task, e.g. by\neither presenting tones at no delay (systole, s+) or a fixed offset (diastole, s-).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Author: Nicolas Legrand <nicolas.legrand@cfin.au.dk>\n# Licence: Apache-2.0\n\nimport itertools\nimport time\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nfrom systole import serialSim\nfrom systole.plots import plot_circular, plot_events\nfrom systole.recording import Oximeter\nfrom systole.utils import norm_triggers, to_angles"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Recording\nFor the purpose of demonstration, here we simulate data acquisition through\nthe pulse oximeter using pre-recorded signal.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ser = serialSim()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If you want to allow online data acquisition, you should uncomment the\nfollowing lines and provide the reference of the COM port where the pulse\noximeter is plugged in.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "```python\nimport serial\nser = serial.Serial('COM4')  # Change this value according to your setup\n```\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create an Oximeter instance, initialize recording and record for 10 seconds\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "oxi = Oximeter(serial=ser, sfreq=75, add_channels=4).setup()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Play the different tones at different intervals following the detection of the\nsystolic peaks, using the length of the previous RR interval as reference. To kee the\nexecution of the example simple we do no pla atual sounds (e.g. using Psychopy) but\nthis step can easily be added by adding new lines for each time window.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "systoleTime1, systoleTime2, systoleTime3 = None, None, None\ntstart = time.time()\nwhile time.time() - tstart < 20:\n\n    # Check if there are new data to read\n    if oxi.serial.inWaiting() >= 5:\n\n        # Read one incoming data point - Convert bytes into list of int\n        paquet = list(oxi.serial.read(5))\n\n        if oxi.check(paquet):  # Data consistency\n            oxi.add_paquet(paquet[2])  # Add new data point\n\n        # T + 0\n        if oxi.peaks[-1] == 1:\n            oxi.channels[\"Channel_0\"][-1] = 1\n            systoleTime1 = time.time()\n            systoleTime2 = time.time()\n            systoleTime3 = time.time()\n\n        # T + 1/4\n        if systoleTime1 is not None:\n            if time.time() - systoleTime1 >= ((oxi.instant_rr[-1] / 4) / 1000):\n                oxi.channels[\"Channel_1\"][-1] = 1\n                systoleTime1 = None\n\n        # T + 2/4\n        if systoleTime2 is not None:\n            if time.time() - systoleTime2 >= (((oxi.instant_rr[-1] / 4) * 2) / 1000):\n                oxi.channels[\"Channel_2\"][-1] = 1\n                systoleTime2 = None\n\n        # T + 3/4\n        if systoleTime3 is not None:\n            if time.time() - systoleTime3 >= (((oxi.instant_rr[-1] / 4) * 3) / 1000):\n                oxi.channels[\"Channel_3\"][-1] = 1\n                systoleTime3 = None"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Events\nWe now have triggers for the presentation of different notes relative to the systolic \npeaks detected in the PPG signal in real-time. We can visualize the distribution\nof these events using the `plot_events` function that will take the cleaned triggers\nchannels as input. We can see in the plot below that the events are nicely aligned\nwith the PPG signal.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "_, axs = plt.subplots(2, 1, figsize=(18, 6), sharex=True)\n\n# Loop across channels\nevents_labels = [\"C\", \"E\", \"G\", \"Bfl\"]\npalette = itertools.cycle(sns.color_palette(\"deep\"))\nfor i in np.arange(3, -1, -1):\n    \n    # Clan the triggers vector so we just have one trigger / event\n    triggers = norm_triggers(\n        triggers=oxi.channels[f\"Channel_{i}\"], n=10, direction=\"higher\"\n        )\n\n    # Plot the event according to its duration (0.1 second)\n    plot_events(\n        triggers=triggers, tmin=0.0, tmax=0.1, sfreq=75, \n        labels=events_labels[i], ax=axs[0], palette=[next(palette)])\n    axs[0].set_xlabel(\"\")\n\noxi.plot_raw(ax=axs[1])\n\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Cardiac cycle\nAnother way to have a look at this event distribution is to look at when they occur\nrelative to the cardiac cycle. Considering that one full cycle (beat to beat) can be\nrepresented as the full revolution around a trigonometric circle, it is possible to\nshow the occurrence of an event relative to this cardiac phase.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Creat a list of angular values for each conditions\nangles = []\nx = np.asarray(oxi.peaks)\nfor ev in oxi.channels:\n    events = norm_triggers(\n        np.asarray(oxi.channels[ev]), threshold=1, n=10, direction=\"higher\"\n    )\n    angles.append(to_angles(np.where(x)[0], np.where(events)[0]))\n\n# Plot the cardiac phase for each condition using different density representation\n_, axs = plt.subplots(1, 3, figsize=(13, 5), subplot_kw=dict(projection=\"polar\"))\nfor i, density in enumerate([\"height\", \"area\", \"alpha\"]):\n    plot_circular(angles, density=density, ax=axs[i])\n    axs[i].set_title(f\"Density method = {density}\")"
      ]
    }
  ],
  "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.11.16"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}