Add some graphs, add frequency spectra comparison

compare between commercial measurements from Dr. Gobmaier GmbH and ours.
Turns out we agree!

3 files changed, 1214 insertions, 26 deletions

diff --git a/lab-windows/dsss_experiments-ber.ipynb b/lab-windows/dsss_experiments-ber.ipynb
index f490f36..d4a2992 100644
--- a/lab-windows/dsss_experiments-ber.ipynb
+++ b/lab-windows/dsss_experiments-ber.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -47,7 +47,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 37,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -73,7 +73,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 38,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -90,7 +90,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 39,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -104,7 +104,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 40,
    "metadata": {},
    "outputs": [
     {
@@ -124,12 +124,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 77,
    "metadata": {},
    "outputs": [],
    "source": [
-    "def generate_test_signal(duration, nbits=6, signal_amplitude=2.0e-3, decimation=10, seed=0):\n",
-    "    test_data = np.random.RandomState(seed=seed).randint(0, 2 * (2**nbits), duration)\n",
+    "def generate_test_signal(duration, nbits=6, signal_amplitude=2.0e-3, decimation=10, seed=0, data=None):\n",
+    "    test_data = np.random.RandomState(seed=seed).randint(0, 2 * (2**nbits), duration) if data is None else data\n",
     "    \n",
     "    signal = np.repeat(modulate(test_data, nbits) * 2.0 - 1, decimation) * signal_amplitude\n",
     "    noise = np.resize(mains_noise, len(signal))\n",
@@ -139,7 +139,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 42,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -153,13 +153,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 43,
    "metadata": {},
    "outputs": [
     {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "<ipython-input-42-882fdfbdc9fa>:4: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
+      "  fig, ax = plt.subplots()\n"
+     ]
+    },
+    {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "ac41374377284feaa0d8dc20bcb4a341",
+       "model_id": "cddadc4204f54789905875fcfb8e4126",
        "version_major": 2,
        "version_minor": 0
       },
@@ -177,7 +185,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 44,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -186,7 +194,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 45,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -366,13 +374,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 217,
    "metadata": {},
    "outputs": [
     {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "<ipython-input-217-c761a61291a7>:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
+      "  fig, ax = plt.subplots(figsize=(12, 9))\n"
+     ]
+    },
+    {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "b4952032e6ed4d62b23ddf7889a0efb9",
+       "model_id": "42da8a07c0e241f2b1d0f3358e68935b",
        "version_major": 2,
        "version_minor": 0
       },
@@ -382,19 +398,57 @@
      },
      "metadata": {},
      "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "[<matplotlib.lines.Line2D at 0x7efd98cd17f0>]"
+      ]
+     },
+     "execution_count": 217,
+     "metadata": {},
+     "output_type": "execute_result"
     }
    ],
    "source": [
     "fig, ax = plt.subplots(figsize=(12, 9))\n",
     "\n",
-    "params = dict(\n",
-    "        nbits=8,\n",
-    "        signal_amplitude=0.00015,\n",
-    "        decimation=10,\n",
-    "        threshold_factor=6.0,\n",
-    "        power_avg_width=2.5,\n",
-    "        max_lookahead=6.5)\n",
+    "decimation = 100\n",
+    "extra_dec = 10\n",
+    "nbits = 5\n",
+    "signal_amplitude=3e-3,\n",
+    "\n",
+    "seed = 42\n",
+    "\n",
+    "test_data, signal = generate_test_signal(duration, nbits, signal_amplitude, decimation, seed, data=np.array([4,5,6,7] * 8))\n",
+    "\n",
+    "#sosh = sig.butter(6, 1/(2**nbits), btype='highpass', output='sos', fs=decimation)\n",
+    "#sosl = sig.butter(6, 1.0, btype='lowpass', output='sos', fs=decimation)\n",
+    "#filtered = sig.sosfilt(sosh, sig.sosfilt(sosl, signal))\n",
+    "#filtered = sig.sosfilt(sosh, signal)\n",
+    "\n",
+    "cor_an1 = correlate(signal, nbits=nbits, decimation=decimation)\n",
+    "#cor_an2 = correlate(filtered, nbits=nbits, decimation=decimation)\n",
+    "#cor_an2 = correlate(sig.decimate(signal, 9), nbits=nbits, decimation=decimation//9)\n",
+    "cor_an2 = correlate(sig.decimate(signal, extra_dec), nbits=nbits, decimation=int(round(decimation/extra_dec)))\n",
+    "#ax.plot(cor_an[2])\n",
+    "#ax.matshow(sig.cwt(cor_an[2], sig.ricker, np.arange(1, 64)), aspect='auto')\n",
+    "cwt_ed1 = sig.cwt(cor_an1[2], sig.ricker, np.arange(1, 130))\n",
+    "cwt_ed2 = sig.cwt(cor_an2[2], sig.ricker, np.arange(1, 130))\n",
+    "#for f in [0.73, 1.0]:\n",
+    "#    ax.plot(cwt_ed[int(round(f * decimation))], label=f'{f}')\n",
+    "\n",
+    "#ax.plot(signal)\n",
+    "#ax.twiny().plot(sig.decimate(signal, 9), color='orange')\n",
+    "#ax.twiny().plot(sig.decimate(signal, 10), color='orange')\n",
+    "\n",
+    "#ax.matshow(cwt_ed2, aspect='auto')\n",
+    "ax.plot(cwt_ed1[int(round(0.73 * decimation))], label=f'unfiltered')\n",
+    "ax.twinx().twiny().plot(cwt_ed2[int(round(0.73 * decimation / extra_dec))], label=f'filtered', color='orange')\n",
     "\n",
+    "#ax.legend()\n",
+    "#ax.matshow(cor_an[2:4], aspect='auto')\n",
+    "#ser, br = run_ser_test(**params, sample_duration=100, print=print, seed=seed, ax=ax) #, debug_range=(16100, 16700))\n",
     "#seed = 0xcbb3b8cf\n",
     "#for seed in range(10):\n",
     "#    ser, br = run_ser_test(**params, sample_duration=32, print=print, seed=seed) #, debug_range=(16100, 16700))\n",
@@ -612,21 +666,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 103,
+   "execution_count": 36,
    "metadata": {},
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "<ipython-input-103-3e101b35c97a>:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
+      "<ipython-input-36-8e813d331cd8>:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
       "  fig, ((ax, cbar_ax), (intercept_ax, empty)) = plt.subplots(2, 2, figsize=(12, 9), gridspec_kw={'width_ratios': [1, 0.05]})\n"
      ]
     },
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "84e1acbb31914f5f8b3a031654fd7733",
+       "model_id": "c05f590a7c0040628c5c5d4032bb7c8e",
        "version_major": 2,
        "version_minor": 0
       },
@@ -805,7 +859,7 @@
     "    std = data[:,1]\n",
     "    \n",
     "    ax.set_xlim([min(x), max(x)])\n",
-    "    l = ax.plot(x, y, label='SER @ a=0.5', color='orange')\n",
+    "    l = ax.plot(x, y, label='Amplitude at SER=0.5', color='orange')\n",
     "    \n",
     "    x, y, std = zip(*[ (le_x, le_y, le_std) for le_x, le_y, le_std in zip(x, y, std) if le_y is not None ])\n",
     "    y, std = np.array(y), np.array(std)\n",
@@ -815,6 +869,7 @@
     "    ax.fill_between([-1, min(ser_valid)], 0, 1, facecolor='red', alpha=0.2, transform=trans, zorder=1)\n",
     "    ax.fill_between([max(ser_valid), max(ser_valid)*10], 0, 1, facecolor='red', alpha=0.2, transform=trans)\n",
     "    ax.set_ylim([min(y)*0.9, max(y)*1.1])\n",
+    "    ax.grid()\n",
     "    return l\n",
     "\n",
     "l1 = plot_intercepts(intercept_ax)\n",
@@ -835,6 +890,124 @@
    "source": [
     "#fig.savefig('dsss_prototype_symbol_error_rate_5-8_bit.svg')"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 35,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "<ipython-input-35-cafaa6062c72>:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
+      "  fig, ((ax, cbar_ax), (intercept_ax, empty)) = plt.subplots(2, 2, figsize=(12, 9), gridspec_kw={'width_ratios': [1, 0.05]})\n"
+     ]
+    },
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "fa93fb31355840148f941dab8b60f51c",
+       "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": [
+      "<ipython-input-35-cafaa6062c72>:16: RuntimeWarning: divide by zero encountered in log10\n",
+      "  cm_func = lambda x: cmap(np.log10(x - min(decimations)) / (np.log10(max(decimations)) - np.log10(min(decimations))))\n"
+     ]
+    }
+   ],
+   "source": [
+    "fig, ((ax, cbar_ax), (intercept_ax, empty)) = plt.subplots(2, 2, figsize=(12, 9), gridspec_kw={'width_ratios': [1, 0.05]})\n",
+    "empty.axis('off')\n",
+    "#fig.tight_layout()\n",
+    "\n",
+    "results = []\n",
+    "for fn in [\n",
+    "    '/mnt/c/Users/jaseg/shared/dsss_experiments_res-2020-02-20-14-10-13.json',\n",
+    "    '/mnt/c/Users/jaseg/shared/dsss_experiments_res-2020-02-20-13-21-57.json',\n",
+    "    '/mnt/c/Users/jaseg/shared/dsss_experiments_res-2020-02-20-13-23-47.json',\n",
+    "]:\n",
+    "    with open(fn, 'r') as f:\n",
+    "        results += json.load(f)\n",
+    "\n",
+    "decimations = [decimation for (_nbits, thf, _reps, _points, _duration, decimation), series in results]\n",
+    "cmap = matplotlib.cm.viridis\n",
+    "cm_func = lambda x: cmap(np.log10(x - min(decimations)) / (np.log10(max(decimations)) - np.log10(min(decimations))))\n",
+    "\n",
+    "decimation_sers = {}\n",
+    "for (nbits, thf, reps, points, duration, decimation), series in results:\n",
+    "    data = [ [ mean for mean, _std, _msg in reps if mean is not None ] for _amp, reps in series ]\n",
+    "    amps = [ amp for amp, _reps in series ]\n",
+    "    sers = np.array([ np.mean(values) for values in data ])\n",
+    "    stds = np.array([ np.std(values)  for values in data ])\n",
+    "    decimation_sers[decimation] = list(zip(amps, sers, stds))\n",
+    "    \n",
+    "    l, = ax.plot(amps, np.clip(sers, 0, 1), label=f'decimation={decimation}', color=cm_func(decimation))\n",
+    "    ax.fill_between(amps, np.clip(sers + stds, 0, 1), np.clip(sers - stds, 0, 1), facecolor=l.get_color(), alpha=0.2)\n",
+    "    ax.axhline(0.5, color='gray', ls=(0, (3, 4)), lw=0.8)\n",
+    "ax.grid()\n",
+    "ax.set_xlabel('Amplitude in mHz')\n",
+    "ax.set_ylabel('Symbol error rate')\n",
+    "\n",
+    "norm = matplotlib.colors.Normalize(vmin=np.log10(min(decimations)), vmax=np.log10(max(decimations)))\n",
+    "cb1 = matplotlib.colorbar.ColorbarBase(cbar_ax, cmap=cmap, norm=norm, orientation='vertical', label=\"Decimation\", ticks=[np.log10(d) for d in decimations])\n",
+    "cb1.ax.set_yticklabels([f'{d:.1f}' for d in decimations])\n",
+    "\n",
+    "def plot_intercepts(ax, SER_TH = 0.5):\n",
+    "    intercepts = {}\n",
+    "    for dec, sers in decimation_sers.items():\n",
+    "        last_ser, last_amp, last_std = 0, 0, 0\n",
+    "        for amp, ser, std in sorted(sers):\n",
+    "            if last_ser > SER_TH and ser < SER_TH:\n",
+    "                icp = last_amp + (SER_TH - last_ser) / (ser - last_ser) * (amp - last_amp)\n",
+    "                ic_std = abs(last_amp - amp) / 2# np.sqrt(np.mean(last_std**2 + std**2))\n",
+    "                intercepts[dec] = (icp, ic_std)\n",
+    "                break\n",
+    "            last_amp, last_ser = amp, ser\n",
+    "        else:\n",
+    "            intercepts[dec] = None, None\n",
+    "    \n",
+    "    ser_valid = [dec for dec, (ser, _std) in intercepts.items() if ser is not None]\n",
+    "    #ax.axvline(min(ser_valid), color='red')\n",
+    "    #ax.axvline(max(ser_valid), color='red')\n",
+    "    \n",
+    "    x = sorted(intercepts.keys())\n",
+    "    data = np.array([ intercepts[dec] for dec in x ])\n",
+    "    y = data[:,0]\n",
+    "    std = data[:,1]\n",
+    "    \n",
+    "    ax.set_xlim([min(x), max(x)])\n",
+    "    l = ax.plot(x, y, label='Amplitude at SER=0.5', color='orange')\n",
+    "    ax.legend(loc=3)\n",
+    "    ax.grid()\n",
+    "    \n",
+    "    x, y, std = zip(*[ (le_x, le_y, le_std) for le_x, le_y, le_std in zip(x, y, std) if le_y is not None ])\n",
+    "    y, std = np.array(y), np.array(std)\n",
+    "    ax.fill_between(x, y-std, y+std, color=l[0].get_color(), alpha=0.3)\n",
+    "    \n",
+    "    trans = matplotlib.transforms.blended_transform_factory(ax.transData, ax.transAxes)\n",
+    "    ax.fill_between([-1, min(ser_valid)], 0, 1, facecolor='red', alpha=0.2, transform=trans, zorder=1)\n",
+    "    ax.fill_between([max(ser_valid), max(ser_valid)*10], 0, 1, facecolor='red', alpha=0.2, transform=trans)\n",
+    "    ax.set_ylim([min(y)*0.9, max(y)*1.1])\n",
+    "    ax.set_xscale('log')\n",
+    "    ax.xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(lambda x, _: '{:g}'.format(x)))\n",
+    "    ax.set_xticks([1, 2, 5, 10, 20, 50])\n",
+    "    ax.set_xlim([1, 60])\n",
+    "    return l\n",
+    "\n",
+    "l1 = plot_intercepts(intercept_ax)\n"
+   ]
   }
  ],
  "metadata": {
diff --git a/lab-windows/grid_frequency_spectra.ipynb b/lab-windows/grid_frequency_spectra.ipynb
new file mode 100644
index 0000000..3d375ef
--- /dev/null
+++ b/lab-windows/grid_frequency_spectra.ipynb
@@ -0,0 +1,340 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import csv\n",
+    "\n",
+    "import numpy as np\n",
+    "from matplotlib import pyplot as plt\n",
+    "import scipy.fftpack"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data = np.genfromtxt('data/Netzfrequenz_Sekundenwerte_2012_KW37.csv', delimiter=',')[1:,1:]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "36f1f4d7970e41afa4737b6f63b7c449",
+       "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"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "[<matplotlib.lines.Line2D at 0x7f952a6e4580>]"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(data[:3600*24, 0])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.02051102806199375"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.std(data[:,0])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "d7fe0512f4254efeb15235a5617ef064",
+       "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"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "(1e-06, 0.5)"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# Number of samplepoints\n",
+    "N = len(data[:,0])\n",
+    "# sample spacing\n",
+    "T = 1.0\n",
+    "x = np.linspace(0.0, N*T, N)\n",
+    "yf = scipy.fftpack.fft(data[:,0])\n",
+    "xf = np.linspace(0.0, 1.0/(2.0*T), N//2)\n",
+    "\n",
+    "yf = 2.0/N * np.abs(yf[:N//2])\n",
+    "\n",
+    "#yf = sum(yf[s::10] for s in range(10)) / 10\n",
+    "#xf = sum(xf[s::10] for s in range(10)) / 10\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.loglog(xf, yf)\n",
+    "ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: f'{1/x:.1f}'))\n",
+    "ax.set_xlabel('T in s')\n",
+    "ax.set_ylabel('Amplitude Δf')\n",
+    "ax.grid()\n",
+    "ax.set_xlim([1/1000000, 0.5])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "fd94bc97a276400db0539b703c4eeeac",
+       "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"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "(1.6666666666666667e-05, 0.5)"
+      ]
+     },
+     "execution_count": 23,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# Number of samplepoints\n",
+    "N = len(data[:,0])\n",
+    "# sample spacing\n",
+    "T = 1.0\n",
+    "x = np.linspace(0.0, N*T, N)\n",
+    "yf = scipy.fftpack.fft(data[:,0])\n",
+    "xf = np.linspace(0.0, 1.0/(2.0*T), N//2)\n",
+    "\n",
+    "yf = 2.0/N * np.abs(yf[:N//2])\n",
+    "\n",
+    "average_from = lambda val, start, average_width: np.hstack([val[:start], [ np.mean(val[i:i+average_width]) for i in range(start, len(val), average_width) ]])\n",
+    "\n",
+    "average_width = 20\n",
+    "average_start = 100\n",
+    "yf = average_from(yf, average_start, average_width)\n",
+    "xf = average_from(xf, average_start, average_width)\n",
+    "yf = average_from(yf, 300, average_width)\n",
+    "xf = average_from(xf, 300, average_width)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.loglog(xf, yf)\n",
+    "ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: f'{1/x:.1f}'))\n",
+    "ax.set_xlabel('T in s')\n",
+    "ax.set_ylabel('Amplitude Δf')\n",
+    "ax.grid()\n",
+    "ax.set_xlim([1/60000, 0.5])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "TypeError",
+     "evalue": "object of type <class 'float'> cannot be safely interpreted as an integer.",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
+      "\u001b[0;32m~/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/function_base.py\u001b[0m in \u001b[0;36mlinspace\u001b[0;34m(start, stop, num, endpoint, retstep, dtype, axis)\u001b[0m\n\u001b[1;32m    116\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 117\u001b[0;31m         \u001b[0mnum\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moperator\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    118\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mTypeError\u001b[0m: 'float' object cannot be interpreted as an integer",
+      "\nDuring handling of the above exception, another exception occurred:\n",
+      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
+      "\u001b[0;32m<ipython-input-7-75728c9461c4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mys\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconvolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mys\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mones\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'valid'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mxs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.0\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      8\u001b[0m \u001b[0;31m#xs = np.linspace(len(data)/2, 1, len(data)/2)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36mlinspace\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
+      "\u001b[0;32m~/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/function_base.py\u001b[0m in \u001b[0;36mlinspace\u001b[0;34m(start, stop, num, endpoint, retstep, dtype, axis)\u001b[0m\n\u001b[1;32m    117\u001b[0m         \u001b[0mnum\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moperator\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    118\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 119\u001b[0;31m         raise TypeError(\n\u001b[0m\u001b[1;32m    120\u001b[0m             \u001b[0;34m\"object of type {} cannot be safely interpreted as an integer.\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    121\u001b[0m                 .format(type(num)))\n",
+      "\u001b[0;31mTypeError\u001b[0m: object of type <class 'float'> cannot be safely interpreted as an integer."
+     ]
+    }
+   ],
+   "source": [
+    "ys = scipy.fftpack.fft(data[:,0])\n",
+    "ys = 2.0/len(data) * np.abs(ys[:len(data)//2])\n",
+    "s = 60\n",
+    "\n",
+    "ys = np.convolve(ys, np.ones((s,))/s, mode='valid')\n",
+    "\n",
+    "xs = np.linspace(0, 1.0/2.0, len(data)/2)\n",
+    "#xs = np.linspace(len(data)/2, 1, len(data)/2)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.loglog(xs[s//2:-s//2+1], ys)\n",
+    "ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: 1/x))\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ys = scipy.fftpack.fft(data[:,0])\n",
+    "ys = 2.0/len(data) * np.abs(ys[:len(data)//2])\n",
+    "s = 1\n",
+    "\n",
+    "ys = np.convolve(ys, np.ones((s,))/s, mode='valid')\n",
+    "\n",
+    "xs = np.linspace(0, 1.0/2.0, len(data)/2)\n",
+    "#xs = np.linspace(len(data)/2, 1, len(data)/2)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.loglog(xs[s//2:-s//2+1 if s > 1 else None], ys)\n",
+    "ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: 1/x))\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ys = scipy.fftpack.fft(data[:,0])\n",
+    "ys = 2.0/len(data) * np.abs(ys[:len(data)//2])\n",
+    "s = 1\n",
+    "\n",
+    "ys = np.convolve(ys, np.ones((s,))/s, mode='valid')\n",
+    "\n",
+    "xs = np.linspace(0, 1.0/2.0, len(data)/2)\n",
+    "\n",
+    "ys *= 2*np.pi*xs\n",
+    "#xs = np.linspace(len(data)/2, 1, len(data)/2)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.loglog(xs[s//2:-s//2+1 if s > 1 else None], ys)\n",
+    "ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: 1/x))\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ys = scipy.fftpack.fft(data[:,0])\n",
+    "ys = 2.0/len(data) * np.abs(ys[:len(data)//2])\n",
+    "s = 30\n",
+    "\n",
+    "ys = np.convolve(ys, np.ones((s,))/s, mode='valid')\n",
+    "\n",
+    "xs = np.linspace(0, 1.0/2.0, len(data)/2)\n",
+    "\n",
+    "ys *= 2*np.pi*xs[s//2:-s//2+1]\n",
+    "\n",
+    "#xs = np.linspace(len(data)/2, 1, len(data)/2)\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(9,5))\n",
+    "ax.loglog(xs[s//2:-s//2+1], ys)\n",
+    "ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: 1/x))\n",
+    "ax.grid()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "1/0.0628"
+   ]
+  }
+ ],
+ "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.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/lab-windows/grid_scope.ipynb b/lab-windows/grid_scope.ipynb
new file mode 100644
index 0000000..9f53906
--- /dev/null
+++ b/lab-windows/grid_scope.ipynb
@@ -0,0 +1,675 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 104,
+   "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",
+    "\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('/mnt/c/Users/jaseg/shared/waveform-raspi.sqlite3')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Run 000: 2020-01-31 18:05:24 - 2020-02-01 00:13:45 (  6:08:21.589,  22126080sp)\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": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "par = lambda *rs: 1/sum(1/r for r in rs) # resistor parallel calculation\n",
+    "\n",
+    "# FIXME: These are for the first prototype only!\n",
+    "vmeas_source_impedance = 330e3\n",
+    "vmeas_source_scale = 0.5\n",
+    "\n",
+    "vcc = 15.0\n",
+    "vmeas_div_high = 27e3\n",
+    "vmeas_div_low = par(4.7e3, 10e3)\n",
+    "vmeas_div_voltage = vcc * vmeas_div_low / (vmeas_div_high + vmeas_div_low)\n",
+    "vmeas_div_impedance = par(vmeas_div_high, vmeas_div_low)\n",
+    "\n",
+    "#vmeas_overall_factor = vmeas_div_impedance / (vmeas_source_impedance + vmeas_div_impedance)\n",
+    "v0 = 1.5746\n",
+    "v100 = 2.004\n",
+    "vn100 = 1.1452\n",
+    "\n",
+    "adc_vcc = 3.3 # V\n",
+    "adc_fullscale = 4095\n",
+    "\n",
+    "adc_val_to_voltage_factor = 1/adc_fullscale * adc_vcc\n",
+    "\n",
+    "adc_count_to_vmeas = lambda x: (x*adc_val_to_voltage_factor - v0) / (v100-v0) * 100"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "f960454ab93244db97e039ae7ebd02ee",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "HBox(children=(FloatProgress(value=0.0, max=691440.0), HTML(value='')))"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "limit = n_records\n",
+    "record_size = 32\n",
+    "skip_dropped_sections = False\n",
+    "\n",
+    "data = np.zeros(limit*record_size)\n",
+    "data[:] = np.nan\n",
+    "\n",
+    "last_seq = None\n",
+    "write_index = 0\n",
+    "for i, (seq, chunk) in tqdm(enumerate(db.execute(\n",
+    "        'SELECT seq, data FROM measurements WHERE run_id = ? ORDER BY rx_ts LIMIT ? OFFSET ?',\n",
+    "        (last_run, limit, n_records-limit))), total=n_records):\n",
+    "    \n",
+    "    if last_seq is None or seq == (last_seq + 1)%0xffff:\n",
+    "        last_seq = seq\n",
+    "        idx = write_index if skip_dropped_sections else i\n",
+    "        data[idx*record_size:(idx+1)*record_size] = np.frombuffer(chunk, dtype='<H')\n",
+    "        write_index += 1\n",
+    "        \n",
+    "    elif seq > last_seq:\n",
+    "        last_seq = seq\n",
+    "        # nans = np.empty((record_size,))\n",
+    "        # nans[:] = np.nan\n",
+    "        # data = np.append(data, nans) FIXME\n",
+    "        \n",
+    "data = (data * adc_val_to_voltage_factor - v0) / (v100-v0) * 100"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "227.68691180713367"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "data_not_nan = data[~np.isnan(data)]\n",
+    "np.sqrt(np.mean(np.square(data_not_nan)))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "264a9f8478e449289c1592c9595dc6cc",
+       "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"
+    }
+   ],
+   "source": [
+    "fig, (top, bottom) = plt.subplots(2, figsize=(9,6))\n",
+    "fig.tight_layout(pad=3, h_pad=0.1)\n",
+    "\n",
+    "range_start, range_len = -300, 60 # [s]\n",
+    "\n",
+    "data_slice = data[ int(range_start * sampling_rate) : int((range_start + range_len) * sampling_rate) ]\n",
+    "\n",
+    "top.grid()\n",
+    "top.plot(np.linspace(0, range_len, int(range_len*sampling_rate)), data_slice, lw=1.0)\n",
+    "top.set_xlim([range_len/2-0.25, range_len/2+0.25])\n",
+    "mean = np.mean(data_not_nan)\n",
+    "rms = np.sqrt(np.mean(np.square(data_not_nan - mean)))\n",
+    "peak = np.max(np.abs(data_not_nan - mean))\n",
+    "top.axhline(mean, color='red')\n",
+    "bbox = {'facecolor': 'black', 'alpha': 0.8, 'pad': 2}\n",
+    "top.text(0, mean, f'mean: {mean:.3f}', color='white', bbox=bbox)\n",
+    "top.text(0.98, 0.2, f'V_RMS: {rms:.3f}', transform=top.transAxes, color='white', bbox=bbox, ha='right')\n",
+    "top.text(0.98, 0.1, f'V_Pk: {peak:.3f}', transform=top.transAxes, color='white', bbox=bbox, ha='right')\n",
+    "\n",
+    "bottom.grid()\n",
+    "bottom.specgram(data_slice, Fs=sampling_rate)\n",
+    "top.set_ylabel('U [V]')\n",
+    "bottom.set_ylabel('F [Hz]')\n",
+    "bottom.set_xlabel('t [s]')\n",
+    "None"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fs = sampling_rate # Hz\n",
+    "ff = 50 # Hz\n",
+    "\n",
+    "analysis_periods = 10\n",
+    "window_len = fs * analysis_periods/ff\n",
+    "nfft_factor = 4\n",
+    "sigma = window_len/8 # samples\n",
+    "\n",
+    "f, t, Zxx = signal.stft(data,\n",
+    "            fs = fs,\n",
+    "            window=('gaussian', sigma),\n",
+    "            nperseg = window_len,\n",
+    "            nfft = window_len * nfft_factor)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "7c2382eb8e124ef9b29546ecaed3e155",
+       "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"
+    }
+   ],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(9, 3))\n",
+    "fig.tight_layout(pad=2, h_pad=0.1)\n",
+    "\n",
+    "ax.pcolormesh(t[-200:-100], f[:250], np.abs(Zxx[:250,-200:-100]))\n",
+    "ax.set_title(f\"Run {last_run}\", pad=-20, color='white')\n",
+    "ax.grid()\n",
+    "ax.set_ylabel('f [Hz]')\n",
+    "ax.set_ylim([30, 75]) # Hz\n",
+    "ax.set_xlabel('simulation time t [s]')\n",
+    "None"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "d29e8ee5bf7f4e94a1749239aec24d6d",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "HBox(children=(FloatProgress(value=0.0, max=221260.0), HTML(value='')))"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "f_t = t\n",
+    "\n",
+    "n_f, n_t = Zxx.shape\n",
+    "# start, stop = 180, 220\n",
+    "# start, stop = 90, 110\n",
+    "# start, stop = 15, 35\n",
+    "# bounds_f = slice(start // 4 * nfft_factor, stop // 4 * nfft_factor)\n",
+    "f_min, f_max = 30, 70 # Hz\n",
+    "bounds_f = slice(np.argmax(f > f_min), np.argmin(f < f_max))\n",
+    "\n",
+    "\n",
+    "f_mean = np.zeros(Zxx.shape[1])\n",
+    "for le_t in tnrange(1, Zxx.shape[1] - 1):\n",
+    "    frame_f = f[bounds_f]\n",
+    "    frame_step = frame_f[1] - frame_f[0]\n",
+    "    time_step = f_t[1] - f_t[0]\n",
+    "    #if t == 10:\n",
+    "    #    axs[-1].plot(frame_f, frame_Z)\n",
+    "    frame_Z = np.abs(Zxx[bounds_f, le_t])\n",
+    "    # frame_f = f[180:220]\n",
+    "    # frame_Z = np.abs(Zxx[180:220, 40])\n",
+    "    # frame_f = f[15:35]\n",
+    "    # frame_Z = np.abs(Zxx[15:35, 40])\n",
+    "    # plt.plot(frame_f, frame_Z)\n",
+    "\n",
+    "    # peak_f = frame_f[np.argmax(frame)]\n",
+    "    # plt.axvline(peak_f, color='red')\n",
+    "\n",
+    "#         def gauss(x, *p):\n",
+    "#             A, mu, sigma, o = p\n",
+    "#             return A*np.exp(-(x-mu)**2/(2.*sigma**2)) + o\n",
+    "\n",
+    "    def gauss(x, *p):\n",
+    "        A, mu, sigma = p\n",
+    "        return A*np.exp(-(x-mu)**2/(2.*sigma**2))\n",
+    "\n",
+    "    f_start = frame_f[np.argmax(frame_Z)]\n",
+    "    A_start = np.max(frame_Z)\n",
+    "    p0 = [A_start, f_start, 1.]\n",
+    "    try:\n",
+    "        coeff, var = optimize.curve_fit(gauss, frame_f, frame_Z, p0=p0)\n",
+    "        # plt.plot(frame_f, gauss(frame_f, *coeff))\n",
+    "        #print(coeff)\n",
+    "        A, mu, sigma, *_ = coeff\n",
+    "        f_mean[le_t] = mu\n",
+    "    except Exception:\n",
+    "        f_mean[le_t] = np.nan"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "47c40f28b5a34a94b4a631fd62107521",
+       "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"
+    }
+   ],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(9, 5), sharex=True)\n",
+    "fig.tight_layout(pad=2.2, h_pad=0, w_pad=1)\n",
+    "\n",
+    "label = f'Run {last_run}'\n",
+    "ax.plot(f_t[1:-1], f_mean[1:-1])\n",
+    "\n",
+    "#     b, a = signal.butter(3,\n",
+    "#                          1/5, # Hz\n",
+    "#                          btype='lowpass',\n",
+    "#                          fs=1/time_step)\n",
+    "#     filtered = signal.lfilter(b, a, f_mean[1:-1], axis=0)\n",
+    "#     ax.plot(f_t[1:-1], filtered)\n",
+    "\n",
+    "ax.set_title(label, pad=-20)\n",
+    "ax.set_ylabel('f [Hz]')\n",
+    "ax.grid()\n",
+    "if not label in ['off_frequency', 'sweep_phase_steps']:\n",
+    "    ax.set_ylim([49.90, 50.10])\n",
+    "    var = np.var(f_mean[~np.isnan(f_mean)][1:-1])\n",
+    "    ax.text(0.5, 0.08, f'σ²={var * 1e3:.3g} mHz²', transform=ax.transAxes, ha='center', color='white', bbox=bbox)\n",
+    "    ax.text(0.5, 0.15, f'σ={np.sqrt(var) * 1e3:.3g} mHz', transform=ax.transAxes, ha='center', color='white', bbox=bbox)\n",
+    "\n",
+    "#         ax.text(0.5, 0.2, f'filt. σ²={np.var(filtered) * 1e3:.3g} mHz', transform=ax.transAxes, ha='center')\n",
+    "else:\n",
+    "    f_min, f_max = min(f_mean[1:-1]), max(f_mean[1:-1])\n",
+    "    delta = f_max - f_min\n",
+    "    ax.set_ylim(f_min - delta * 0.1, f_max + delta * 0.3)\n",
+    "\n",
+    "for i in np.where(np.isnan(f_mean))[0]:\n",
+    "    ax.axvspan(f_t[i], f_t[i+1], color='lightblue')\n",
+    "\n",
+    "ax.set_xlabel('recording time t [s]')\n",
+    "None"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "4a4cb62296df496bad37d93547d3c26a",
+       "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"
+    }
+   ],
+   "source": [
+    "f_copy = np.copy(f_mean[1:-1])\n",
+    "f_copy[np.isnan(f_copy)] = np.mean(f_copy[~np.isnan(f_copy)])\n",
+    "b, a = signal.cheby2(7, 86, 100, 'low', output='ba', fs=1000)\n",
+    "filtered = signal.lfilter(b, a, f_copy)\n",
+    "\n",
+    "b2, a2 = signal.cheby2(3, 30, 1, 'high', output='ba', fs=1000)\n",
+    "filtered2 = signal.lfilter(b2, a2, filtered)\n",
+    "\n",
+    "fig, (ax2, ax1) = plt.subplots(2, figsize=(9,7))\n",
+    "ax1.plot(f_t[1:-1], f_copy, color='lightgray')\n",
+    "ax1.set_ylim([49.90, 50.10])\n",
+    "ax1.grid()\n",
+    "formatter = matplotlib.ticker.FuncFormatter(lambda s, x: str(datetime.timedelta(seconds=s)))\n",
+    "ax1.xaxis.set_major_formatter(formatter)\n",
+    "zoom_offx = 7000 # s\n",
+    "zoom_len = 300 # s\n",
+    "ax1.set_xlim([zoom_offx, zoom_offx + zoom_len])\n",
+    "\n",
+    "ax1.plot(f_t[1:-1], filtered, color='orange')\n",
+    "ax1r = ax1.twinx()\n",
+    "ax1r.plot(f_t[1:-1], filtered2, color='red')\n",
+    "ax1r.set_ylim([-0.015, 0.015])\n",
+    "ax1.set_title(f'Zoomed trace ({datetime.timedelta(seconds=zoom_len)})', pad=-20)\n",
+    "\n",
+    "\n",
+    "ax2.set_title(f'Run {last_run}')\n",
+    "ax2.plot(f_t[1:-1], f_copy, color='orange')\n",
+    "\n",
+    "ax2r = ax2.twinx()\n",
+    "ax2r.set_ylim([-0.1, 0.1])\n",
+    "ax2r.plot(f_t[1:-1], filtered2, color='red')\n",
+    "#ax2.plot(f_t[1:-1], filtered, color='orange', zorder=1)\n",
+    "ax2.set_ylim([49.90, 50.10])\n",
+    "ax2.set_xlim([0, f_t[-2]])\n",
+    "ax2.grid()\n",
+    "formatter = matplotlib.ticker.FuncFormatter(lambda s, x: str(datetime.timedelta(seconds=s)))\n",
+    "ax2.xaxis.set_major_formatter(formatter)\n",
+    "\n",
+    "ax2.legend(handles=[\n",
+    "    patches.Patch(color='lightgray', label='Raw frequency'),\n",
+    "    patches.Patch(color='orange', label='low-pass filtered'),\n",
+    "    patches.Patch(color='red', label='band-pass filtered')])\n",
+    "\n",
+    "#ax2r.spines['right'].set_color('red')\n",
+    "ax2r.yaxis.label.set_color('red')\n",
+    "#ax2r.tick_params(axis='y', colors='red')\n",
+    "\n",
+    "#ax1r.spines['right'].set_color('red')\n",
+    "ax1r.yaxis.label.set_color('red')\n",
+    "#ax1r.tick_params(axis='y', colors='red')\n",
+    "\n",
+    "ax1.set_ylabel('f [Hz]')\n",
+    "ax1r.set_ylabel('band-pass Δf [Hz]')\n",
+    "ax2.set_ylabel('f [Hz]')\n",
+    "ax2r.set_ylabel('band-pass Δf [Hz]')\n",
+    "\n",
+    "# Cut out first 10min of filtered data to give filters time to settle\n",
+    "rms_slice = filtered2[np.where(f_t[1:] > 10*60)[0][0]:]\n",
+    "rms = np.sqrt(np.mean(np.square(rms_slice)))\n",
+    "ax1.text(0.5, 0.1, f'RMS (band-pass): {rms*1e3:.3f}mHz', transform=ax1.transAxes, color='white', bbox=bbox, ha='center')\n",
+    "None"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "chunk_size = 256\n",
+    "#\n",
+    "#with open('filtered_freq.bin', 'wb') as f:\n",
+    "#    for chunk in range(0, len(rms_slice), chunk_size):\n",
+    "#        out_data = rms_slice[chunk:chunk+chunk_size]\n",
+    "#        f.write(struct.pack(f'{len(out_data)}f',  *out_data))\n",
+    "#        \n",
+    "#with open('raw_freq.bin', 'wb') as f:\n",
+    "#    for chunk in range(0, len(f_copy), chunk_size):\n",
+    "#        out_data = f_copy[chunk:chunk+chunk_size]\n",
+    "#        f.write(struct.pack(f'{len(out_data)}f',  *out_data))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 81,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "<ipython-input-81-a751e13723ea>:17: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
+      "  fig, ax = plt.subplots(figsize=(9,5))\n"
+     ]
+    },
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "3809a1a83b5844e3906f1d74bcd15b5d",
+       "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": [
+      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part\n",
+      "  return array(a, dtype, copy=False, order=order)\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "5.0"
+      ]
+     },
+     "execution_count": 81,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "data = f_copy\n",
+    "ys = scipy.fftpack.fft(data)\n",
+    "ys = scipy.fftpack.fftshift(ys)\n",
+    "#ys = 2.0/len(data) * np.abs(ys[:len(data)//2])\n",
+    "#s = 3\n",
+    "\n",
+    "#ys = np.convolve(ys, np.ones((s,))/s, mode='valid')\n",
+    "\n",
+    "#xs = np.linspace(0, 5, len(data)//2)\n",
+    "xs = np.linspace(-5, 5, len(data))\n",
+    "\n",
+    "#ys *= 2*np.pi*xs[s//2:-s//2+1]\n",
+    "#ys *= xs\n",
+    "\n",
+    "#xs = np.linspace(len(data)/2, 1, len(data)/2)\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(9,5))\n",
+    "#ax.loglog(xs[s//2:-s//2+1], ys)\n",
+    "#ax.loglog(xs[s//2:-s//2+1], ys)\n",
+    "#ax.loglog(xs, ys)\n",
+    "#ys[len(xs)//2] = 0\n",
+    "#ax.set_yscale('log')\n",
+    "ax.plot(xs, ys)\n",
+    "#ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: f'{1/x:.1f}'))\n",
+    "ax.grid()\n",
+    "#plt.show()\n",
+    "xs[-1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 156,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "<ipython-input-156-ddde6af5dee1>:20: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
+      "  fig, ax = plt.subplots()\n"
+     ]
+    },
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "d137ae59ed7947ce8e3f7295e102f2f0",
+       "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"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "(1.6666666666666667e-05, 0.5)"
+      ]
+     },
+     "execution_count": 156,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# Number of samplepoints\n",
+    "N = len(data)\n",
+    "# sample spacing\n",
+    "T = 1.0 / 10.0\n",
+    "x = np.linspace(0.0, N*T, N)\n",
+    "yf = scipy.fftpack.fft(data)\n",
+    "xf = np.linspace(0.0, 1.0/(2.0*T), N//2)\n",
+    "\n",
+    "yf = 2.0/N * np.abs(yf[:N//2])\n",
+    "\n",
+    "average_from = lambda val, start, average_width: np.hstack([val[:start], [ np.mean(val[i:i+average_width]) for i in range(start, len(val), average_width) ]])\n",
+    "\n",
+    "average_width = 6\n",
+    "average_start = 20\n",
+    "yf = average_from(yf, average_start, average_width)\n",
+    "xf = average_from(xf, average_start, average_width)\n",
+    "yf = average_from(yf, 200, average_width)\n",
+    "xf = average_from(xf, 200, average_width)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.loglog(xf, yf)\n",
+    "ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: f'{1/x:.1f}'))\n",
+    "ax.set_xlabel('T in s')\n",
+    "ax.set_ylabel('Amplitude Δf')\n",
+    "ax.grid()\n",
+    "ax.set_xlim([1/60000, 0.5])"
+   ]
+  }
+ ],
+ "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.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}