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)"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"cell_type": "code",
|
|
|
|
|
"execution_count": 5,
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"outputs": [
|
|
|
|
|
{
|
|
|
|
|
"data": {
|
|
|
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEWCAYAAACOv5f1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABw30lEQVR4nO29d5xcZ33v//6e6X22qndLcpdsycaUxKabEgw3BbhgauJLaMnNDQEC4TokpHBDgDT4ATEloSUQMBBMMcG02NiSu2XLlq22Ktun95nn98c5Z3fKmSLtrnZX+7xfL700+5zzzHx3duZznvN9vkWUUmg0Go1mZWEstgEajUajOfdo8ddoNJoViBZ/jUajWYFo8ddoNJoViBZ/jUajWYFo8ddoNJoViBZ/jeYMEZE3iMjP5zD/NhF5/XzapNGcKVr8NcsWEfmfIrJPRDIicsoS1Wcttl31iMjNIvKv9WNKqRcppT6/WDZpNKDFX7NMEZE/AD4G/AWwCtgI/BNwwxk+j7uXMY3mfEOLv2bZISIx4IPA25RS/6GUyiqlykqpbyul3iUiPhH5mIictP59TER81tzrRGRERN4tIqeBz1qr86+JyL+KSAp4g4jEROSfrTuKEyLy5yLiamPPx0XkuIikRGS/iPyKNX498MfAK627kwes8TtE5Letx4aIvF9EjorImIh8wfr9EJHNIqJE5PUickxEJkTkfXWve7V155MSkVER+duFe9c15xta/DXLkacDfuAbbY6/D7gG2A3sAq4G3l93fDXQD2wCbrLGbgC+BsSBLwKfByrABcAVwAuA327zevdYr9UPfAn4dxHxK6W+h3ln8lWlVFgptcth7husf88GtgJh4B+aznkWsBN4LvABEbnIGv848HGlVBTYBvxbG/s0mha0+GuWIwPAhFKq0ub4a4APKqXGlFLjwJ8CN9YdrwH/VylVVErlrbE7lVLfVErVgCjwIuD3rbuKMeCjwKucXkwp9a9KqUmlVEUp9RHAhynWvfAa4G+VUk8ppTLAe4FXNbme/lQplVdKPQA8gHlBAygDF4jIoFIqo5S6q8fX1Gi0+GuWJZPAYAff/FrgaN3PR60xm3GlVKFpzvG6x5sAD3BKRBIikgD+P2DY6cVE5P+IyKMikrTOjQGDPf4uTra6MfcxbE7XPc5h3h0AvBnYATwmIveIyEt7fE2NRou/ZllyJ1AAXt7m+ElMAbfZaI3ZOJWyrR87DhSBQaVU3PoXVUpd0jzJ8u+/G/gtoE8pFQeSgHR4rW62VoDRLvNQSj2hlHo15kXpr4GviUio2zyNBrT4a5YhSqkk8AHgH0Xk5SISFBGPiLxIRD4MfBl4v4gMicigde6/dnrOpuc/BfwA+IiIRK1N2W0icq3D6RFMsR4H3CLyAUy3kc0osFlE2n3Xvgz8bxHZIiJhZvcI2rm0ZhCR14rIkOWqSljD1V5+R41Gi79mWaKU+lvgDzA3cscxV+tvB74J/DmwD3gQeAi41xo7E14HeIEDwDTmZvAah/O+D9wGPI7psinQ6EL6d+v/SRG512H+LcC/AD8FDlvz39GjjdcDj4hIBnPz91UO7iyNxhHRzVw0Go1m5aFX/hqNRrMC0eKv0Wg0KxAt/hqNRrMC0eKv0Wg0K5BlU8BqcHBQbd68ebHN0Gg0mmXF/v37J5RSQ83jy0b8N2/ezL59+xbbDI1Go1lWiMhRp3Ht9tFoNJoViBZ/jUajWYFo8ddoNJoVyLLx+Ws0moWlXC4zMjJCoaArRCxH/H4/69evx+Px9HS+Fn+NRgPAyMgIkUiEzZs3IyLdJ2iWDEopJicnGRkZYcuWLT3N0W4fjUYDQKFQYGBgQAv/MkREGBgYOKO7Ni3+Go1mBi38y5cz/dtp8desKHKZJL/86l+RSkwutikazaKixV+zonjwGx/haY/+JY98/S8X2xTNHPjkJz/JF77whZbxI0eOcOmlly6CRSZveMMb+NrXvtbxnM997nOcPDnbWO63f/u3OXDgwEKb1oLe8NWsKIInfgFAbFT3Ol+uVCoV3vKWt5zT13O73W1/PlM+97nPcemll7J2rdlW+jOf+cycbTwbel75i8gtIjImIg/XjfWLyA9F5Anr/766Y+8VkUMiclBEXlg3vkdEHrKO/Z1oJ6PmHDJYPAbAmpJjxrtmkfmzP/szLrzwQp7//Ofz6le/mr/5m78B4LrrruOP//iPufbaa/n4xz/OzTffPHNs//797Nq1i6c//en84z/+Y9vn/vCHP8xll13Grl27eM973gPA/fffzzXXXMPll1/OK17xCqanpx1fr/nn/fv3c+2117Jnzx5e+MIXcurUqZbX++AHP8hVV13FpZdeyk033YRSiq997Wvs27eP17zmNezevZt8Ps911103U7rmy1/+MpdddhmXXnop7373u2eeKxwO8773vY9du3ZxzTXXMDratcVzV87k8vU54B+A+nut9wA/Ukr9lYi8x/r53SJyMfAq4BJgLXC7iOxQSlWBTwA3AXcB38VsRXfbXH8RjaYbhVyG1bVxcvjokxTJ6QlifYOLbdaS5E+//QgHTqbm9TkvXhvl//7aJW2P79u3j69//evcd999VCoVrrzySvbs2TNzPJFI8JOf/ASAm2++eWb8jW98I3//93/Ptddey7ve9S7H577tttv45je/yS9/+UuCwSBTU1MAvO51r5uZ+4EPfIA//dM/5WMf+1jL633729+e+blcLnPttddy6623MjQ0xFe/+lXe9773ccsttzS85tvf/nY+8IEPAHDjjTfyne98h9/4jd/gH/7hH/ibv/kb9u7d23D+yZMnefe7383+/fvp6+vjBS94Ad/85jd5+ctfTjab5ZprruFDH/oQf/RHf8SnP/1p3v/+9/fwrren55W/UuqnwFTT8A3A563HnwdeXjf+FaVUUSl1GDgEXC0ia4CoUupOZfaP/ELdHI1mQRkbOYQhiieCuwBInNar/6XEz3/+c2644QYCgQCRSIRf+7Vfazj+yle+smVOMpkkkUhw7bXXAqbIOnH77bfzxje+kWAwCEB/f3/L3Ne//vX89Kc/bft69s8HDx7k4Ycf5vnPfz67d+/mz//8zxkZGWl5zR//+Mc87WlP47LLLuO//uu/eOSRRzr+/vfccw/XXXcdQ0NDuN1uXvOa18zY4/V6eelLXwrAnj17OHLkSMfn6oW5+vxXKaVOASilTonIsDW+DnNlbzNijZWtx83jjojITZh3CWzcuHGOpmpWOplJ89Y8N3ApjNxNanwELtrTZdbKpNMKfaHo1k88FAo5zunFc9zreZ1ez/5ZKcUll1zCnXfe2XZuoVDgrW99K/v27WPDhg3cfPPNXWPwO/3+Ho9nxn6Xy0WlUun112jLQkX7OL3LqsO4I0qpTyml9iql9g4NtZSj1mjOiELiNAD+DbvNn6eOL6I1mmae9axn8e1vf5tCoUAmk+E///M/u86Jx+PEYjF+/vOfA/DFL37R8bwXvOAF3HLLLeRyOQCmpqaIxWL09fXxs5/9DIB/+Zd/mbkL6MTOnTsZHx+fEf9yudyyqreFfnBwkEwm0xABFIlESKfTLc/7tKc9jZ/85CdMTExQrVb58pe/3JM9Z8tcV/6jIrLGWvWvAcas8RFgQ91564GT1vh6h3GNZsEpp8xNsuHtV8OdUEnqj95S4qqrruJlL3sZu3btYtOmTezdu5dYLNZ13mc/+1ne9KY3EQwGeeELX+h4zvXXX8/999/P3r178Xq9vPjFL+Yv/uIv+PznP89b3vIWcrkcW7du5bOf/WzX1/N6vXzta1/jne98J8lkkkqlwu///u9zySWzd0vxeJzf+Z3f4bLLLmPz5s1cddVVM8fe8IY38Ja3vIVAINBw97BmzRr+8i//kmc/+9kopXjxi1/MDTfc0NWes0Yp1fM/YDPwcN3P/w94j/X4PcCHrceXAA8APmAL8BTgso7dA1yDeRdwG/DiXl57z549SqOZC//9mT9Q1Q/EVLlUVNkPDKk7/+l/LbZJS4oDBw4stgkqnU4rpZTKZrNqz549av/+/Yts0fLC6W8I7FMOmtrzyl9
|
|
|
|
|
"text/plain": [
|
|
|
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
"metadata": {
|
|
|
|
|
"needs_background": "light"
|
|
|
|
|
},
|
|
|
|
|
"output_type": "display_data"
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"name": "stdout",
|
|
|
|
|
"output_type": "stream",
|
|
|
|
|
"text": [
|
|
|
|
|
"ids 872 871\n",
|
|
|
|
|
"ts -2.58e-08 -2.58e-08 =: -2.58e-08\n",
|
|
|
|
|
"avail ts -2e-07 - 1.996e-07 : 1.9979999999999997e-10\n",
|
|
|
|
|
"grid size 2000 2000 2000\n"
|
|
|
|
|
]
|
|
|
|
|
}
|
|
|
|
|
],
|
|
|
|
|
"source": [
|
|
|
|
|
"#################\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"grid_size = 2000\n",
|
|
|
|
|
"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))"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"cell_type": "code",
|
|
|
|
|
"execution_count": 6,
|
|
|
|
|
"metadata": {},
|
|
|
|
|
"outputs": [
|
|
|
|
|
{
|
|
|
|
|
"name": "stdout",
|
|
|
|
|
"output_type": "stream",
|
|
|
|
|
"text": [
|
|
|
|
|
"-2.58e-08 = -2.58e-08 \\pm (-1e-10, 1e-10)\n"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
{
|
|
|
|
|
"data": {
|
|
|
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7gAAAEGCAYAAABPZjl3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAC2s0lEQVR4nOz9d3hkV5XvD3935RwVW7Fzu7O726ExxgYbkwyMwWBgAswwYGC4mMt9mfBy78CP++PeeRlgGA9mMAwDA8OMSWOiAYPB2ODUDm23O6qDpFYuqXJO+/1j15FKUoWTSlKr1+d59Ki7ap9zdkl1Svu711rfxTjnIAiCIAiCIAiCIIhLHcNqT4AgCIIgCIIgCIIg9IAELkEQBEEQBEEQBLEuIIFLEARBEARBEARBrAtI4BIEQRAEQRAEQRDrAhK4BEEQBEEQBEEQxLrAtNoTaAVtbW18cHBwtadBECvG+VAKmUIJdrMRm9qdqz0dgiAIgiAIgmgZzzzzzCznvL3Wc+tS4A4ODuLpp59e7WkQxIpxx72P48RkHDu7Pfj2nYdXezoEQRAEQRAE0TIYYyP1nqMUZYJYB3z7zsPY2e1Z7WkQBEEQBEEQxKqyLiO4BHE5EkrkVnsKBEEQBEEQBLGqkMAliHVCKEkClyAIgiAIgri8IYFLEARBEARBEAShkkKhgLGxMWSz2dWeyrrDZrOht7cXZrNZ9jEkcAmCIAiCIAiCIFQyNjYGt9uNwcFBMMZWezrrBs455ubmMDY2ho0bN8o+btVMphhjfYyx3zDGTjLGjjPG7qoxhjHG7maMnWWMvcAYO7AacyUIgiAIgiAIgqhFNptFMBgkcaszjDEEg0HFkfHVdFEuAvgfnPMrAFwL4C8YYzuXjHkNgK2Vr/cC+OeVnSJBXBrcce/jSOeKqz0NgiAIgiCIyxISt61Bzc911QQu53ySc/5s5d8JACcB9CwZ9kYA3+CCJwD4GGPdKzxVgrgkMBsNaHdZV3sa2kmHgZHHgfjkas9EF0pljhfHYzg+EUO5zFd7OgRBEARBEOuaNdEHlzE2COBKAE8ueaoHwMWq/49huQiWzvFextjTjLGnQ6FQS+ZJEGsZs8mAdvclLnAzEeDZfwPOPww8900gOrraM9LM//z1PXjgxWE8eHwavzw5vdrTIQiCIAhiHWI0GrF//37s3r0br3/96xGNRnW/xo033oinn35a9/PqzaoLXMaYC8D3AXyYcx5f+nSNQ2qGQDjnX+acH+KcH2pvb9d7mgRBrATnfwuUi8DBdwJWN3D2VwC/dKOePz3zKH4yfi9Kzqdx1WAAJybimIqRwyJBEARBEPpit9tx9OhRvPjiiwgEArjnnntWe0oolUqrct1VFbiMMTOEuP0W5/y/agwZA9BX9f9eABMrMTeCIFaYdBiYOQn0XgV4NgADLwES00D4/GrPTDVffuErcBoD+OBVb8dVG/2wmY14aji82tMiCIIgCGIdc/jwYYyPjwMAzp07h1e/+tU4ePAgrr/+epw6dWr+8WuvvRZXXXUV/vZv/xYulwsA8PDDD+PWW2+dP9cHP/hBfP3rX192jfe///04dOgQdu3ahY9//OPzjw8ODuKTn/wkXvrSl+K73/1uC19lfVatTRATFcNfBXCSc/65OsN+BOCDjLH7AFwDIMY5Xx+FeQRBLGbmhPi+oWKW3rkbOPcbYPpFILh59ealkuHoJM6njuJ1fX8Mt9UBANjd48GzI1FkCyXYzMZVniFBEARBEHrz8OkZhBI5Xc/Z7rbixu0dssaWSiU89NBDePe73w0AeO9734svfelL2Lp1K5588kl84AMfwK9//WvcdddduOuuu/D2t78dX/rSlxTP6VOf+hQCgQBKpRJuuukmvPDCC9i7dy8A0bv2d7/7neJz6sVqRnCvA/DHAF7BGDta+XotY+x9jLH3VcY8AOA8gLMAvgLgA6s0V4JY03z7zsPY2e1Z7WloY+Yk4OsDbJXXYTACbduAubNA6dJziP72iZ8A4HjrjtfPP7at040y5zg7k1y9iREEQRAEse7IZDLYv38/gsEgwuEwXvnKVyKZTOKxxx7DW97yFuzfvx933nknJidFrPDxxx/HW97yFgDAO97xDsXX+853voMDBw7gyiuvxPHjx3HixIn55+644w59XpRKVi2Cyzn/HWrX2FaP4QD+YmVmRBCXNnrvFq4o2RiQmgW23LT48fbtwOTzQHTkkoviPjH1GALmHhzYsGP+sQ63FR67GednU9jd413F2REEQRAE0QrkRlr1RqrBjcViuPXWW3HPPffgXe96F3w+H44ePSr7PCaTCeVyef7/tXrQXrhwAZ/5zGdw5MgR+P1+vOtd71o0zul0anotWll1kymCIPQhlMwhlLxERW5kWHz3Dy5+3NcvIrnRkZWekSayhTyGk8ewK3Bo0eOMMfQHHBiLpKllEEEQBEEQuuP1enH33XfjM5/5DOx2OzZu3DhfC8s5x/PPPw8AuPbaa/H9738fAHDffffNHz8wMIATJ04gl8shFovhoYceWnaNeDwOp9MJr9eL6elp/OxnP1uBVyYfErgEQaw+kRHA4gCcSxzQjWZhOBW5tATuwyNHUOQ5XNdzeNlzfQE7coUyZi/VzQhi7VDMAad+Cjz3rbXZUmv8GeCZfxN9rWW6oR+fiOG+p0bx2LlZ8DXmoD6XzOHvfvt9vO0H78P58PhqT4cgCKIuV155Jfbt24f77rsP3/rWt/DVr34V+/btw65du/DDH/4QAPD5z38en/vc53D11VdjcnISXq/ILOvr68Nb3/pW7N27F3/4h3+IK6+8ctn59+3bhyuvvBK7du3Cn/3Zn+G6665b0dfXjFVLUSYIgpgnPg54+wBWo2rB1w+MPCYW86ZLo8/vkxPPAQBeMXjtsud6/cJw6mIkgw6PbUXnRawzhh4Epk+IzaEXvgNc/R7AtkZS32eHgDMPipr68w8DZhuwYfkiqZqxSBoPHp+G22bCk+fDsBgNODQYWJn5NiFfLONHz0+glAvidPwZ3PnL/4afv+W7MBrILI4giLVBMrnY3+PHP/7x/L9//vOfLxvf09ODJ554Aowx3HfffTh0aCHr7NOf/jQ+/elPLzvm4Ycfnv93LWdlABgeHlY28RZAEVyCIFaXfBrIREWkthaeHhH9SUyt6LS0cCpyHF5TF7rdwWXPuawmuG0mzMSpHy6hgdg4MPUi0H8tcOUfAbwsNoLWAuUyMPRLwNUBXH0n4O0Fhn/f0CyOc46HT4fgsZvxJ4cHsandiaeGw8gWVqeH4lKeG40gmi7gTw5djTt3fhRT2SF869hPV3taBEEQqnnmmWewf/9+7N27F1/84hfx2c9+drWnpBskcAliHXDHvY8jnbv0nIYBAIlK5y93V+3npceT0yszHx0YSZ7BgHt73ec7PTZMk8AltDDxnEjh7z8M2P1A9z5g6hhQyKz2zITzeTYGDL4UMJrE91wCCJ2se8gLUyMYjU7j6sEALCYDDm8OIlco4/hEbAUnXptymePYeAwDQQf6Ag68e/+b4DdvwH+e/vfVnhpBEIRqrr/+ejz//PN44YUX8Mgjj2DLli2rPSXdIIFLEOsEs9GAdtelkcK7CCky66ojcC1OkeaYuDRaYI/HZ5AozuAK/866YzrcVkTShTUTnSIuMUoFIRY7dwEmi3isex9QLgGhU6s7NwCYegGwuoDgVvF//6C4h2fqz+2Lz92L7898CIPt4vV0uG3o9Nhwemr1W2qNhNNIZIvY2+sDAJhNJry857UYy5zEUPji6k6OIAiCWAYJXIJYJ5hNBrS7L0WBOwk4AqJGrx7urksmRfnxMeFOeGXnnrpjOj02FMs5jITnVmpaxHoiMiLSfdursgRcnYAj2FBErgilAhC5ALRtBwyVJQZjYq6RC0BheeZCuVzGC+HHsc17JdxW+/zj27vcmI5nEU3nV2r2NTkfSsJiMmAw6Jh/7M3bXwcA+N7JB1ZrWgRBEEQdSOASBLG6JKYAd3fjMa4uIB2uuTheaxydeREAw7V9e+uOcdjy+ObUH+E7Z763chMj1g/hcyL119u/8Bhjold
|
|
|
|
|
"text/plain": [
|
|
|
|
|
"<Figure size 1152x288 with 1 Axes>"
|
|
|
|
|
]
|
|
|
|
|
},
|
|
|
|
|
"metadata": {
|
|
|
|
|
"needs_background": "light"
|
|
|
|
|
},
|
|
|
|
|
"output_type": "display_data"
|
|
|
|
|
}
|
|
|
|
|
],
|
|
|
|
|
"source": [
|
|
|
|
|
"##########################\n",
|
|
|
|
|
"##########################\n",
|
|
|
|
|
"calc_shift, err = lag_gridsearch(grid, signal_1.sample_rate, sig_2, sig_1)\n",
|
|
|
|
|
"print(signal_2.t_0, \"=\", calc_shift, \"\\pm\", err)\n",
|
|
|
|
|
"\n",
|
|
|
|
|
"if True:\n",
|
|
|
|
|
" signal_2.periodic = False\n",
|
|
|
|
|
"\n",
|
|
|
|
|
" 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",
|
|
|
|
|
"\n",
|
|
|
|
|
" \n",
|
|
|
|
|
" 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",
|
|
|
|
|
"\n",
|
|
|
|
|
" ax.set_xlabel(\"Time (ns)\")\n",
|
|
|
|
|
" ax.set_ylabel(\"Amplitude\")\n",
|
|
|
|
|
" ax.legend()\n",
|
|
|
|
|
"\n",
|
|
|
|
|
" plt.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.10.2"
|
|
|
|
|
}
|
|
|
|
|
},
|
|
|
|
|
"nbformat": 4,
|
|
|
|
|
"nbformat_minor": 4
|
|
|
|
|
}
|
