1 files changed, 153 insertions, 0 deletions
diff --git a/lab-windows/gps_clock_jitter_analysis.ipynb b/lab-windows/gps_clock_jitter_analysis.ipynb
new file mode 100644
index 0000000..fbd540a
--- /dev/null
+++ b/lab-windows/gps_clock_jitter_analysis.ipynb
@@ -0,0 +1,153 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import math\n",
+    "import sqlite3\n",
+    "import struct\n",
+    "import datetime\n",
+    "import scipy.fftpack\n",
+    "from scipy import signal as sig\n",
+    "from scipy import optimize as opt\n",
+    "\n",
+    "import matplotlib\n",
+    "from matplotlib import pyplot as plt\n",
+    "from matplotlib import patches\n",
+    "import numpy as np\n",
+    "from scipy import signal, optimize\n",
+    "from tqdm.notebook import tnrange, tqdm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "db = sqlite3.connect('data/waveform_1pps_debug.sqlite3')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Run 000: 2020-03-31 16:58:00 - 2020-03-31 16:58:36 (  0:00:36.029,     36512sp)\n",
+      "Run 001: 2020-03-31 16:58:51 - 2020-03-31 17:05:19 (  0:06:27.729,    392608sp)\n",
+      "Run 002: 2020-03-31 17:07:02 - 2020-03-31 17:41:34 (  0:34:32.105,     37024sp)\n",
+      "Run 003: 2020-03-31 18:50:05 - 2020-03-31 18:50:43 (  0:00:37.576,     38048sp)\n",
+      "Run 004: 2020-03-31 18:54:08 - 2020-03-31 19:14:32 (  0:20:24.104,   1239424sp)\n"
+     ]
+    }
+   ],
+   "source": [
+    "for run_id, start, end, count in db.execute('SELECT run_id, MIN(rx_ts), MAX(rx_ts), COUNT(*) FROM measurements GROUP BY run_id'):\n",
+    "    foo = lambda x: datetime.datetime.fromtimestamp(x/1000)\n",
+    "    start, end = foo(start), foo(end)\n",
+    "    print(f'Run {run_id:03d}: {start:%Y-%m-%d %H:%M:%S} - {end:%Y-%m-%d %H:%M:%S} ({str(end-start)[:-3]:>13}, {count*32:>9d}sp)')\n",
+    "last_run, n_records = run_id, count\n",
+    "sampling_rate = 1000.0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "2712791d61b549ecb531f886d88e1d54",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
+     ]
+    }
+   ],
+   "source": [
+    "histogram = np.array(db.execute('SELECT gps_1pps, COUNT(*) FROM measurements WHERE gps_1pps != -1 AND run_id = ? GROUP BY gps_1pps', (last_run,)).fetchall())\n",
+    "hist_plot = histogram.astype(float)[1:-1]\n",
+    "hist_plot[:, 0] *= 2 / 5 * 2\n",
+    "hist_plot[:, 1] /= (1000 / 32)\n",
+    "\n",
+    "f_nom = 19.440e6\n",
+    "\n",
+    "font = {'family' : 'normal',\n",
+    "        'weight' : 'normal',\n",
+    "        'size'   : 10}\n",
+    "matplotlib.rc('font', **font)\n",
+    "fig, ax = plt.subplots(figsize=(5, 4))\n",
+    "ax.grid()\n",
+    "# We have a bug that causes our measurements to occassionally be out by +/- 65534 counts. For now, fix this by simply throwing away these (very obviously invalid) bins.\n",
+    "ax.bar(hist_plot[:,0] - f_nom , hist_plot[:, 1])\n",
+    "\n",
+    "def gauss(x, *p):\n",
+    "    A, mu, sigma = p\n",
+    "    return A*np.exp(-(x-mu)**2/(2.*sigma**2))\n",
+    "\n",
+    "gauss_x = np.linspace(np.min(hist_plot[:,0]), np.max(hist_plot[:,0]), 10000)\n",
+    "coeff, var_matrix = opt.curve_fit(gauss, hist_plot[:,0], hist_plot[:,1], p0=[np.max(hist_plot[:,1]), np.mean(hist_plot[:,0]), 1])\n",
+    "hist_fit = gauss(gauss_x, *coeff)\n",
+    "ax.plot(gauss_x - f_nom, hist_fit, color='orange')\n",
+    "_A, mu, sigma = coeff\n",
+    "bbox_props = dict(fc='white', alpha=0.8, ec='none')\n",
+    "ax.annotate(f'σ² = {sigma**2 * 1e3:.1f} mHz ({sigma**2 / f_nom * 1e9:.2f} ppb)\\nμ = {mu-f_nom:+.1f} Hz ({(mu-f_nom)/f_nom * 1e6:+.2f} ppm)', xy=[0.6, 0.5], xycoords='figure fraction', bbox=bbox_props)\n",
+    "ax.set_xlabel('$f - f_{nom}$ [Hz]')\n",
+    "ax.set_ylabel('# observations')\n",
+    "\n",
+    "#ax.set_title('OCXO frequency derivation relative to GPS 1pps')\n",
+    "fig.savefig('fig_out/ocxo_freq_stability.eps', format='eps')"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "labenv",
+   "language": "python",
+   "name": "labenv"
+  },
+  "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.8.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}