{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# AMPD - Automatic Multiscale Peak Detection\n", "\n", "This notebook provides an end-to-end pipeline for processing and analyzing fNIRS data collected during a finger-tapping task. The primary goal is to identify peaks in the time series data using an **Optimized AMPD** algorithm.\n", "\n", "The **AMPD** algorithm is a multiscale peak detection technique that is especially effective for periodic and quasi-periodic signals, such as heart beats, even in the presence of noise. By analyzing the signal at multiple scales, the algorithm can reliably detect local maxima while minimizing false positives. This method is based on the work by **[Scholkmann et al. 2012](https://doi.org/10.3390/a5040588)**\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2025-06-06T16:13:01.264907Z", "iopub.status.busy": "2025-06-06T16:13:01.264733Z", "iopub.status.idle": "2025-06-06T16:13:02.909160Z", "shell.execute_reply": "2025-06-06T16:13:02.908696Z" } }, "outputs": [], "source": [ "import cedalion.nirs\n", "from cedalion import units\n", "from cedalion.sigproc import quality\n", "from cedalion.sigproc.frequency import freq_filter\n", "import cedalion.xrutils as xrutils\n", "from cedalion.datasets import get_fingertapping_snirf_path\n", "import time\n", "import numpy as np\n", "import xarray as xr\n", "from cedalion.sigproc.physio import ampd\n", "import matplotlib.pyplot as plt\n", "\n", "xr.set_options(display_max_rows=3, display_values_threshold=50)\n", "np.set_printoptions(precision=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Loading raw CW-NIRS data from a SNIRF file\n", "This notebook uses a finger-tapping dataset in BIDS layout provided by [Rob Luke](https://github.com/rob-luke/BIDS-NIRS-Tapping). It can can be downloaded via `cedalion.datasets`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load amplitude data from the snirf file and extract the first 60 seconds for further processing" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2025-06-06T16:13:02.911929Z", "iopub.status.busy": "2025-06-06T16:13:02.911365Z", "iopub.status.idle": "2025-06-06T16:13:03.020603Z", "shell.execute_reply": "2025-06-06T16:13:03.020176Z" } }, "outputs": [], "source": [ "path_to_snirf_file = get_fingertapping_snirf_path()\n", "\n", "recordings = cedalion.io.read_snirf(path_to_snirf_file)\n", "rec = recordings[0] # there is only one NirsElement in this snirf file...\n", "amp = rec[\"amp\"] # ... which holds amplitude data\n", "\n", "# restrict to first 60 seconds and fill in missing units\n", "amp = amp.sel(time=amp.time < 60)\n", "times = amp.time.values * 1000\n", "# print(amp.time.values[-1] / 60, len(times))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Following are utility methods for normalizing, filtering and plotting the signal" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2025-06-06T16:13:03.022937Z", "iopub.status.busy": "2025-06-06T16:13:03.022746Z", "iopub.status.idle": "2025-06-06T16:13:03.026654Z", "shell.execute_reply": "2025-06-06T16:13:03.026281Z" } }, "outputs": [], "source": [ "\n", "# collection of utility functions\n", "\n", "def normalize(sig):\n", " min_val = np.min(sig)\n", " max_val = np.max(sig)\n", " return (sig - min_val) / (max_val - min_val)\n", "\n", "def filter_signal(amplitudes):\n", " return freq_filter(amplitudes, 0.5 * units.Hz, 3 * units.Hz, 2)\n", "\n", "def plot_peaks(signal, s_times, s_peaks, label, title='peaks'):\n", " fig, ax = plt.subplots(1, 1, figsize=(24, 8))\n", " ax.plot(s_times, signal, label=label)\n", "\n", " for ind, peak in enumerate(s_peaks):\n", " if peak > 0:\n", " ax.axvline(x=peak, color='black', linestyle='--', linewidth=1)\n", "\n", " plt.title(title)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### This is the amplitude data structure" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2025-06-06T16:13:03.028592Z", "iopub.status.busy": "2025-06-06T16:13:03.028278Z", "iopub.status.idle": "2025-06-06T16:13:03.041646Z", "shell.execute_reply": "2025-06-06T16:13:03.041160Z" } }, "outputs": [ { "data": { "text/html": [ "
<xarray.DataArray (channel: 28, wavelength: 2, time: 469)> Size: 210kB\n", "<Quantity([[[0.0914 0.091 0.091 ... 0.0903 0.0902 0.0899]\n", " [0.1857 0.1864 0.1837 ... 0.1849 0.185 0.1847]]\n", "\n", " [[0.2275 0.2297 0.2261 ... 0.2241 0.2243 0.2257]\n", " [0.6355 0.6377 0.6298 ... 0.6223 0.6237 0.6272]]\n", "\n", " [[0.1065 0.1066 0.1053 ... 0.1065 0.1062 0.1056]\n", " [0.2755 0.2762 0.2727 ... 0.2737 0.2742 0.276 ]]\n", "\n", " ...\n", "\n", " [[0.2028 0.1997 0.2005 ... 0.1998 0.2007 0.2026]\n", " [0.4666 0.4554 0.4562 ... 0.4482 0.4511 0.4541]]\n", "\n", " [[0.4885 0.4802 0.4818 ... 0.5005 0.5036 0.5045]\n", " [0.8458 0.826 0.826 ... 0.8386 0.8441 0.8475]]\n", "\n", " [[0.6305 0.6284 0.6287 ... 0.6373 0.638 0.6392]\n", " [1.2286 1.2206 1.219 ... 1.2232 1.2259 1.2278]]], 'dimensionless')>\n", "Coordinates: (3/6)\n", " * time (time) float64 4kB 0.0 0.128 0.256 0.384 ... 59.65 59.78 59.9\n", " samples (time) int64 4kB 0 1 2 3 4 5 6 7 ... 462 463 464 465 466 467 468\n", " ... ...\n", " * wavelength (wavelength) float64 16B 760.0 850.0\n", "Attributes:\n", " data_type_group: unprocessed raw
<xarray.DataArray (channel: 28, wavelength: 2, time: 469)> Size: 210kB\n", "array([[[0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0]],\n", "\n", " [[0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0]],\n", "\n", " [[0, 0, 0, ..., 1, 1, 1],\n", " [0, 0, 0, ..., 0, 0, 0]],\n", "\n", " ...,\n", "\n", " [[0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0]],\n", "\n", " [[0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0]],\n", "\n", " [[0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0]]])\n", "Coordinates: (3/6)\n", " * time (time) float64 4kB 0.0 0.128 0.256 0.384 ... 59.65 59.78 59.9\n", " samples (time) int64 4kB 0 1 2 3 4 5 6 7 ... 462 463 464 465 466 467 468\n", " ... ...\n", " * wavelength (wavelength) float64 16B 760.0 850.0\n", "Attributes:\n", " data_type_group: unprocessed raw