Add some graphs, add frequency spectra comparison

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

1 files changed, 199 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": {