2022-06-14 17:03:00 +02:00
|
|
|
{
|
|
|
|
|
"cells": [
|
|
|
|
|
{
|
|
|
|
|
"cell_type": "markdown",
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"source": [
|
|
|
|
|
"# Perform a grid search to find the timing error"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"cell_type": "code",
|
|
|
|
|
"execution_count": 1,
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"outputs": [],
|
|
|
|
|
"source": [
|
|
|
|
|
"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",
|
|
|
|
|
"\n",
|
|
|
|
|
"import numpy as np\n",
|
|
|
|
|
"rng = np.random.default_rng(12345)\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"import matplotlib.pyplot as plt\n",
|
2022-07-11 11:20:03 +02:00
|
|
|
"from scipy import signal\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
"import copy\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"from lib.signals import DigitisedSignal\n",
|
|
|
|
|
"import lib.util as util\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"km = 1e3\n",
|
|
|
|
|
"ns = 1e9\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"# monkey patch correlation_lags into signal if it does not exist\n",
|
2022-07-11 11:20:03 +02:00
|
|
|
"if not hasattr(signal, 'correlation_lags'):\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
" 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",
|
2022-07-11 11:20:03 +02:00
|
|
|
" signal.correlation_lags = correlation_lags"
|
2022-06-14 17:03:00 +02:00
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"cell_type": "code",
|
|
|
|
|
"execution_count": 2,
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"outputs": [
|
|
|
|
|
{
|
|
|
|
|
"data": {
|
|
|
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACVuklEQVR4nOyddZhc1d2A3zPu6xt3IQLxBHd3t1IopYWihQqFflSoQFugghWKFSheiru7xggRQhJim2yS1XGf8/1x7/id2dnI7rK57/PsMztXzpx7Z/f87s+FlBIdHR0dHZ1SGHp7Ajo6Ojo6fRtdUOjo6OjolEUXFDo6Ojo6ZdEFhY6Ojo5OWXRBoaOjo6NTFl1Q6Ojo6OiURRcUOjo6Ojpl0QWFTr9ECLFGCHFIL8/BKoS4VwixVgjhF0IsEEIcmbN/pBBCCiECOT+/LjPeSCHES0KIDiHEJiHEbUIIU89cjc7OjP5HpqOz4zAB64H9gXXAUcATQojdpJRrco6rllImKhjvn8AWYBBQDbwOXAzcsh3nrKNThK5R6OxUCCFqhBAvCCFa1CfzF4QQQ3P2jxJCvKdqAG8IIW4XQjy0NZ8lpQxKKa+VUq6RUqaklC8Aq4GZWzn9UcATUsqIlHIT8AowWetAIcS5QogPhBA3qde5ukCbOVcI8Y16nauFEGdt5Zx0dgJ0QaGzs2EA/g2MAIYDYeC2nP2PAJ8BdcC1wNm5JwshFgkhvrM1HyyEGACMB5YU7ForhGgSQvxbCFFfZoibgTOEEA4hxBDgSBRhUYrdgeVAPXADcK9QcKJoIUdKKd3AXsDCrbkmnZ0DXVDo7FRIKduklP+TUoaklH7gOhTTEEKI4cBs4DdSypiU8gPguYLzp0gpH+nu5wohzMDDwANSyq/Uza3q541A0TLc6jGleBdFg/ABTcBc4Jkyx6+VUt4tpUwCD6CYrAao+1LArkIIu5SyWUpZKLx0dDLogkJnp0J9Gv+X6mD2Ae8B1UIIIzAYaJdShnJOWd+NsV/OcUqflbPdAPwHiAGXprdLKQNSyrlSyoSUcrO67zAhhEdjbAPwKvAU4ETREmqAv5SZ0qacz0pfk0tKGQROBy4EmoUQLwohJlR6nTo7H7qg0NnZ+BmwC7C7lNID7KduF0AzUCuEcOQcP6zSgaWUR0opXerPwwBCCAHci/Ikf7KUMl5uiJy5FFKrzuU2KWVUStmGYkI7qtL5Fcz1VSnloShaxlfA3Vszjs7OgS4odPozZiGELefHhGLeCQOdQoha4Lfpg6WUa1HMOdcKISxCiD2BY7dxDncAE4FjpZTh3B1CiN2FELsIIQxCiDoUv8E7Ukpv4SBSylYUR/hFQgiTEKIa+B7wRXcnJIQYIIQ4TvVVRIEAkOzuODo7D7qg0OnPvIQiFNI/1wL/AOwo/oFPKHYGnwXsCbQBfwQeR1lMARBCLKk0QkgIMQL4ETAN2KRhlhqtfr4fWKx+zpk55/+fEOLlnCFPAo4AWoCVQAL4SSVzKcCAolltBNpRfDQXb8U4OjsJQm9cpKNTGiHE48BXUsrfdnmwjk4/RdcodHRyEELMFkKMUc1BRwDHUz6ySEen36NnZuvo5DMQJbKoDiUE9SIp5YLenZKOTu+im550dHR0dMqim550dHR0dMrSL01P9fX1cuTIkb09DR0dHZ1vDfPmzWuVUjZo7euXgmLkyJHMnTu3t6eho6Oj861BCLG21D7d9KSjo6OjUxZdUOjo6OjolEUXFDo6Ojo6ZemXPgodHZ2dk3g8TlNTE5FIpLen0mex2WwMHToUs9lc8Tm6oNDR0ek3NDU14Xa7GTlyJErhXp1cpJS0tbXR1NTEqFGjKj6v10xPQohhQoi3hRDL1EJrl2scI4QQtwghVqqdxWb0xlx1dHS+HUQiEerq6nQhUQIhBHV1dd3WuHpTo0gAP5NSzhdCuIF5QojXpZRLc445Ehin/uyOUrJ5956fqo6OzrcFXUiUZ2vuT69pFGr7xfnq735gGTCk4LDjgQelwiconcgG9fBUdbTwb4Z590Nq+7Qx+HhVGx+vatsuY+no6Gxf+kTUkxBiJDAd+LRg1xDyW1E2USxM0mNcIISYK4SY29LSskPmqZPDMxfB85fD0me3eShvKM6Zd3/CmXd/Qjim98/R+XZjNBqZNm0au+66K8ceeyydnZ3b/TMOOOCAHk0q7nVBIYRwAf8DrpBS+gp3a5yiWcVQSnmXlHKWlHJWQ4NmFrrO9iKVhLUfKb+venObh1t13yPMvfUs7nnyd3z94tvbPJ6OTm9it9tZuHAhixcvpra2lttvv723p0QyuW0PYL0qKIQQZhQh8bCU8imNQ5rI71k8FKUrl05v4m2ChNrVc/OSbR4utnAR9SEvh6z6nPCbuqDQ6T/sueeebNiwAYBVq1ZxxBFHMHPmTPbdd1+++uqrzPY99tiD2bNn85vf/AaXywXAO++8wzHHHJMZ69JLL+X+++8v+oyLLrqIWbNmMXnyZH7722x/rZEjR/L73/+effbZh//+97/bdB295szOaTq/TEr5txKHPQdcKoR4DMWJ7ZVSNvfUHHVK0L5KeW2YCO2rt3k4w6pVtLhrscZjiFWrtnk8HR2A3z2/hKUbC40U28akwR5+e+zkio5NJpO8+eab/OAHPwDgggsu4M4772TcuHF8+umnXHzxxbz11ltcfvnlXH755Zx55pnceeed3Z7TddddR21tLclkkoMPPphFixYxZcoUQMmZ+OCDD7o9ZiG9qVHsDZwNHCSEWKj+HCWEuFAIcaF6zEvANyj9ge9G7+vbN2j/RnkdtR9EOiEa2KbhXE1raBkwjNbaAVibN2z7/HR0epFwOMy0adOoq6ujvb2dQw89lEAgwEcffcSpp57KtGnT+NGPfkRzs/LM+/HHH3PqqacC8J3vfKfbn/fEE08wY8YMpk+fzpIlS1i6NBs4evrpp2+Xa+o1jUJK+QHaPojcYyRwSc/MSKdifBtBGGHITPX9BmjYZauHq27fwsaJ04gHfLhbNm2nSers7FT65L+9SfsovF4vxxxzDLfffjvnnnsu1dXVLFy4sOJxTCYTqVQq814r92H16tXcdNNNfP7559TU1HDuuefmHed0OrfpWtL0ujNb51tIsAWc9VA9XHnvXV/++C74fMgkWnabSXDoCIjFtsMEdXR6n6qqKm655RZuuukm7HY7o0aNyvgKpJR88cUXAOyxxx7873//A+Cxxx7LnD9ixAiWLl1KNBrF6/Xy5pvFgSM+nw+n00lVVRWbN2/m5Zdf3iHXogsKne4TaAFnI7gHZN9vJaFYgsuP+glrTj+XuT+5loN/8E9iiVTXJ37bCHfCw6fB8le6f+7GhfDgCVmTH/DBilbOue8z2oPdF6yfXPJ/fHbwSchUP7zPfYzp06czdepUHnvsMR5++GHuvfdepk6dyuTJk3n2WSW0/B//+Ad/+9vfmDNnDs3NzVRVVQEwbNgwTjvtNKZMmcJZZ53F9OnTi8afOnUq06dPZ/LkyZx33nnsvffeO+Q69FpPOt0nuEXRKBx1yvvQ1ifKtfqVha7eZSGiCojOUIxGj22bp9mnWPgIrHgVWr+GXY7o3rnv3QjfvA2f3QNHXA/AnQ+9wwcRG09/tJIfHDqp4qFCHV72+OefAPj6pXcYf8xB3ZuLTpcEAvk+u+effz7z+yuvFD8oDBkyhE8++QQhBI899hizZs3K7Lvhhhu44YYbis555513Mr9rRUIBrFmzpnsTL4MuKHS6T7AF6saC1QMG0zYJitD7H7Dg5jNZM+lBAhY7dz31R/wnDKFxj2nbb759gXUfK6+dayEWBEuFtmMpoUlNrGr6DICIL8D9fzgdk0yx4M0D4NDKQ4qXfbWe1JBJzN6wlPY33gVdUPQ68+bN49JLL0VKSXV1Nffdd19vT6kI3fSk0z2kVE1PDSCEolVsg6AIr99ITcSPq76GmmSEw1Z8Qnjlmu03375C2mwkU9CyvPLzgi0QUB38W5aBlGxasBSTVLSvhjUrujWNBQkHp571F9rsHgxLtj0HRmfb2Xffffniiy9YtGgR7733HmPHju3tKRW
|
|
|
|
|
"text/plain": [
|
|
|
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
"metadata": {
|
|
|
|
|
"needs_background": "light"
|
|
|
|
|
},
|
|
|
|
|
"output_type": "display_data"
|
|
|
|
|
}
|
|
|
|
|
],
|
|
|
|
|
"source": [
|
|
|
|
|
"# Setup signal\n",
|
|
|
|
|
"########################\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"signal_t_length = 200e-9\n",
|
|
|
|
|
"signal_sample_rate = 5e9 # Hz\n",
|
|
|
|
|
"signal_time = util.sampled_time(\n",
|
|
|
|
|
" sample_rate=signal_sample_rate,\n",
|
|
|
|
|
" start=0,\n",
|
|
|
|
|
" end=signal_t_length\n",
|
|
|
|
|
")\n",
|
|
|
|
|
"double_signal_time = np.concatenate((signal_time, signal_time[-1] + signal_time)) \\\n",
|
|
|
|
|
" - (signal_time[-1] - signal_time[0])/2\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"fc = 50e6# MHz\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"if True:\n",
|
|
|
|
|
"\n",
|
|
|
|
|
" shift = -25.8e-9\n",
|
|
|
|
|
" signal_pulse1_center = 125e-9 # 100 ns\n",
|
|
|
|
|
" \n",
|
|
|
|
|
" if not True:\n",
|
2022-07-11 11:20:03 +02:00
|
|
|
" signal_1 = DigitisedSignal(signal.gausspulse(signal_time - signal_pulse1_center, fc=fc),\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
" signal_sample_rate)\n",
|
2022-07-11 11:20:03 +02:00
|
|
|
" signal_2 = DigitisedSignal(signal.gausspulse(signal_time - signal_pulse1_center, fc=fc), \n",
|
2022-06-14 17:03:00 +02:00
|
|
|
" signal_sample_rate)\n",
|
|
|
|
|
" signal_2.t_0 = shift\n",
|
|
|
|
|
" else:\n",
|
|
|
|
|
" signal_1 = DigitisedSignal(np.sin(2*np.pi*(signal_time)*fc/5.4) + np.sin(2*np.pi*(signal_time)*fc),\n",
|
|
|
|
|
" signal_sample_rate)\n",
|
|
|
|
|
"\n",
|
|
|
|
|
" signal_2 = copy.copy(signal_1)\n",
|
|
|
|
|
" signal_2.t_0 = shift\n",
|
|
|
|
|
" \n",
|
|
|
|
|
" signal_1.periodic = False\n",
|
|
|
|
|
" signal_2.periodic = False\n",
|
|
|
|
|
" del signal_pulse1_center\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"del fc, signal_t_length, signal_sample_rate\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"# Pre calculated values\n",
|
|
|
|
|
"sig_1 = signal_1(double_signal_time)\n",
|
|
|
|
|
"sig_2 = signal_2(double_signal_time)\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"if True:\n",
|
|
|
|
|
" fig, ax = plt.subplots()\n",
|
|
|
|
|
" ax.plot(double_signal_time * ns, sig_1, label='Regular')\n",
|
|
|
|
|
" ax.plot((double_signal_time) * ns, sig_2, label=\"Shifted\")\n",
|
|
|
|
|
" ax.plot((double_signal_time) * ns, signal_2(double_signal_time + shift), 'r--', label=\"Cheat\")\n",
|
|
|
|
|
" ax.set_xlabel(\"Time (ns)\")\n",
|
|
|
|
|
" ax.set_ylabel(\"Amplitude\")\n",
|
|
|
|
|
" ax.set_title(\"Lag: {:.4g} ns\".format(shift*ns))\n",
|
|
|
|
|
" plt.legend()\n",
|
|
|
|
|
" plt.show()"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"cell_type": "code",
|
|
|
|
|
"execution_count": 3,
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"outputs": [
|
|
|
|
|
{
|
|
|
|
|
"data": {
|
|
|
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACY40lEQVR4nOyddZgd1d2A33Pd967GPcHiQiBAILi709JSSoHiFKh8VeoFSikFSinuxTVQ3IolgUCMkIQkRDab9et+vj9mrs/cvbHNksz7PPvcvSNnzswm5zc/F1JKDAwMDAwM9DBt7wkYGBgYGPRtDEFhYGBgYFARQ1AYGBgYGFTEEBQGBgYGBhUxBIWBgYGBQUUMQWFgYGBgUBFDUBj0KYQQUggxegvH+I0Q4sGtNacerjVLCLG2N65lYLC9MASFwVZHCHGWEGKuECIkhGgWQrwkhNhve89LCyHEvUKI32/D8bdY8G2FOTQJIR4RQqwXQnQLIf4nhNirYP8sIURG/Xtlf75bYbx9hBAfCyGCQojP++rf1mDrYQgKg62KEOJHwE3AH4F+wFDgNuD47TitnR0PMAeYCtQB9wEvCiE8Bcesl1J6Cn7u0xpICFEHPAdcD/iB64DnhRC12/IGDLYvhqAw2GoIIWqA3wIXSymfklKGpZRJKeXzUspr1GOmCyE+EEJ0qdrGLUIIm854TiHEX4UQq9U34ffUbWXmHiHEKiHEITrjPC6E2KCO8Y4QYqy6/XzgW8CP1bfo59XtA4UQTwohWoUQK4UQl5XM6V4hRKcQYjGw52Y+q1FCiDeEEO1CiDYhxENCCH/B/ilCiE/Vt/bHhRD/2VzNR0r5lZTyRills5QyLaW8A7ABu27GcPsALVLKx9WxHgRagZO0DlbNgI8JIe5X72WREGJawf6fCCHWqfuWCiEO3px7NNi2GILCYGsyA3AAT1c4Jg1cCTSoxx8MXKRz7A0ob8H7oLwJ/xjIbMa8XgLGAE3AJ8BDAOqC+RBwnfoWfawQwgQ8D3wGDFLnd4UQ4nB1rF8Do9SfwwFdE00PCOBPwEBgd2AI8BsAVXA+DdyLct+PACcWnawI2s0y+QghJqEIiuUFm5uEEC2qYPybEMJdYd5CY9u4Cpc8DngURQN5DrhFnceuwCXAnlJKL8rzXLVJN2PQKxiCwmBrUg+0SSlTegdIKedJKT+UUqaklKuAfwEHlB6nLtjnApdLKdepb6/vSynjmzopKeXdUsqgeu5vgImq9qPFnkCjlPK3UsqElPIr4N/AGer+04A/SCk7pJRrgJs3dT7qnJZLKV+VUsallK3AjeSfw96ABbhZ1cieAj4uOd8vpXxvU68rhPABDwDXSim71c1fAJOAAcBBKML5Rp0h3gcGCiHOFEJYVV/GKMBV4bLvSSlnSynT6rUnqtvTgB3YQwhhlVKuklKu2NR7Mtj2GILCYGvSDjQIISx6BwghdhFCvKCaggIovowGjUMbULSTLVo4hBBmIcSfhRAr1OutKhhfi2EoC2FX9gf4PxR/CygawJqC41dv5ryahBCPqmaXAPBgwZwGAutkccXONWWD6I+9qMApPbNguxNFW/pQSvmn7HYp5QYp5WIpZUZKuRJFcztFa2wpZTuKv+lHQAtwBPAaUCnya0PB7xHAIYSwSCmXA1egCO+N6vMYWO19GvQehqAw2Jp8AMSAEyoc80+UN9gxUkofyiJcasoAaFPHGqWxL0zBG6wQwgw06lzvLJSF7RCgBhiePU39LC2fvAZYqb6xZ3+8Usqj1P3NKGaiLEN1rtsTf1KvPUF9Dt8umFMzMEgIUfhchlAlUsqxBU7pdwGEEHbgGWAdcEFPQ6D9N8mO/7aUck8pZR1wNoqv42O943uY68NSyv1QBLQE/rI54xhsWwxBYbDVUE0ZvwJuFUKcIIRwqeaJI4UQ16mHeYEAEBJC7Ab8UGesDHA3cKPqXDYLIWaoC96XKG+lRwshrMAvUEwYWniBOIq240LRYAppAUYWfP8YCKhOVqd63XFCiKzT+jHgZ0KIWiHEYODSKh6NTQjhKPgxq/MKAV1CiEHANQXHf4BilrlECGERQhwPTK/iOpqoz+gJIAp8R322hftnCSGGCoUhwJ+BZyuMN1n9u/pQ/EhrpZT/3Yx57SqEOEj9m8bU+aU3dRyDbY8hKAy2KlLKG1HMEr9AiYZZg+KwfEY95GqUt/wgiu3/PxWGuxpYgBLa2YHytmlSBdJFwJ0ob8hh9E0f96OYh9YBi4EPS/bfhWIj7xJCPKPa0Y9FsdmvRNFs7kTRRgCuVcdbCbyCYnPviUUoi2D253vqOFOAbuBF4KnswVLKBEoU0feBLhRt4wUUgQdAqVmpB/YBjgEOQxFMpWapKSjCKYzig1gIFEZ63S6EuL1gvB+jPJc1KH6NIkf7JmBHEUptKOapJhQN06CPIYzGRQYGfR8hxEfA7VLKe7b3XAx2PgyNwsCgDyKEOEAI0V81PX0XmAC8vL3nZbBzohudYmBgsF3ZFcUf4kGJ/DpFStm8fadksLNimJ4MDAwMDCpimJ4MDAwMDCqyQ5qeGhoa5PDhw7f3NAwMDAy+McybN69NSqmZj7RDCorhw4czd+7c7T0NAwMDg28MQgjdKgOG6cnAwMDAoCKGoDAwMDAwqIghKAwMDAwMKrJD+igMDAy+2SSTSdauXUssFtveU9nhcDgcDB48GKvVWvU5hqAwMDDoc6xduxav18vw4cMpLqJrsCVIKWlvb2ft2rWMGDGi6vO2m+lJCDFECPGmEGKJWj//co1jhBDiZiHEcqE0cZ+yPeZqYGDQu8RiMerr6w0hsZURQlBfX7/Jmtr21ChSwFVSyk+EEF5gnhDiVSnl4oJjjkRpYTkG2Aull8FevT9VAwOD3sYQEtuGzXmu202jUBu9f6L+HgSWoPQoLuR44H6p8CHgF0IM6OWpGmgRaoVPH4TM1mkfMGdVBx+v7NgqYxkYGGxd+kTUkxBiODAZ+Khk1yCKW0CupVyYZMc4XwgxVwgxt7W1dZvM06CAZ34Iz14MS57f4qG6I0lOvf0DTvvXB8SSRt8ag76B2Wxm0qRJjBs3jmOPPZaurq6tfo1Zs2Z9I5KDt7ugEEJ4gCeBK6SUgdLdGqdoVjGUUt4hpZwmpZzW2KjXFdNgq5BJw6r3lN+Xv7bFwy157jWefOBq3rn9+6y8+5EtHs/AYGvgdDqZP38+CxcupK6ujltvvXV7T4l0evu8SG1XQaG2aHwSeEhK+ZTGIWsp7hU8GFjfG3MzqEBgHaSiyu8tC7d4uMjcT5m6/guGdrcQf+/9LR7PwGBrM2PGDNatWwfAihUrOOKII5g6dSozZ87kiy++yG3fe++92XPPPfnVr36Fx+MB4K233uKYY47JjXXJJZdw7733ll3jhz/8IdOmTWPs2LH8+te/zm0fPnw4v/3tb9lvv/14/PHHt+Fd6rPdnNlq4/i7gCVq+0wtnkPpG/woihO726jJ3wdoX6F8Nu4GHSu3eDjzl0uJWWy0e2oxrV61xeMZ7Fhc+/wiFq8vNTZsGXsM9PHrY8dWdWw6neb111/n+9//PgDnn38+t99+O2PGjOGjjz7ioosu4o033uDyyy/n8ssv58wzz+T222/vYdRy/vCHP1BXV0c6nebggw/m888/Z8KECYCS+/Dee+9t8phbi+2pUewLnA0cJISYr/4cJYS4UAhxoXrMbOArYDlKf+WLttNcDQrp+Er5HD4TYl0QD23RcM6vV7GxcRCtDQNxNeu1vjYw6F2i0SiTJk2ivr6ejo4ODj30UEKhEO+//z6nnnoqkyZN4oILLqC5WXl3/eCDDzj11FMBOOusszb5eo899hhTpkxh8uTJLFq0iMWL8wGgp59++ta5qc1ku2kUUsr30PZBFB4jgYt7Z0YGVRNYD8IMg/eEOf9WTFGNu272cK7ONiIN/QjWNDLkc8P0ZFBMtW/+W5usj6K7u5tjjjmGW2+9lXPOOQe/38/8+fOrHsdisZDJZHLftXIYVq5cyQ033MCcOXOora3lnHPOKTrO7XZv0b1sKdvdmW3wDSS8EdwN4FfdR91rKh/fA+02N4F
|
|
|
|
|
"text/plain": [
|
|
|
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
"metadata": {
|
|
|
|
|
"needs_background": "light"
|
|
|
|
|
},
|
|
|
|
|
"output_type": "display_data"
|
|
|
|
|
}
|
|
|
|
|
],
|
|
|
|
|
"source": [
|
|
|
|
|
"## Using scipy\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"if True:\n",
|
2022-07-11 11:20:03 +02:00
|
|
|
" scipy_corr = signal.correlate(sig_1, sig_2, mode='same')\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
" scipy_id = scipy_corr.argmax()\n",
|
|
|
|
|
" calc_shift = ((len(scipy_corr)-1)/2 -scipy_id)/signal_1.sample_rate\n",
|
|
|
|
|
"\n",
|
|
|
|
|
" if not True:\n",
|
|
|
|
|
" signal_2.periodic = True\n",
|
|
|
|
|
"\n",
|
|
|
|
|
" fig, ax = plt.subplots()\n",
|
|
|
|
|
" ax.plot(double_signal_time * ns, sig_1, label='Regular')\n",
|
|
|
|
|
" ax.plot((double_signal_time) * ns, sig_2, label=\"Shifted\")\n",
|
|
|
|
|
" ax.plot((double_signal_time) * ns, signal_2(double_signal_time + calc_shift), 'r--', label='Recovered')\n",
|
|
|
|
|
" ax.set_xlabel(\"Time (ns)\")\n",
|
|
|
|
|
" ax.set_ylabel(\"Correlation\")\n",
|
|
|
|
|
" ax.set_title(\"Calculated Lag: {:.4g} ns\".format(calc_shift*ns))\n",
|
|
|
|
|
" ax.legend()\n",
|
|
|
|
|
" plt.show()"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"cell_type": "markdown",
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"source": [
|
|
|
|
|
"## new code"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"cell_type": "code",
|
|
|
|
|
"execution_count": 4,
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"outputs": [],
|
|
|
|
|
"source": [
|
|
|
|
|
"def lag_gridsearch(grid, sample_rate, signal_data, reference):\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\n",
|
|
|
|
|
" \n",
|
|
|
|
|
" corrs = grid_correlate(grid, signal_data, reference)\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, (err_min, err_plus)\n",
|
|
|
|
|
" \n",
|
|
|
|
|
"\n",
|
|
|
|
|
"def grid_correlate(grid, x, reference):\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\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",
|
2022-07-11 11:20:03 +02:00
|
|
|
" lags = signal.correlation_lags(in1_len, in2_len, mode=mode)\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
"\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",
|
|
|
|
|
" return np.linspace(start, end, grid_size, dtype=int, endpoint=False)"
|
|
|
|
|
]
|
|
|
|
|
},
|
2022-07-11 15:08:03 +02:00
|
|
|
{
|
|
|
|
|
"cell_type": "markdown",
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"source": [
|
|
|
|
|
"## Debug and test functions"
|
|
|
|
|
]
|
|
|
|
|
},
|
2022-06-14 17:03:00 +02:00
|
|
|
{
|
|
|
|
|
"cell_type": "code",
|
|
|
|
|
"execution_count": 5,
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"outputs": [
|
|
|
|
|
{
|
|
|
|
|
"data": {
|
2022-07-11 15:08:03 +02:00
|
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEWCAYAAACOv5f1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABylklEQVR4nO29d5xcV3n//37mTu/bd9WLLRdZtmzLxrTYBmzTgiHVhNAS4q8JhDT4UsPXISGFQEiHH04IJQRMIKHbgBOag40tg9yLZFtlpe270/vM+f1x7+xOuVNk7Wp2tef9euml2XPvmXl2duZzz33OU0QphUaj0WjWF45eG6DRaDSa048Wf41Go1mHaPHXaDSadYgWf41Go1mHaPHXaDSadYgWf41Go1mHaPHXaE4SEXmDiNx5CvNvE5HXL6dNGs3JosVfs2YRkV8Tkf0ikhKRCUtUn9dru2oRkZtF5N9qx5RSL1FKfbpXNmk0oMVfs0YRkT8A/gb4M2AE2AL8E3D9ST6Ps5sxjeZMQ4u/Zs0hIhHgA8BblFL/qZRKK6WKSqmvK6XeISIeEfkbETlh/fsbEfFYc68SkXEReaeITAL/aq3OvyQi/yYiCeANIhIRkX+x7iiOi8ifiojRwp6/FZFjIpIQkftE5PnW+IuB9wC/at2d3G+Nf19E3mQ9dojI+0TkiIhMi8hnrN8PEdkmIkpEXi8iR0VkVkTeW/O6l1t3PgkRmRKRv165d11zpqHFX7MWeTbgBf6rxfH3AlcAe4GLgMuB99UcHwX6ga3AjdbY9cCXgCjwOeDTQAk4C7gYuBZ4U4vXu9d6rX7g34H/EBGvUup2zDuTW5VSQaXURTZz32D9uxrYAQSBf2g453nAOcALgfeLyHnW+N8Cf6uUCgM7gS+2sE+jaUKLv2YtMgDMKqVKLY6/BviAUmpaKTUD/DHw2prjFeD/KaXySqmsNXaXUuorSqkKEAZeAvyedVcxDXwUuMHuxZRS/6aUmlNKlZRSHwE8mGLdDa8B/lop9ZRSKgW8G7ihwfX0x0qprFLqfuB+zAsaQBE4S0QGlVIppdTdXb6mRqPFX7MmmQMG2/jmNwBHan4+Yo1VmVFK5RrmHKt5vBVwARMiEhORGPD/AcN2LyYifygij4pI3Do3Agx2+bvY2erE3MeoMlnzOIN5dwDwm8Au4DERuVdEXt7la2o0Wvw1a5K7gBzwyhbHT2AKeJUt1lgVu1K2tWPHgDwwqJSKWv/CSqndjZMs//47gV8B+pRSUSAOSJvX6mRrCZjqMA+l1EGl1KsxL0p/CXxJRAKd5mk0oMVfswZRSsWB9wP/KCKvFBG/iLhE5CUi8iHg88D7RGRIRAatc/+t3XM2PP8E8B3gIyIStjZld4rIlTanhzDFegZwisj7Md1GVaaAbSLS6rv2eeD3RWS7iARZ2iNo5dJaRER+XUSGLFdVzBoud/M7ajRa/DVrEqXUXwN/gLmRO4O5Wn8r8BXgT4H9wAPAg8BPrbGT4XWAG3gEWMDcDB6zOe/bwG3AE5gumxz1LqT/sP6fE5Gf2sz/JPBZ4IfA09b83+nSxhcDD4tICnPz9wYbd5ZGY4voZi4ajUaz/tArf41Go1mHaPHXaDSadYgWf41Go1mHaPHXaDSadciaKWA1ODiotm3b1mszNBqNZk1x3333zSqlhhrH14z4b9u2jf379/faDI1Go1lTiMgRu3Ht9tFoNJp1iBZ/jUajWYdo8ddoNJp1yJrx+Ws0mpWlWCwyPj5OLqcrRKxFvF4vmzZtwuVydXW+Fn+NRgPA+Pg4oVCIbdu2ISKdJ2hWDUop5ubmGB8fZ/v27V3N0W4fjUYDQC6XY2BgQAv/GkREGBgYOKm7Ni3+Go1mES38a5eT/dtp8desLwppuOcWyMV7bYlG01O0+GvWFT/78ofhW2+n8uN/7LUpmlPg4x//OJ/5zGeaxg8fPswFF1zQA4tM3vCGN/ClL32p7Tmf+tSnOHFiqbHcm970Jh555JGVNq0JveGrWVf4j98JQO7g9/G/4D09tkbzTCiVStx0002n9fWcTmfLn0+WT33qU1xwwQVs2GC2lf7nf/7nU7bxmdD1yl9EPiki0yLyUM1Yv4h8V0QOWv/31Rx7t4gcEpHHReS6mvFLReRB69jfiXYyak4j/dmjADjmnuixJRo7/uRP/oRzzz2Xa665hle/+tV8+MMfBuCqq67iPe95D1deeSV/+7d/y80337x47L777uOiiy7i2c9+Nv/4j63v6D70oQ+xZ88eLrroIt71rncBcODAAa644gouvPBCXvWqV7GwsGD7eo0/33fffVx55ZVceumlXHfddUxMTDS93gc+8AEuu+wyLrjgAm688UaUUnzpS19i//79vOY1r2Hv3r1ks1muuuqqxdI1n//859mzZw8XXHAB73znOxefKxgM8t73vpeLLrqIK664gqmpji2eO3Iyl69PAf8A1N5rvQv4b6XUX4jIu6yf3yki5wM3ALuBDcAdIrJLKVUGPgbcCNwNfAuzFd1tp/qLaDQdKWYZKE+RVh4ChQXIxsAX7bVVq5I//vrDPHIisazPef6GMP/v53e3PL5//36+/OUv87Of/YxSqcQll1zCpZdeung8Fovxgx/8AICbb755cfyNb3wjf//3f8+VV17JO97xDtvnvu222/jKV77CT37yE/x+P/Pz8wC87nWvW5z7/ve/nz/+4z/mb/7mb5pe7+tf//riz8VikSuvvJKvfvWrDA0Nceutt/Le976XT37yk3Wv+da3vpX3v//9ALz2ta/lG9/4Br/0S7/EP/zDP/DhD3+Yffv21Z1/4sQJ3vnOd3LffffR19fHtddey1e+8hVe+cpXkk6nueKKK/jgBz/I//2//5dbbrmF973vfV28663peuWvlPohMN8wfD3waevxp4FX1ox/QSmVV0o9DRwCLheRMSCslLpLmf0jP1MzR6NZWeLjOFD8pHKe+XNysrf2aOq48847uf766/H5fIRCIX7+53++7viv/uqvNs2Jx+PEYjGuvPJKwBRZO+644w7e+MY34vf7Aejv72+a+/rXv54f/vCHLV+v+vPjjz/OQw89xDXXXMPevXv50z/9U8bHx5te83vf+x7Petaz2LNnD//zP//Dww8/3Pb3v/fee7nqqqsYGhrC6XTymte8ZtEet9vNy1/+cgAuvfRSDh8+3Pa5uuFUff4jSqkJAKXUhIgMW+MbMVf2VcatsaL1uHHcFhG5EfMugS1btpyiqZr1TikxiRN4TLbzAg6QWxjHO3xur81albRboa8UnfqJBwIB2zndeI67Pa/d61V/Vkqxe/du7rrrrpZzc7kcv/3bv83+/fvZvHkzN998c8cY/Ha/v8vlWrTfMAxKpVK3v0ZLVirax+5dVm3GbVFKfUIptU8ptW9oqKkctUZzUiTnLL/s6EUAHD/6VA+t0TTyvOc9j69//evkcjlSqRTf/OY3O86JRqNEIhHuvNPcyP/c5z5ne961117LJz/5STKZDADz8/NEIhH6+vr40Y9+BMBnP/vZxbuAdpxzzjnMzMwsin+xWGxa1VeFfnBwkFQqVRcBFAqFSCaTTc/7rGc9ix/84AfMzs5SLpf5/Oc/35U9z5RTXflPiciYteofA6at8XFgc815m4AT1vgmm3GNZsVJz0/QB2zd/SyYhrmJI+zstVGaRS677DJe8YpXcNFFF7F161b27dtHJBLpOO9f//Vf+Y3f+A38fj/XXXed7TkvfvGLOXDgAPv27cPtdvPSl76UP/uzP+PTn/40N910E5lMhh07dvCv//qvHV/P7XbzpS99ibe97W3E43FKpRK/93u/x+7dS3dL0WiU3/qt32LPnj1s27aNyy67bPHYG97wBm666SZ8Pl/d3cPY2Bh//ud/ztVXX41Sipe+9KVcf/31He15xiiluv4HbAMeqvn5r4B3WY/fBXzIerwbuB/wANuBpwDDOnYvcAXmXcBtwEu7ee1LL71UaTSnwsFb363K74+oB4/OqtT7h9TdH/s/vTZpVfHII4/02gSVTCaVUkql02l16aWXqvvuu6/HFq0t7P6GwH5lo6l
|
2022-06-14 17:03:00 +02:00
|
|
|
"text/plain": [
|
|
|
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
"metadata": {
|
|
|
|
|
"needs_background": "light"
|
|
|
|
|
},
|
|
|
|
|
"output_type": "display_data"
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"name": "stdout",
|
|
|
|
|
"output_type": "stream",
|
|
|
|
|
"text": [
|
2022-07-11 15:08:03 +02:00
|
|
|
"ids 87 871\n",
|
|
|
|
|
"ts -2.62e-08 -2.58e-08 =: -2.58e-08\n",
|
|
|
|
|
"avail ts -2e-07 - 1.978e-07 : 1.989e-09\n",
|
|
|
|
|
"grid size 200 2000 2000\n"
|
2022-06-14 17:03:00 +02:00
|
|
|
]
|
|
|
|
|
}
|
|
|
|
|
],
|
|
|
|
|
"source": [
|
|
|
|
|
"#################\n",
|
2022-07-11 15:08:03 +02:00
|
|
|
"grid_size = 200\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
"my_mode = 'same'\n",
|
|
|
|
|
"grid = correlation_grid(grid_size, in1_len=len(sig_1), in2_len=len(sig_2), mode=my_mode)\n",
|
|
|
|
|
"t_grid = grid/signal_1.sample_rate\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"scipy_mode = my_mode\n",
|
2022-07-11 11:20:03 +02:00
|
|
|
"scipy_corr = signal.correlate(sig_2, sig_1, mode=scipy_mode)\n",
|
|
|
|
|
"scipy_lags = signal.correlation_lags(len(sig_2), len(sig_1), mode=scipy_mode)\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
"corrs = grid_correlate(grid, sig_2, sig_1)\n",
|
|
|
|
|
" \n",
|
|
|
|
|
"\n",
|
|
|
|
|
"scipy_id = np.argmax(scipy_corr)\n",
|
|
|
|
|
"my_id = np.argmax(corrs)\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"fig, ax = plt.subplots()\n",
|
|
|
|
|
"ax.set_title(\"Correlations\")\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"if not True:\n",
|
|
|
|
|
" x_grid = grid\n",
|
|
|
|
|
" scipy_grid = scipy_lags\n",
|
|
|
|
|
" ax.set_xlabel(\"Index Lag\")\n",
|
|
|
|
|
"else:\n",
|
|
|
|
|
" x_grid = t_grid * 1e9\n",
|
|
|
|
|
" scipy_grid = scipy_lags / signal_1.sample_rate * 1e9\n",
|
|
|
|
|
" ax.set_xlabel(\"Time [ns]\")\n",
|
|
|
|
|
" \n",
|
|
|
|
|
"ax.plot(x_grid, corrs, label='grid correlation')\n",
|
|
|
|
|
"ax.plot(scipy_grid, scipy_corr, label='scipy')\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"ax.legend()\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"plt.show();\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"if scipy_mode == 'same':\n",
|
|
|
|
|
" scipy_corr_grid = scipy_corr[grid]\n",
|
|
|
|
|
"else:\n",
|
|
|
|
|
" scipy_corr_grid = scipy_corr[grid]\n",
|
|
|
|
|
" \n",
|
|
|
|
|
"\n",
|
|
|
|
|
"if len(grid) <= 10:\n",
|
|
|
|
|
" print(\"G\", grid)\n",
|
|
|
|
|
" print(\"TG\", grid/signal_1.sample_rate)\n",
|
|
|
|
|
" print(\"C\", corrs)\n",
|
|
|
|
|
" print(\"SC\", scipy_corr_grid)\n",
|
|
|
|
|
" \n",
|
|
|
|
|
" \n",
|
|
|
|
|
"print(\"ids\", my_id, scipy_id)\n",
|
|
|
|
|
"print(\"ts\", t_grid[my_id], scipy_lags[np.argmax(scipy_corr)]/signal_1.sample_rate, '=:', signal_2.t_0)\n",
|
|
|
|
|
"print(\"avail ts\", t_grid[0], '-', t_grid[-1] ,':', (t_grid[-1] - t_grid[0])/len(t_grid))\n",
|
|
|
|
|
"print(\"grid size\", len(t_grid), len(sig_2), len(sig_1))"
|
|
|
|
|
]
|
|
|
|
|
},
|
2022-07-11 15:08:03 +02:00
|
|
|
{
|
|
|
|
|
"cell_type": "markdown",
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"source": [
|
|
|
|
|
"## Showcase functions"
|
|
|
|
|
]
|
|
|
|
|
},
|
2022-06-14 17:03:00 +02:00
|
|
|
{
|
|
|
|
|
"cell_type": "code",
|
|
|
|
|
"execution_count": 6,
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"outputs": [
|
|
|
|
|
{
|
|
|
|
|
"name": "stdout",
|
|
|
|
|
"output_type": "stream",
|
|
|
|
|
"text": [
|
2022-07-11 15:08:03 +02:00
|
|
|
"Known time offset: -2.58e-08\n",
|
|
|
|
|
"Calculated offset: -2.62e-08\n",
|
|
|
|
|
"Uncertainty on offset: (-1e-09, 1e-09)\n"
|
2022-06-14 17:03:00 +02:00
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"data": {
|
2022-07-11 15:08:03 +02:00
|
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7gAAAEGCAYAAABPZjl3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAC9IklEQVR4nOz9eZgcZ3nvjX+e3vfu2WekGWm0epFtyZa8CGNsdghLArFjQs4hzoYJyYmTnGwc3jcQfi9nSQghBBMMYUlOSGzsxGEJAQcTwmZsS95lW9YuzWj26X1f6vfH0zVrr1XVM9L4+VzXXKPprq56eqa7Vd+67/v7FZqmoVAoFAqFQqFQKBQKxcWObb0XoFAoFAqFQqFQKBQKhRUogatQKBQKhUKhUCgUig2BErgKhUKhUCgUCoVCodgQKIGrUCgUCoVCoVAoFIoNgRK4CoVCoVAoFAqFQqHYEDjWewGdoLe3VxsdHV3vZSgUHefkTHrZz9v7/Ou0EoVCoVAoFAqFYm04fPjwrKZpfbXu25ACd3R0lEOHDq33MhSKjnP7PY8s+/m+Ow+u00oUCoVCoVAoFIq1QQhxpt59qkVZobiIue/Og0rUKhQKhUKhUCgUVZTAVSg2EPcfOsf9h86t9zIUCoVCoVAoFIp1QQlchWID8cDhMR44PLbey1AoFAqFQqFQKNaFDTmDq1AoFAqFQqFQKBRrQbFYZGxsjFwut95L2XB4PB6Gh4dxOp0tP0YJXIVCoVAoFAqFQqEwyNjYGMFgkNHRUYQQ672cDYOmaczNzTE2Nsa2bdtafty6tSgLIUaEEP8hhHhBCHFECHFXjW2EEOKTQojjQohnhBDXrMdaFQqFQqFQKBQKhaIWuVyOnp4eJW4tRghBT09P25Xx9azgloD/rmnaE0KIIHBYCPHvmqY9v2SbNwO7ql/XA39d/a5QKFgdE6RQKBQKhUKhWHuUuO0MRn6v6yZwNU2bACaq/04KIV4ANgNLBe5PA3+naZoG/EQIERFCDFUfq1AoVnDr/uH1XoJ5MvMwcxS6RiE0tN6rMU25ovHCRAIh4LLBEDab+g9QoVAoFAqFolNcEC7KQohR4Grg0RV3bQaWZp6MVW+rtY/3CiEOCSEOzczMdGSdCsWFzm0HRrjtwMh6L8M42Sg88bdw8nvw5P+F2Nn1XpFpvn1kkn9/foqHjkzx7y9MrfdyFAqFQqFQbEDsdjv79u3jiiuu4G1vexuxWMzyY9xyyy0cOnTI8v1azboLXCFEAPgn4Lc1TUusvLvGQ7Ra+9E07bOaph3QNO1AX1+f1ctUKBRrwcn/hEoJ9v8iuINw/Dug1XzLXxScm8/w1RP/xKT9Qa7ZGub58wkm48phUaFQKBQKhbV4vV6eeuopnnvuObq7u7n77rvXe0mUy+V1Oe66ClwhhBMpbr+sado/19hkDFhajhoGzq/F2hQKxRqTmYfpF2D4Wghtgq2vgOQUzJ9c75UZ5pGTMzyV+gozxZe4YXsvHqedx07Pr/eyFAqFQqFQbGAOHjzI+Pg4ACdOnOBNb3oT+/fv56abbuLFF19cuP2GG27g2muv5Y//+I8JBAIAfO973+Otb33rwr5+8zd/ky996UurjvHrv/7rHDhwgD179vChD31o4fbR0VE+8pGP8MpXvpL777+/g8+yPus2gyvkxPDngRc0Tft4nc2+BvymEOJepLlUXM3fKhQblOnq+P2mqln6wBVw4j9g6jno2bF+6zJIKl/isclHyZSj/MKl78btsHPF5hBPnImRK5bxOO3rvUSFQqFQKBQW872j08wk85busy/o5pZL+lvatlwu8/DDD/Mrv/IrALz3ve/lM5/5DLt27eLRRx/l/e9/P9/97ne56667uOuuu/j5n/95PvOZz7S9po9+9KN0d3dTLpd57WtfyzPPPMNVV10FyOzaH/7wh23v0yrWs4J7I/BfgdcIIZ6qfv2UEOJ9Qoj3Vbf5JnASOA58Dnj/Oq1Vobggue/Og9x358H1XoY1TL8AkRHwhOTPNjv07oa541Aure/aDHBsKsnJzI/wO/y8cviVAOweCFLRNI5Pp9Z5dQqFQqFQKDYS2WyWffv20dPTw/z8PK9//etJpVL8+Mc/5rbbbmPfvn3ceeedTEzIWuEjjzzCbbfdBsC73/3uto/3la98hWuuuYarr76aI0eO8Pzziz7Bt99+uzVPyiDr6aL8Q2rP2C7dRgN+Y21WpFBc/Nx/SHqyXXRGU7k4pGdh52uX3953CUw8DbEzF10V99RsivOFJ7lp5CbcdjcA/UE3Ia+Tk7NprtgcXucVKhQKhUKhsJpWK61Wo8/gxuNx3vrWt3L33Xdzxx13EIlEeOqpp1rej8PhoFKpLPxcK4P21KlTfOxjH+Pxxx+nq6uLO+64Y9l2fr/f1HMxy7qbTCkUCut44PAYDxweW+9ltE/0tPzeNbr89sgWWcmNnVnrFZmiXNE4MnOMdDnKDUM3LNwuhGBLt4+xaIZK5eI1z1IoFAqFQnFhEg6H+eQnP8nHPvYxvF4v27ZtW5iF1TSNp59+GoAbbriBf/qnfwLg3nvvXXj81q1bef7558nn88TjcR5++OFVx0gkEvj9fsLhMFNTU/zbv/3bGjyz1lECV6FQrD/RM+DygX+FA7rdKQ2noheXwJ2IZzmbeQaAGzbdsOy+kW4v+WKF2ZS18zmKlyGlPLz4r/Dkly/MSK3xw3D4b+HMIy27oR85H+fex87y4xOzaBeYg/pcKs+fff9rvPWBd/HHP/wTsqXsei9JoVAoanL11Vezd+9e7r33Xr785S/z+c9/nr1797Jnzx6++tWvAvCJT3yCj3/841x33XVMTEwQDsvOspGREX7u536Oq666il/4hV/g6quvXrX/vXv3cvXVV7Nnzx5++Zd/mRtvvHFNn18z1q1FWaFQKBZIjEN4BESNqYXIFjjzY3ky73Cv/doMMBHPMVV8kQHfIJsDy6O7h7t8AJyLZukPedZjeYqNwrGHYOp5eXHoma/Adb8Gnguk9X32GLz0kJypP/k9cHpg0+qTpKWMRTM8dGSKoMfBoyfncdltHBjtXpv1NqFQqvC1p8+TyboplDQePPFPIMp85MaPrPfSFAqFAoBUarm/x9e//vWFf3/rW99atf3mzZv5yU9+ghCCe++9lwMHDizc96d/+qf86Z/+6arHfO9731v4dy1nZYDTp0+3t/AOoCq4CoVifSlkIBuTldpahDbL6k9yck2XZYapRI754gmu7L1i1X0Bt4Ogx8F0QuXhKkwQH4fJ52DLDXD1fwGtIi8EXQhUKnDs3yHQD9fdCeFhOP2jhmZxmqbxvaMzhLxO3nNwlO19fh47PU+uuD4Ziit58myUWKbIL197E5941efY438rXz3+VU7HT6/30hQKhcIQhw8fZt++fVx11VV8+tOf5s///M/Xe0mWoQSuQnERc/s9j3D7PY+s9zLMkawmfwUHa9+v356aWpv1WMCZ6Bzx0iR7evfUvH8g5GFKCVyFGc4/KVv4txwEbxcM7YXJZ6F4AbTNzh2XxnGjrwS7Q37PJ2HmhboPmUzkmEnmuW60G5fDxsEdPeSLFY6cj6/hwmtTqWg8Ox5na4+PkW4flw4GuanvVmzCwd+/8PfrvTyFQqEwxE033cTTTz/NM888w/e//3127ty53kuyDCVwFYoNxK37h7l1//B6L6M99MpsoI7Adfllm2Py4ojAzhRKnErKEPU9PbUFbn/QTTRTvGCqU4qLjHJRisWBPeBwyduG9kKlDDMvru/aACafAXcAenbJn7tG5Xt4uv7anj+f4Fj2IQIBKWj7gx4GQh6OTq5/pNaZ+QzJXImrhiMA2GyCa7eMMuzez7+f+Q7linofKxQKxYWEErgKxQbitgMjF19EUHICfN1yRq8ewcGLpkV5KpFntnACgMt7Lq+5zUDIw+nsI7z5n99APL/+FSrFRUb0jGz37btk8bbAAPh6GorINaFchOgp6L0EbNVTDCHkWqOnoLi6c0HTNJ6bHOcH0Xt46MyiE+clg0GmEjlimcJarb4mJ2dSuBw2Rnt8C7ft6g+wzfsK5nNzPDH9xDq
|
2022-06-14 17:03:00 +02:00
|
|
|
"text/plain": [
|
|
|
|
|
"<Figure size 1152x288 with 1 Axes>"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
"metadata": {
|
|
|
|
|
"needs_background": "light"
|
|
|
|
|
},
|
|
|
|
|
"output_type": "display_data"
|
|
|
|
|
}
|
|
|
|
|
],
|
|
|
|
|
"source": [
|
2022-07-11 15:08:03 +02:00
|
|
|
"grid_size = 200\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"grid = correlation_grid(grid_size, in1_len=len(sig_1), in2_len=len(sig_2), mode='same')\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
"calc_shift, err = lag_gridsearch(grid, signal_1.sample_rate, sig_2, sig_1)\n",
|
2022-07-11 15:08:03 +02:00
|
|
|
"\n",
|
|
|
|
|
"print(\"Known time offset:\", signal_2.t_0)\n",
|
|
|
|
|
"print(\"Calculated offset:\", calc_shift)\n",
|
|
|
|
|
"print(\"Uncertainty on offset:\", err)\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
"\n",
|
|
|
|
|
"if True:\n",
|
|
|
|
|
" signal_2.periodic = False\n",
|
|
|
|
|
"\n",
|
2022-07-11 15:08:03 +02:00
|
|
|
"fig, ax = plt.subplots(figsize=(16,4))\n",
|
|
|
|
|
"ax.plot(double_signal_time * ns, sig_1, alpha=0.5, label='Regular')\n",
|
|
|
|
|
"ax.plot(double_signal_time * ns, sig_2, alpha=0.5, label='Shifted')\n",
|
|
|
|
|
"ax.plot((double_signal_time) * ns, signal_2(double_signal_time + calc_shift), ls=(0, (10,)), label='Recovered')\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
"\n",
|
|
|
|
|
"\n",
|
2022-07-11 15:08:03 +02:00
|
|
|
"ax.axvline(calc_shift * ns)\n",
|
|
|
|
|
"for i, _ in enumerate(err):\n",
|
|
|
|
|
" ax.axvline((calc_shift + err[i]) * ns, ls=(5*i, (4, 4)))\n",
|
2022-06-14 17:03:00 +02:00
|
|
|
"\n",
|
2022-07-11 15:08:03 +02:00
|
|
|
"ax.set_xlabel(\"Time (ns)\")\n",
|
|
|
|
|
"ax.set_ylabel(\"Amplitude\")\n",
|
|
|
|
|
"ax.legend()\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"plt.show()"
|
2022-06-14 17:03:00 +02:00
|
|
|
]
|
|
|
|
|
}
|
|
|
|
|
],
|
|
|
|
|
"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.10.2"
|
|
|
|
|
}
|
|
|
|
|
},
|
|
|
|
|
"nbformat": 4,
|
|
|
|
|
"nbformat_minor": 4
|
|
|
|
|
}
|
