tutorial_audioproblems_interpeak_detection_estrategies.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpeak detection estrategies\n",
    "this notebook demostrates some of the intersample detection technuques: **throw resampling** and throw **parabolic interpolation**.  The accuracy of both methods can be tested on real time by shifting a sinc function from the sampling point and evaluating the error introduced by both systems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipywidgets as wg\n",
    "import essentia.standard as es\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "matplotlib.use('TkAgg')\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "\n",
    "duration = 10 # s\n",
    "fs = 1 # hz \n",
    "k = 1. # amplitude\n",
    "oversamplingFactor = 4 # factor of oversampling for the real signal\n",
    "nSamples = fs * duration\n",
    "\n",
    "time = np.arange(-nSamples/2, nSamples/2,\n",
    "                 2 ** -oversamplingFactor, dtype='float')\n",
    "samplingPoints = time[::2 ** oversamplingFactor]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shifted_sinc(x, k, offset):\n",
    "    xShifted = x - offset\n",
    "    y = np.zeros(len(xShifted))\n",
    "    for idx, i in enumerate(xShifted):\n",
    "        if not i: \n",
    "            y[idx] = k\n",
    "        else:\n",
    "            y[idx] = (k * np.sin(np.pi * i) / (np.pi * i))\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def resampleStrategy(y, fs, quality=0, oversampling=4):\n",
    "    yResample = es.Resample(inputSampleRate=fs,\n",
    "                            outputSampleRate=fs*oversampling, \n",
    "                            quality=quality)(y.astype(np.float32))\n",
    "    \n",
    "    tResample = np.arange(np.min(samplingPoints), np.max(samplingPoints) \n",
    "                          + 1, 1. / (fs * oversampling))\n",
    "    tResample = tResample[:len(yResample)]        \n",
    "    \n",
    "    # getting the stimated peaks\n",
    "    yResMax = np.max(yResample)\n",
    "    tResMax = tResample[np.argmax(yResample)]\n",
    "    \n",
    "    return yResample, tResample, yResMax, tResMax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parabolicInterpolation(y, threshold=.6):\n",
    "    # todo plot the parabol maybe\n",
    "    positions, amplitudes = es.PeakDetection(threshold=threshold)\\\n",
    "                                        (y.astype(np.float32))\n",
    "       \n",
    "    pos = int(positions[0] * (len(y-1)))\n",
    "    a = y[pos - 1]\n",
    "    b = y[pos]\n",
    "    c = y[pos + 1]\n",
    "\n",
    "    tIntMax = samplingPoints[pos] + (a - c) / (2 * (a - 2 * b + c))\n",
    "    yIntMax = b - ((a - b) ** 2) / (8 * (a - 2 * b + c))\n",
    "    return tIntMax, yIntMax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process():\n",
    "    \n",
    "    ## Processing\n",
    "    \n",
    "    # \"real\" sinc\n",
    "    yReal = shifted_sinc(time, k, offset.value)\n",
    "    \n",
    "    # sampled sinc\n",
    "    y = shifted_sinc(samplingPoints, k, offset.value)\n",
    "    \n",
    "    \n",
    "    # Resample strategy\n",
    "    yResample, tResample, yResMax, tResMax = \\\n",
    "        resampleStrategy(y, fs, quality=0, oversampling=4)\n",
    "    \n",
    "    # Parabolic Interpolation extrategy\n",
    "    tIntMax, yIntMax = parabolicInterpolation(y)\n",
    "    \n",
    "    \n",
    "    \n",
    "    ## Plotting\n",
    "    ax.clear()\n",
    "    plt.title('Interpeak detection estrategies')\n",
    "    ax.grid(True)\n",
    "    ax.grid(xdata=samplingPoints)\n",
    "    \n",
    "    \n",
    "    ax.plot(time, yReal, label='real signal')\n",
    "    yRealMax = np.max(yReal)\n",
    "    \n",
    "    sampledLabel = 'sampled signal. Error:{:.3f}'\\\n",
    "                   .format(np.abs(np.max(y) - yRealMax))\n",
    "    ax.plot(samplingPoints, y, label=sampledLabel, ls='-.',\n",
    "         color='r', marker='x', markersize=6, alpha=.7)\n",
    "\n",
    "    ax.plot(tResample, yResample, ls='-.',\n",
    "                 color='y', marker='x', alpha=.7)\n",
    "\n",
    "    resMaxLabel = 'Resample Peak. Error:{:.3f}'\\\n",
    "                  .format(np.abs(yResMax - yRealMax))\n",
    "    ax.plot(tResMax, yResMax, label= resMaxLabel, \n",
    "            color='y', marker = 'x', markersize=12)\n",
    "\n",
    "    intMaxLabel = 'Interpolation Peak. Error:{:.3f}'\\\n",
    "                  .format(np.abs(yIntMax - yRealMax))\n",
    "    ax.plot(tIntMax, yIntMax, label= intMaxLabel, \n",
    "            marker = 'x', markersize=12)\n",
    "    \n",
    "    \n",
    "    fig.legend()\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "464181c336284d2ab623cbacc726f1a5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>FloatSlider</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "FloatSlider(value=0.0, max=1.0, min=-1.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "offset = wg.FloatSlider()\n",
    "offset.max = 1\n",
    "offset.min = -1\n",
    "offset.step = .1\n",
    "display(offset)\n",
    "fig, ax = plt.subplots()\n",
    "process()\n",
    "\n",
    "def on_value_change(change):\n",
    "    process()\n",
    "    \n",
    "offset.observe(on_value_change, names='value')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}