m-thesis-introduction/simulations/08_find_beacon_phase_delay.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Find Sine delay in Beacon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import signal\n",
    "\n",
    "\n",
    "import os\n",
    "import sys\n",
    "# Append parent directory to import path so lib can be found\n",
    "sys.path.append(os.path.dirname(os.path.abspath(os.getcwd())))\n",
    "from lib.util import *\n",
    "\n",
    "\n",
    "# monkey patch correlation_lags into signal if it does not exist\n",
    "if not hasattr(signal, 'correlation_lags'):\n",
    "    def correlation_lags(in1_len, in2_len, mode='full'):\n",
    "        r\"\"\"\n",
    "        Calculates the lag / displacement indices array for 1D cross-correlation.\n",
    "        Parameters\n",
    "        ----------\n",
    "        in1_size : int\n",
    "            First input size.\n",
    "        in2_size : int\n",
    "            Second input size.\n",
    "        mode : str {'full', 'valid', 'same'}, optional\n",
    "            A string indicating the size of the output.\n",
    "            See the documentation `correlate` for more information.\n",
    "        See Also\n",
    "        --------\n",
    "        correlate : Compute the N-dimensional cross-correlation.\n",
    "        Returns\n",
    "        -------\n",
    "        lags : array\n",
    "            Returns an array containing cross-correlation lag/displacement indices.\n",
    "            Indices can be indexed with the np.argmax of the correlation to return\n",
    "            the lag/displacement.\n",
    "        Notes\n",
    "        -----\n",
    "        Cross-correlation for continuous functions :math:`f` and :math:`g` is\n",
    "        defined as:\n",
    "        .. math::\n",
    "            \\left ( f\\star g \\right )\\left ( \\tau \\right )\n",
    "            \\triangleq \\int_{t_0}^{t_0 +T}\n",
    "            \\overline{f\\left ( t \\right )}g\\left ( t+\\tau \\right )dt\n",
    "        Where :math:`\\tau` is defined as the displacement, also known as the lag.\n",
    "        Cross correlation for discrete functions :math:`f` and :math:`g` is\n",
    "        defined as:\n",
    "        .. math::\n",
    "            \\left ( f\\star g \\right )\\left [ n \\right ]\n",
    "            \\triangleq \\sum_{-\\infty}^{\\infty}\n",
    "            \\overline{f\\left [ m \\right ]}g\\left [ m+n \\right ]\n",
    "        Where :math:`n` is the lag.\n",
    "        Examples\n",
    "        --------\n",
    "        Cross-correlation of a signal with its time-delayed self.\n",
    "        >>> from scipy import signal\n",
    "        >>> from numpy.random import default_rng\n",
    "        >>> rng = default_rng()\n",
    "        >>> x = rng.standard_normal(1000)\n",
    "        >>> y = np.concatenate([rng.standard_normal(100), x])\n",
    "        >>> correlation = signal.correlate(x, y, mode=\"full\")\n",
    "        >>> lags = signal.correlation_lags(x.size, y.size, mode=\"full\")\n",
    "        >>> lag = lags[np.argmax(correlation)]\n",
    "        \"\"\"\n",
    "\n",
    "        # calculate lag ranges in different modes of operation\n",
    "        if mode == \"full\":\n",
    "            # the output is the full discrete linear convolution\n",
    "            # of the inputs. (Default)\n",
    "            lags = np.arange(-in2_len + 1, in1_len)\n",
    "        elif mode == \"same\":\n",
    "            # the output is the same size as `in1`, centered\n",
    "            # with respect to the 'full' output.\n",
    "            # calculate the full output\n",
    "            lags = np.arange(-in2_len + 1, in1_len)\n",
    "            # determine the midpoint in the full output\n",
    "            mid = lags.size // 2\n",
    "            # determine lag_bound to be used with respect\n",
    "            # to the midpoint\n",
    "            lag_bound = in1_len // 2\n",
    "            # calculate lag ranges for even and odd scenarios\n",
    "            if in1_len % 2 == 0:\n",
    "                lags = lags[(mid-lag_bound):(mid+lag_bound)]\n",
    "            else:\n",
    "                lags = lags[(mid-lag_bound):(mid+lag_bound)+1]\n",
    "        elif mode == \"valid\":\n",
    "            # the output consists only of those elements that do not\n",
    "            # rely on the zero-padding. In 'valid' mode, either `in1` or `in2`\n",
    "            # must be at least as large as the other in every dimension.\n",
    "\n",
    "            # the lag_bound will be either negative or positive\n",
    "            # this let's us infer how to present the lag range\n",
    "            lag_bound = in1_len - in2_len\n",
    "            if lag_bound >= 0:\n",
    "                lags = np.arange(lag_bound + 1)\n",
    "            else:\n",
    "                lags = np.arange(lag_bound, 1)\n",
    "        return lags\n",
    "\n",
    "    signal.correlation_lags = correlation_lags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lag_gridsearch(grid, sample_rate, reference, signal_data):\n",
    "    \"\"\"\n",
    "    Return the best time shift found when doing a grid search.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lag_grid - ndarray\n",
    "        The array specifying the grid that is to be searched.\n",
    "    sample_rate - float\n",
    "        Sample rate of signal_data to transform index to time.\n",
    "    signal_data - ndarray\n",
    "        The real signal to find the time shift for.\n",
    "    reference - ndarray\n",
    "        Real signal to use as reference to obtain lag.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    lag : ndarray\n",
    "        The best time shift obtained\n",
    "    err : tuple\n",
    "        Difference to the previous and next time shift from lag, resp.\n",
    "    \"\"\"\n",
    "\n",
    "    assert signal_data.shape >= reference.shape, str(signal_data.shape) + \" \" + str(reference.shape)\n",
    "    \n",
    "    corrs = grid_correlate(grid, reference, signal_data)\n",
    "    \n",
    "    idx = np.argmax(corrs)\n",
    "\n",
    "    lag = grid[idx]/sample_rate\n",
    "    \n",
    "    err_min =  (grid[idx-1]-grid[idx])/(2*sample_rate)\n",
    "    err_plus = (grid[idx+1]-grid[idx])/(2*sample_rate)\n",
    "\n",
    "    return lag, np.array([err_min, err_plus])\n",
    "    \n",
    "\n",
    "def grid_correlate(grid, reference, x):\n",
    "    \"\"\"\n",
    "    Determine correlation between x and reference using grid as \n",
    "    the lags to be used for the correlation.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    grid - ndarray\n",
    "        The array specifying the grid that is to be searched.\n",
    "    x - ndarray\n",
    "        The real signal to find the time shift for.\n",
    "    reference - ndarray\n",
    "        Real signal to use as reference to obtain lag.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    corrs - ndarray\n",
    "        The correlations along grid.\n",
    "    \"\"\"\n",
    "    grid = np.asarray(grid)\n",
    "    x = np.asarray(x)\n",
    "    reference = np.asarray(reference)\n",
    "\n",
    "    assert x.shape >= reference.shape, str(signal_data.shape) + \" \" + str(reference.shape)\n",
    "    \n",
    "    reference = np.pad(reference, (0,len(x)-len(reference)), 'constant', constant_values=0)\n",
    "    \n",
    "    ref_conj = np.conjugate(reference)\n",
    "    \n",
    "    corrs = np.array([np.dot(np.roll(ref_conj, lag), x) for lag in grid], dtype=np.float64)\n",
    "    \n",
    "    return corrs\n",
    "\n",
    "def correlation_grid(grid_size=None, in1_len=None, in2_len = None, end = None, start=None, mode='full'):\n",
    "    \"\"\"\n",
    "    Abuse correlation_lags to determine the endpoints of the grid.\n",
    "    \"\"\"\n",
    "    \n",
    "    if in1_len is not None or in2_len is not None:\n",
    "        if in2_len is None:\n",
    "            in2_len = in1_len\n",
    "        elif in1_len is None:\n",
    "            in1_len = in2_len\n",
    "\n",
    "        lags = signal.correlation_lags(in1_len, in2_len, mode=mode)\n",
    "\n",
    "        max_lag = max(lags)\n",
    "        min_lag = min(lags)\n",
    "    else:\n",
    "        max_lag = np.inf\n",
    "        min_lag = -np.inf\n",
    "\n",
    "    if end is None:\n",
    "        end = max_lag\n",
    "    elif end > max_lag:\n",
    "        raise ValueError(\"Grid end is too high\")\n",
    "\n",
    "    if start is None:\n",
    "        start = min_lag\n",
    "    elif start < min_lag:\n",
    "        raise ValueError(\"Grid start is too low\")\n",
    "        \n",
    "    if grid_size is None:\n",
    "        grid_size = end - start\n",
    "\n",
    "    return np.linspace(start, end, grid_size, dtype=int, endpoint=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def beacon_time_delay(samplerate, ref_beacon, beacon):\n",
    "    grid = correlation_grid(in1_len=len(ref_beacon), in2_len=len(beacon), mode='full')\n",
    "    time_lag, errs = lag_gridsearch(grid, samplerate, ref_beacon, beacon)\n",
    "\n",
    "    return time_lag, errs\n",
    "\n",
    "def beacon_phase_delay(samplerate, f_beacon, ref_beacon, beacon):\n",
    "    time_delay, errs = beacon_time_delay(samplerate, ref_beacon, beacon)\n",
    "\n",
    "    phase = time2phase(time_delay, f_beacon)\n",
    "    phase_err = time2phase(errs, f_beacon)\n",
    "    \n",
    "    return phase, phase_err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "us = 1e3 # ns\n",
    "ns = 1/us # us\n",
    "\n",
    "samplerate = 500 # MHz\n",
    "timelength = 0.2 # us\n",
    "\n",
    "\n",
    "beacon_time = np.arange(0, timelength, 1/samplerate)\n",
    "f_beacon = 100 # MHz\n",
    "\n",
    "\n",
    "beacon_init_phase = phase2time(42/ns, f_beacon)\n",
    "beacon_phase_offset = 1.2*np.pi\n",
    "\n",
    "beacons = np.array([\n",
    "    sin_delay(f_beacon, beacon_time, phase=-beacon_init_phase),\n",
    "    sin_delay(f_beacon, beacon_time, phase=-beacon_init_phase-beacon_phase_offset)\n",
    "])\n",
    "\n",
    "ref_time = np.arange(0, 1/f_beacon, 1/samplerate)\n",
    "ref_beacon = sin_delay(f_beacon, ref_time, t_delay=0, phase=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reference Delay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cheat Time Delay [ns]: 6.387242824454661 \\pm [0. 0.]\n",
      "Calculated Time Delay [ns]: 5.999999999999994 \\pm [-1.  1.]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAEWCAYAAABmAMpDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd3gU1frHvycBktBD76HaEAhNismKWOjNhoJib9dyr9efV1G5oBd7AVEQA3KVoCKKIiV0EhABadIJJLSQEEJIIAmQQMr7++PdQyab3c3MTtmFO5/n2QeyM3Pms2Vm55x53/MKIoKNjY2NjY2NjY2NjY2NTaAS5G8BGxsbGxsbGxsbGxsbGxtv2B1XGxsbGxsbGxsbGxsbm4DG7rja2NjY2NjY2NjY2NjYBDR2x9XGxsbGxsbGxsbGxsYmoLE7rjY2NjY2NjY2NjY2NjYBjd1xtbGxsbGxsbGxsbGxsQlo7I6rjY2NjY1fEUKMFkKs8MN+Jwgh5li9X29ocRJCJAghnjDB4VohxF9CiDwhxItGt28WQohHhBDr/e1hY2NjY2MOdsfVxsbG5ipCCHFUCJHv7HScFUJsEEI8I4RQdb4XQvQRQqSa7amEiL4jojut3KeNV/4FIIGIahDRFH/L2NjY2NjYAHbH1cbGxuZqZAgR1QAQAeB9AK8C+NqKHQshKlmxHxtTiQCw15cN7c/fxsbGxsYs7I6rjY2NzVUKEeUQ0UIAIwE8LIS4EQCEECFCiI+FEClCiAwhxHQhRJgQohqApQCaCCHOOR9NhBBBQojXhBCHhBBZQoh5Qog6zrZaCiFICPG4ECIFwBrFc48KIY4LIc447/p2F0Lsct4J/kJ6uoZ4Ord9RgiR5Nx2qhBCKJY/JoTY71y2XAgR4e71KzyeEkKcEEKkCyFedlmtihBitvMO9V4hRDfF9vI15wkh9gkhRiiWtRVCrBVC5AghTgshflQsu04IsVIIkS2EOCCEuM/TZySEaOVsJ08IsRJAPZflPZ13zc8KIXYKIfp4aKeNEGKN8/M5LYT4TghR27nsFSHEfJf1PxdCTHbTzhoAtwL4wvn5XyOEqOV8jzKFEMeEEG/KO/jOz+4PIcQkIUQ2gAlu2gwRQkx2fgYnnP8PcS7rI4RIFUK8LIQ45fyMHnXZttx31dP7yZuIz52fS6IQ4jbFglpCiK+d+0gTQkwUQgRX9P45lzcXQvzifA+y5PdX8LHxpvN9OeV8n2o5l8nv38NO/9NCiDe8uNvY2NjYeMHuuNrY2Nhc5RDRZgCpAKKdT30A4BoAkQDaAmgK4N9EdB7AAAAniKi683ECwIsAhgO4BUATAGcATHXZzS0ArgfQT/FcDwDtwB3nyQDeAHA7gPYA7hNC3OJFezCA7gA6AbhPtiuEGA7gdQB3AagP4HcAP1TwFtzq9LgTwGtCiNsVy4YCmAugNoCFAL5QLDsEfs9qAXgLwBwhRGPnsv8AWAEgHEAzAJ87/aoBWAngewANADwAYJoQor0Ht+8BbAN3WP8D4GG5QAjRFMASABMB1AHwfwDmCyHqu2lHAHgP/PlcD6A5SjuRcwD0V3RkK4E/k1jXRoioL/g9fd75+R90vrZaAFqDP+cxAB5VbNYDwGHn633HjdsbAHqCv2+dANwE4E3F8kbO9psCeBzAVCFEuHOZ2++qm324utQDMB7AL8I5yALgWwBFznY6g78PMkfY4/vn7NwuBnAMQEunw1zndo84H7c635/qKPsdAoAoANcCuA3Av4UQ13vxt7GxsbHxBBHZD/thP+yH/bhKHgCOArjdzfObwB0IAeA8gDaKZb0AHHH+vw+AVJdt9wO4TfF3YwCFACqBL+QJQGvFcvlcU8VzWQBGKv6eD+Afzv8/AmC9YhkBiFL8PQ/Aa87/LwXwuGJZEIALACLcvGbpcZ3iuQ8BfO38/wQAqxTLbgCQ7+W93QFgmPP/swHEAGjmss5IAL+7PPcVgPFu2msB7khVUzz3PYA5zv+/CiDWZZvlAB52/j8BwBMeXIcD+Evx91IATzr/PxjAPi+v83K7AIIBXARwg2L50+AcWPnZpVTwnTwEYKDi734Ajiq+b/kAKimWnwJ3dL1+V93s5xEAJwAIxXObATwEoKHzdYQplj0AIL6i98+5z0ylo2K91QD+pvj7WpQ/Npq5+Nxf0XFsP+yH/bAf9qP8w85FsbGxsfnfoCmAbPBdyqoAtglF9C24g+KJCAC/CiFKFM8VgzsDkuNutstQ/D/fzd/VvezzpOL/FxTrRgD4TAjxiWK5AL++Yx7aUrodA9DBy35ChRCViKhICDEGwD/BHRA4HWQo77/Ad0g3CyHOAPiEiGY5/XoIIc4q2q0EN3c34bx7TXynW+nXXPFa7xVCDFEsrwwg3rUhIUQDAFPAd4hrgDv0ZxSrfAvgWQAzADzowccd9QBUQdn39hj4/Za4++yVNHGzfRPF31lEVKT4W37evnxX04iI3OwrAvzepSvaCpLuFbx/zQEcc3H09toqoeyx4em7bGNjY2OjATtU2MbGxuYqRwjRHdzRWA/gNLjT2J6IajsftYhIXkyTmyaOAxigWL82EYUSUZpiHXfbmcFxAE+7uIQR0QYv2zRX/L8F+K6cVwTnzc4A8DyAukRUG8AecMcJRHSSiJ4koibgO5DThBBtnX5rXfyqE9GzbnaTDiDcGV6s9FO+1liXtqoR0ftu2noP/Bl0JKKa4M6pUCxfAKCj4DznwQC+q+g9cHIafAdRmUfcAoCWz/6Em+0r/AxQ8XfVHU2Fomeq2Ndx8B3Xeoq2ahKRDOH29v4dB9BCuJ94yt1rK0LZQRobGxsbGwOwO642NjY2VylCiJpCiMHgfLw5RLSbiErAHbJJzrtMEEI0FULI3NQMAHXlBDNOpgN4x9mZgxCivhBimHWvpAzTAYyVOaPOCXfurWCbcUKIqs5tHgXwYwXrA0A1cEcm07mfRwHcKBcKIe4VQjRz/nnGuW4xOBfyGiHEQ0KIys5Hd3d5jUR0DMBWAG8JIaoIIaIAKO+uzgEwRAjRTwgRLIQIdU5m1My1LfBdwnMAzjpzY19x2VcBgJ/BocibiShFxXsAIioGh2q/I4So4fwO/NPpppYfALzp/N7UA+eoVri9iu+qOxoAeNH5vt8LzleNI6J0cE7yJ87jIsg5IZPMs/b2/m0GDzK8L4So5vwcbla8tpcET7JVHcC7AH70cHfWxsbGxkYHdsfVxsbG5upjkRAiD3yn6A0An6LsZDqvAkgGsEkIkQtgFTg3D0SUCL4YPyx4JtsmAD4DT1y0wtnuJvAkOJZDRL+CJ+yZ63TfA55Qyhtrwa93NYCPiWiFiv3sA/AJgI3gznwHAH8oVukO4E8hxDnwe/N3IjpCRHngSX/uB9+NO+n0DfGwq1Hg9zIbPJnQbIXDcQDDwJNRZYI/z1fg/rf7LQBdAOSAJ3T6xc063zpfh9owYckL4FzTw+C79t8DmKVh+4ngDvouALsBbHc+pwaP31UP/AmeiOs0eKKoe4goy7lsDDjseR94sOFncL424OX9c3beh4AndUoBT3Q20rl4Fvj9XAfgCIAC8PtlY2NjY2MwomwqiI2NjY2NzdWBEKIluDNR2b4DBgghWgBIBNCIiHL97WNjY2NjY6MF+46rjY2NjY3NVY7guqv/BDDX7rTa2NjY2FyJ2LMK29jY2NjYXMU4J3/KAM9429/POjY2NjY2Nj5hhwrb2NjY2NjY2NjY2NjYBDR+DRUWQswSQpwSQuzxp4eNjY2NjY2NjY2NjY1N4OLXO65CCAd4+vnZRHRjRevXq1ePWrZsabqXjY2NjY2NjY2NjY2NjfVs27btNBHVd33erzmuRLTOOeujKlq2bImtW7eaJ2RjY2NjY2NjY2NjY2PjN4QQx9w9H/CzCgshnhJCbBVCbM3MzPS3ji7m7Z2H1NxUf2vY2NjY2NjY2NjY2NhcUQR8x5WIYoioGxF1q1+/3B3jK4ajZ49i5M8jMfj7wcgvzPe3jo2NjY2NjY2NjY2NzRVDwHdcrxaWJi0FAOzM2Im/L/u7n21sbGxsbGx
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "show_originals = True\n",
    "\n",
    "\n",
    "fig, axes = plt.subplots(3,1, sharex=True, figsize=(16,4))\n",
    "axes[0].set_title(\"Determine phase delay for one beacon\")\n",
    "axes[-1].set_xlabel(\"Time [ns]\")\n",
    "\n",
    "## Originals\n",
    "i=0\n",
    "axes[i].set_ylabel(\"Originals\")\n",
    "axes[i].plot(beacon_time/ns, \n",
    "             np.pad(ref_beacon, (0,len(beacons[0])-len(ref_beacon)), 'constant', constant_values=0), \n",
    "             label=\"Reference\",\n",
    "             color='g')\n",
    "axes[i].plot(beacon_time/ns, \n",
    "             beacons[0], \n",
    "             label=\"Beacon 0\",\n",
    "             color='r'\n",
    "            )\n",
    "axes[i].legend(loc='center right')\n",
    "\n",
    "for i in range(1, 3):    \n",
    "    if i == 1:\n",
    "        label = \"Cheat\"\n",
    "        time_delay = phase2time(phase_modulo(beacon_init_phase,0), f_beacon)\n",
    "        time_delay_err = np.array([0,0])\n",
    "\n",
    "        \n",
    "    ## Calculated\n",
    "    else:\n",
    "        label = \"Calculated\"\n",
    "        phase_delay, phase_err = beacon_phase_delay(samplerate, f_beacon, ref_beacon, beacons[0])\n",
    "        time_delay = phase2time(phase_modulo(phase_delay, 0), f_beacon)\n",
    "        time_delay_err = phase2time(phase_err, f_beacon)\n",
    "    \n",
    "    axes[i].set_ylabel(label)\n",
    "    print(\"{} Time Delay [ns]: {} \\pm {}\".format(label, time_delay/ns, time_delay_err/ns))\n",
    "\n",
    "    ### Show original\n",
    "    if show_originals:\n",
    "        axes[i].plot(beacon_time/ns, \n",
    "                     beacons[0],\n",
    "                     color='r', alpha=0.4\n",
    "                    )\n",
    "\n",
    "    axes[i].axvline(time_delay/ns)\n",
    "    axes[i].plot(beacon_time/ns,\n",
    "                 sin_delay(f_beacon, beacon_time, t_delay=time_delay),\n",
    "                 ls = (4, (10, 4)),\n",
    "                 label='sin')\n",
    "    axes[i].plot(ref_time/ns,\n",
    "                 time_roll(ref_beacon, samplerate, time_delay, sample_shift=0),\n",
    "                 ls = (6, (10, 4)),\n",
    "                 label='roll')\n",
    "    axes[i].legend(loc='center right')\n",
    "\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine delay between beacons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAEWCAYAAABmAMpDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydZ3gc1dmGn1e9d8m2ZKvLkg2EVEIvgdBJIOELIUAIAVIIqaQQEhLRSwgBAiGE0AIEQg29g+kETMc7KrasZlmWVr23fb8fZ85qtNoyM2dWXsM817WXrZ2d2XunnzL3IWaGGzdu3Lhx48aNGzdu3LhxE6uJ294Abty4cePGjRs3bty4cePGTbi4BVc3bty4cePGjRs3bty4cRPTcQuubty4cePGjRs3bty4ceMmpuMWXN24cePGjRs3bty4cePGTUzHLbi6cePGjRs3bty4cePGjZuYjltwdePGjRs3bty4cePGjRs3MR234OrGjZuohIjOIaJ/bm8OMyGieiK6IwrLvZWILnR6uW52zBDR/kTkI6JRIjpUcVkH6cvxEdFBTjG6ceNkiIiJaIyILtreLIEhotX6MTRHRKdtbx43btxEjltwdePGja3oF/xRw83zhOHvE5j5YmaO+s0AEZXrN0fyu7cR0aNE9OVof3c0syMUeolo3VLc8BlufkeJyEtEdxFRTrS/N0rpYuYMZn5SZSHM/CwzZwBoNzuPvk9NBxy78SE+m0xEfyGiLiIaIKK/EVGiYfoaInqeiIaIaCMRHRMw/2n6+6NE9CQRFRum5RDRbUTUo7/qA+bdk4jeJKIRIvqAiPY2TCMi+h0RtRPRMBHdTURZhuklRPQQEfUTUScR/SBg2UcR0Uc612tEtPZj8JvziOg/+rHhJaI7jdMNn1tJRB2B7y9BdmXm34WaSETf0LfFOBGtU/kiIjqAiF7Qt1FruM8yc5N+DL2s8p1u3LhZurgFVzdu3NiKfvOdYbh5Psrw3p3bASlHZ9kVwDMAHiSi72wHDjfRya769q0EkAugfvvi7LC53HjsMvNciM+dDeDzAHYGsBrAZwH8HgCIKAHAQwAeBZAH4HsA7iCi1fr0/QBcDOCr+vTNAO4yLPsvANIAlAPYDcBJRHSKPm8egIcB/AlADoDLATxCRLn6vN8GcBKAvQAUA0gF8FfDsu/Qv28ZgCMAXExEB+jLrgFwJ4Af6Mt+BMDD+u/ZkX/zhRDHRCWAKv231yMgzNwJoJ+Idg2ctp3TD+AqAJc6sKwxADcD+JUDy3Ljxk2MxS24unHjJiohQ/dbmm8VPYWIOvTWjB8Q0Rf01oVBIro2YP7vEpGmf/YpIioz873M3M3MV0PcuF1GRHH68oqJ6H4i6iWizUT0kzDs9xJRt15r/xIR7aS//wUSLboJhs9+nYjeC4NUQETP6C0pLxp/BxHV6dP6iaiRiL6hv/89ACcA+LXeevOIvu4eMcy7kYjuMfzdQUSfDrdcfVoyEV2ht95sI6K/E1GqPm1/Eq1UZ+mtQlvlzXWQdXQRgH0AXKszXktE5xHRX/XpiSRaSS/X/04lokl5M05EXyGiDfq2X0dEa8KsQ3+YeRjiJt/YUpZNRDfpvFuI6ELSWxKJqIpEK1mfoTUqxzDvKiJ6QN8v+uR+SERxRPR7ImrT18W/iChbnyb355P19eglot8ZlrkbEa0n0Tq2jYiuDPV7iOg7RPSKvk0G9H3zsIDpLfr+s5mITjCznhzIUQCuYeZ+Zu4FcA2A7+rT6iAKUH9h5jlmfh7AqxCFKznvvcy8gZmnAVwAYF8iqjJMv5yZx5m5FcBNhmXvCWAbM9+rL/sOAL0AvmaY9yZm7mDmUQCXATiOiNKIKAPA/gAuYuYZZn4fwH2GZR8C4GVmfoWZZ/V5SwDst6P+Zn16BYD/MvMwMw8BeBDATgiexwEcLv8w7MvfI9HSvJWIzjJMr9fPh3fo++CHJLrY/lY/LjqI6OAQ32Uqeu+BewB0qSxHX9abzHw7gBbVZblx4yb24hZc3bhxs5T5IoAaAMdB1LD/DsBBEDdZ3yDRagEiOhrAORA3boUQXbnuCrbAMHkAQBGAWhKF10cAvA9xo3oggJ8R0SEh5n1C5ywC8A5EKw2Y+S0AfQCM3ZBPBHB7GI4TIG5iCwC8J5dFROkQLcP/1r/neAB/I6KdmPkf+udk69hRAF4EsI9eoFoBIBGiBQZEVAkgA8AH4Zar81wG0Zr0aQDV+vr4g4F3OYBs/f1TAVxH8y0//uhd/14GcKbOeKbOuL/+kS8A6MZ8oWAPAI3MPECileouAD+D2L6PQ7QwJYVZj9B/ay6AowG8YXj7NgCz+u/5DICDAcguzATgEohCxxoAq6C3RumF20cBtEG0hJUAuFuf7zv66wCIlqwMAAsqVwDsDaAWYn/6g6HwfTWAq5k5C6IF7B6EzxcBNELsI5cDuIlE0iEKT4cxcyZEASdkJQkRfYuIPojwXWfoFRpvE9HXw3yO9Jfx75V64Z1CfH7nMPPCMB1Bpoea18yykyGOVzK8Z3beSNNj/TcDwHUAjiSiXP34+DrEOSxYHoOh4GrIAfryDgZwNi18bvooiHNcLoB3ATwFcf9YAuB8ADeE+C43bty4cTbM7L7cl/tyX0ovAK0ADgp4rx7AHfr/ywEwgBLD9D4Axxn+vh/Az/T/PwHgVMO0OADjAMqCfLdcdkLA+yn6+3tBFAzaA6b/FsAtgaxBlp+jLydb//s3AO7U/5+nc60IMe+tAO42/J0BYA6i8HQcROuP8fM3APijYd4LA6Z3QHRf/CaAfwB4E6Il6BQAD+ufCblciBveMQBVhml7ANis/39/ABPGdQmgB8DuIX7fOgCnGf5OBTAJIB+i2+U5ADr1330eRGsWAJwL4J6A7bsFwP4hvocBDAMY1Ndfg9yXILpFTgFINXz+eAAvhFjW0QDeNfz23sB9R5/2HIAzDH/XApgBkGDY51Yapr8J4Jv6/1/Sf29BwDL3B9AZ8N53AGw0/J2mL3s5gHT9N3/d+PsiHXthjtPP6tsmAaLwMgJgrxCfvRCiRbFQZ/mfziUrTVoA/Fr//8EApgE8pc97IAAvgE/p+8QNAHwAjten3wFRsZQJUdmwCcCUPi1f/83H68s+WZ/3Bn36aQCa9G2QDdH6zgD20Ke/AtGNNkX/vf0QFSaAOFbG9O2QBLEf+gD8dgf/zcUAntXn8UFUXCWF2K7xEMd0bsD5s87wmcshWngBcW58xjDtKACjAOL1vzP1+XPCHLvVJvfP0wCsM/NZE8s6CECryc+ug+E85r7cl/uK3Zfb4urGjZulzDbD/yeC/J2h/78MwNV6N9JBiJtPgqjhNxv52X59ecVyefoyz4Eo9CwIEcUT0aVEtImIhiEKBoBoDQPEDehRerfEb0AUEreG4fDLUFh08+uHuNEsA/DFAKYTIG6YQ0W2aO6r/38dRIvmfvrfiLDcQoiC0duGaU/q78v0sehGKTOO+e0SNsw8AWC9ziMZX4OoPDAyFkO0csr5fBDrKdz2/Swz50AUSK4H8DIRpei/NxHAVsNvugGitRlEVERCZrNF3553YH5brgLQFvB7ZRYw6v9PwMJ9ptvwf+N6OhWiVbuBiN4ioiPD/K4Fy2Hmcf2/Gcw8BlER8QP99z1GRHURlhUyzPwOM/cx8ywzPw7Rqv+1EB+/CKJ17T2IbfhfiIJ7DzPPQFQAHKGznwXRqtypf89zEBUl90Ost1aIQnKnvuyfQBzvzRDPjd5lmLcP4jnRX0CcHw6FKJTJeW/WP78OwAYAL+jvy+knQHSd7YDYT+40LLsBolB4LYCtEPuBxzDvjvqb74Uo2GYCyIIoFAe1pLN4pvl5iG7TxhilTW0Q+79M4Hnay/PPRk/o/5o6R5B4NEGKwc4xM48bN27cyLgFVzdu3MRiOgB8n5lzDK9UZn7NwjKOgWhZaNSXtzlgeZnMHKzL3LcgbiIPgmjdKNffJwBg5i0AXteXfxLCdxMGROFILEAUdvMgnuXqAPBiAFMGM/9Q/zgHWZYsuO6j//9FLC64hluuF+JGcyfDtGwW0iM
      "text/plain": [
       "<Figure size 1152x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "abs_time_delay, time_delay_err = beacon_time_delay(samplerate, beacons[0], beacons[1])\n",
    "time_delay = abs_time_delay % (1/f_beacon)\n",
    "\n",
    "\n",
    "fig, axes = plt.subplots(2, sharex=True, figsize=(16,4))\n",
    "axes[-1].set_xlabel(\"Time [ns]\")\n",
    "axes[0].set_title(\"Time Delay between two Beacons[ns]: {} \\pm {}\".format(time_delay/ns, time_delay_err/ns))\n",
    "\n",
    "\n",
    "## Originals\n",
    "i=0\n",
    "axes[i].set_ylabel(\"Originals\")\n",
    "axes[i].plot(beacon_time/ns, \n",
    "             beacons[0], \n",
    "             label=\"Beacon 0\",\n",
    "             color='g')\n",
    "axes[i].plot(beacon_time/ns, \n",
    "             beacons[1], \n",
    "             label=\"Beacon 1\",\n",
    "             color='r'\n",
    "            )\n",
    "axes[i].axvline(0, color='g', alpha=0.3)\n",
    "axes[i].axvline(time_delay/ns, color='r', alpha=0.3)\n",
    "\n",
    "## Correct Beacon 1\n",
    "i=1\n",
    "axes[i].set_ylabel(\"Corrected Beacon 1\")\n",
    "axes[i].plot(beacon_time/ns, \n",
    "             beacons[0], \n",
    "             #label=\"Beacon 0\",\n",
    "             color='g',\n",
    "            )\n",
    "axes[i].plot(beacon_time/ns,\n",
    "             time_roll(beacons[1], samplerate, -time_delay, sample_shift=0),\n",
    "             ls = (6, (10, 4)),\n",
    "             color='orange',\n",
    "             label='Corrected Beacon 1'\n",
    "            )\n",
    "\n",
    "# put legend between plots\n",
    "lines_labels = [ax.get_legend_handles_labels() for ax in fig.axes]\n",
    "lines, labels = [sum(lol, []) for lol in zip(*lines_labels)]\n",
    "\n",
    "fig.subplots_adjust(hspace=0.3)\n",
    "fig.legend(lines, labels ,loc='right', ncol=3, bbox_to_anchor=(0.55, 0.51))\n",
    "fig.show()"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}