Add BER curve for DSSS experiments

2 files changed, 556 insertions, 35 deletions

diff --git a/lab-windows/dsss_experiments-ber.ipynb b/lab-windows/dsss_experiments-ber.ipynb
new file mode 100644
index 0000000..eda21d7
--- /dev/null
+++ b/lab-windows/dsss_experiments-ber.ipynb
@@ -0,0 +1,505 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 68,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matplotlib import pyplot as plt\n",
+    "import numpy as np\n",
+    "from scipy import signal as sig\n",
+    "import struct\n",
+    "import random\n",
+    "import ipywidgets\n",
+    "import itertools\n",
+    "from multiprocessing import Pool\n",
+    "\n",
+    "import colorednoise\n",
+    "\n",
+    "np.set_printoptions(linewidth=240)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sampling_rate = 10 # sp/s"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 59,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# From https://github.com/mubeta06/python/blob/master/signal_processing/sp/gold.py\n",
+    "preferred_pairs = {5:[[2],[1,2,3]], 6:[[5],[1,4,5]], 7:[[4],[4,5,6]],\n",
+    "                        8:[[1,2,3,6,7],[1,2,7]], 9:[[5],[3,5,6]], \n",
+    "                        10:[[2,5,9],[3,4,6,8,9]], 11:[[9],[3,6,9]]}\n",
+    "\n",
+    "def gen_gold(seq1, seq2):\n",
+    "    gold = [seq1, seq2]\n",
+    "    for shift in range(len(seq1)):\n",
+    "        gold.append(seq1 ^ np.roll(seq2, -shift))\n",
+    "    return gold\n",
+    "\n",
+    "def gold(n):\n",
+    "    n = int(n)\n",
+    "    if not n in preferred_pairs:\n",
+    "        raise KeyError('preferred pairs for %s bits unknown' % str(n))\n",
+    "    t0, t1 = preferred_pairs[n]\n",
+    "    (seq0, _st0), (seq1, _st1) = sig.max_len_seq(n, taps=t0), sig.max_len_seq(n, taps=t1)\n",
+    "    return gen_gold(seq0, seq1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def modulate(data, nbits=5):\n",
+    "    # 0, 1 -> -1, 1\n",
+    "    mask = np.array(gold(nbits))*2 - 1\n",
+    "    \n",
+    "    sel = mask[data>>1]\n",
+    "    data_lsb_centered = ((data&1)*2 - 1)\n",
+    "\n",
+    "    return (np.multiply(sel, np.tile(data_lsb_centered, (2**nbits-1, 1)).T).flatten() + 1) // 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def correlate(sequence, nbits=5, decimation=1, mask_filter=lambda x: x):\n",
+    "    mask = np.tile(np.array(gold(nbits))[:,:,np.newaxis]*2 - 1, (1, 1, decimation)).reshape((2**nbits + 1, (2**nbits-1) * decimation))\n",
+    "\n",
+    "    sequence -= np.mean(sequence)\n",
+    "    \n",
+    "    return np.array([np.correlate(sequence, row, mode='full') for row in mask])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "mean: 49.98625\n"
+     ]
+    }
+   ],
+   "source": [
+    "with open('/mnt/c/Users/jaseg/shared/raw_freq.bin', 'rb') as f:\n",
+    "    mains_noise = np.copy(np.frombuffer(f.read(), dtype='float32'))\n",
+    "    print('mean:', np.mean(mains_noise))\n",
+    "    mains_noise -= np.mean(mains_noise)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 72,
+   "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",
+    "    \n",
+    "    signal = np.repeat(modulate(test_data, nbits) * 2.0 - 1, decimation) * signal_amplitude\n",
+    "    noise = np.resize(mains_noise, len(signal))\n",
+    "    \n",
+    "    return test_data, signal + noise"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nonlinear_distance = lambda x: 100**(2*np.abs(0.5-x%1)) / (np.abs(x)+3)**2\n",
+    "\n",
+    "def plot_distance_func():\n",
+    "    fig, ax = plt.subplots()\n",
+    "    x = np.linspace(-1.5, 5.5, 10000)\n",
+    "    ax.plot(x, nonlinear_distance(x))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 38,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "noprint = lambda *args, **kwargs: None"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 74,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def run_ser_test(sample_duration=128, nbits=6, signal_amplitude=2.0e-3, decimation=10, threshold_factor=4.0, power_avg_width=2.5, max_lookahead=6.5, seed=0, ax=None, print=print):\n",
+    "\n",
+    "    test_data, signal = generate_test_signal(sample_duration, nbits, signal_amplitude, decimation, seed)\n",
+    "    cor_an = correlate(signal, nbits=nbits, decimation=decimation)\n",
+    "\n",
+    "    power_avg_width = int(power_avg_width * (2**nbits - 1) * decimation)\n",
+    "\n",
+    "    bit_period = (2**nbits) * decimation\n",
+    "    peak_group_threshold = 0.1 * bit_period\n",
+    "    \n",
+    "    cwt_res = np.array([ sig.cwt(row, sig.ricker, [0.73 * decimation]).flatten() for row in cor_an ])\n",
+    "    if ax:\n",
+    "        ax.grid()\n",
+    "        ax.plot(cwt_res.T)\n",
+    "        \n",
+    "    th = np.array([ np.convolve(np.abs(row), np.ones((power_avg_width,))/power_avg_width, mode='same') for row in cwt_res ])\n",
+    "\n",
+    "    def compare_th(elem):\n",
+    "        idx, (th, val) = elem\n",
+    "        #print('compare_th:', th.shape, val.shape)\n",
+    "        return np.any(np.abs(val) > th*threshold_factor)\n",
+    "\n",
+    "    peaks = [ list(group) for val, group in itertools.groupby(enumerate(zip(th.T, cwt_res.T)), compare_th) if val ]\n",
+    "    peak_group = []\n",
+    "    for group in peaks:\n",
+    "        pos = np.mean([idx for idx, _val in group])\n",
+    "        pol = np.mean([max(val.min(), val.max(), key=abs) for _idx, (_th, val) in group])\n",
+    "        pol_idx = np.argmax(np.bincount([ np.argmax(np.abs(val)) for _idx, (_th, val) in group ]))\n",
+    "        #print(f'group', pos, pol, pol_idx)\n",
+    "        #for pol, (_idx, (_th, val)) in zip([max(val.min(), val.max(), key=abs) for _idx, (_th, val) in group], group):\n",
+    "        #    print('    ', pol, val)\n",
+    "        if ax:\n",
+    "            ax.axvline(pos, color='cyan', alpha=0.3)\n",
+    "\n",
+    "        if not peak_group or pos - peak_group[-1][1] > peak_group_threshold:\n",
+    "            if peak_group:\n",
+    "                peak_pos = peak_group[-1][3]\n",
+    "                if ax:\n",
+    "                    ax.axvline(peak_pos, color='red', alpha=0.6)\n",
+    "                #ax3.text(peak_pos-20, 2.0, f'{0 if pol < 0 else 1}', horizontalalignment='right', verticalalignment='center', color='black')\n",
+    "\n",
+    "            peak_group.append((pos, pos, pol, pos, pol_idx))\n",
+    "            #ax3.axvline(pos, color='cyan', alpha=0.5)\n",
+    "\n",
+    "        else:\n",
+    "            group_start, last_pos, last_pol, peak_pos, last_pol_idx = peak_group[-1]\n",
+    "\n",
+    "            if abs(pol) > abs(last_pol):\n",
+    "                #ax3.axvline(pos, color='magenta', alpha=0.5)\n",
+    "                peak_group[-1] = (group_start, pos, pol, pos, pol_idx)\n",
+    "            else:\n",
+    "                #ax3.axvline(pos, color='blue', alpha=0.5)\n",
+    "                peak_group[-1] = (group_start, pos, last_pol, peak_pos, last_pol_idx)\n",
+    "\n",
+    "    avg_peak = np.mean(np.abs(np.array([last_pol for _1, _2, last_pol, _3, _4 in peak_group])))\n",
+    "    print('avg_peak', avg_peak)\n",
+    "\n",
+    "    noprint = lambda *args, **kwargs: None\n",
+    "    def mle_decode(peak_groups, print=print):\n",
+    "        peak_groups = [ (pos, pol, idx) for _1, _2, pol, pos, idx in peak_groups ]\n",
+    "        candidates = [ (0, [(pos, pol, idx)]) for pos, pol, idx in peak_groups if pos < bit_period*2.5 ]\n",
+    "\n",
+    "        while candidates:\n",
+    "            chain_candidates = []\n",
+    "            for chain_score, chain in candidates:\n",
+    "                pos, ampl, _idx = chain[-1]\n",
+    "                score_fun = lambda pos, npos, npol: abs(npol)/avg_peak + nonlinear_distance((npos-pos)/bit_period)\n",
+    "                next_candidates = sorted([ (score_fun(pos, npos, npol), npos, npol, nidx) for npos, npol, nidx in peak_groups if pos < npos < pos + bit_period*max_lookahead ], reverse=True)\n",
+    "\n",
+    "                print(f'    candidates for {pos}, {ampl}:')\n",
+    "                for score, npos, npol, nidx in next_candidates:\n",
+    "                    print(f'        {score:.4f} {npos:.2f} {npol:.2f} {nidx:.2f}')\n",
+    "\n",
+    "                nch, cor_len = cor_an.shape\n",
+    "                if cor_len - pos < 1.5*bit_period or not next_candidates:\n",
+    "                    score = sum(score_fun(opos, npos, npol) for (opos, _opol, _oidx), (npos, npol, _nidx) in zip(chain[:-1], chain[1:])) / len(chain)\n",
+    "                    yield score, chain\n",
+    "\n",
+    "                else:\n",
+    "                    print('extending')\n",
+    "                    for score, npos, npol, nidx in next_candidates[:3]:\n",
+    "                        if score > 0.5:\n",
+    "                            new_chain_score = chain_score * 0.9 + score * 0.1\n",
+    "                            chain_candidates.append((new_chain_score, chain + [(npos, npol, nidx)]))\n",
+    "            print('chain candidates:')\n",
+    "            for score, chain in sorted(chain_candidates, reverse=True):\n",
+    "                print('    ', [(score, [(f'{pos:.2f}', f'{pol:.2f}') for pos, pol, _idx in chain])])\n",
+    "            candidates = [ (chain_score, chain) for chain_score, chain in sorted(chain_candidates, reverse=True)[:10] ]\n",
+    "\n",
+    "    res = sorted(mle_decode(peak_group, print=noprint), reverse=True)\n",
+    "    #for i, (score, chain) in enumerate(res):\n",
+    "    #    print(f'Chain {i}@{score:.4f}: {chain}')\n",
+    "    (_score, chain), *_ = res\n",
+    "\n",
+    "    def viz(chain):\n",
+    "        last_pos = None\n",
+    "        for pos, pol, nidx in chain:\n",
+    "            if last_pos:\n",
+    "                delta = int(round((pos - last_pos) / bit_period))\n",
+    "                if delta > 1:\n",
+    "                    print(f'skipped {delta} symbols at {pos}')\n",
+    "                for i in range(delta-1):\n",
+    "                    yield None\n",
+    "            decoded = nidx*2 + (0 if pol < 0 else 1)\n",
+    "            yield decoded\n",
+    "            if ax:\n",
+    "                ax.axvline(pos, color='blue', alpha=0.5)\n",
+    "                ax.text(pos-20, 0.0, f'{decoded}', horizontalalignment='right', verticalalignment='center', color='black')\n",
+    "\n",
+    "            last_pos = pos\n",
+    "\n",
+    "    decoded = list(viz(chain))\n",
+    "    print('decoding [ref|dec]:')\n",
+    "    failures = 0\n",
+    "    for i, (ref, found) in enumerate(itertools.zip_longest(test_data, decoded)):\n",
+    "        print(f'{ref or -1:>3d}|{found or -1:>3d} {\"✔\" if ref==found else \"✘\" if found else \" \"}', end='    ')\n",
+    "        if ref != found:\n",
+    "            failures += 1\n",
+    "        if i%8 == 7:\n",
+    "            print()\n",
+    "    ser = failures/len(test_data)\n",
+    "    print(f'Symbol error rate e={ser}')\n",
+    "    br = sampling_rate / decimation / (2**nbits) * nbits * (1 - ser) * 3600\n",
+    "    print(f'maximum bitrate r={br} b/h')\n",
+    "    return ser, br"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "9b7546a233fb4b6cb6e8f809250ba768",
+       "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": "stdout",
+     "output_type": "stream",
+     "text": [
+      "(63,) (63,)\n",
+      "(63,) (63,)\n",
+      "avg_peak 1.6845488102985742\n",
+      "skipped 2 symbols at 10079.0\n",
+      "skipped 2 symbols at 30870.0\n",
+      "skipped 2 symbols at 42209.5\n",
+      "decoding [ref|dec]:\n",
+      " 44| 44 ✔     47| 47 ✔    117|117 ✔     64| 64 ✔     67| 67 ✔    123|123 ✔     67| 67 ✔    103|103 ✔    \n",
+      "  9|  9 ✔     83| 83 ✔     21| 21 ✔    114|114 ✔     36| 36 ✔     87| 87 ✔     70| -1       88| 88 ✔    \n",
+      " 88| 88 ✔     12| 12 ✔     58| 58 ✔     65| 65 ✔    102|102 ✔     39| 39 ✔     87| 87 ✔     46| 46 ✔    \n",
+      " 88| 88 ✔     81| 81 ✔     37| 37 ✔     25| 25 ✔     77| 77 ✔     72| 72 ✔      9|  9 ✔     20| 20 ✔    \n",
+      "115|115 ✔     80| 80 ✔    115|115 ✔     69| 69 ✔    126|126 ✔     79| 79 ✔     47| 47 ✔     64| 64 ✔    \n",
+      " 82| 82 ✔     99| 99 ✔     88| 88 ✔     49| 49 ✔    115|115 ✔     29| 29 ✔     19| 19 ✔     19| -1      \n",
+      " 14| 14 ✔     39| 39 ✔     32| 32 ✔     65| 65 ✔      9|  9 ✔     57| 57 ✔    127|127 ✔     32| 32 ✔    \n",
+      " 31| 31 ✔     74| 74 ✔    116|116 ✔     23| 23 ✔     35| 35 ✔    126|126 ✔     75| 75 ✔    114|114 ✔    \n",
+      " 55| 55 ✔     28| -1       34| 34 ✔     -1| -1 ✔     -1| -1 ✔     36| 36 ✔     53| 53 ✔      5|  5 ✔    \n",
+      " 38| 38 ✔    104|104 ✔    116|116 ✔     17| 17 ✔     79| 79 ✔      4|  4 ✔    105|105 ✔     42| 42 ✔    \n",
+      " 58| 58 ✔     31| 31 ✔    120|120 ✔      1|  1 ✔     65| 65 ✔    103|103 ✔     41| 41 ✔     57| 57 ✔    \n",
+      " 35| 35 ✔    102|103 ✘    119|119 ✔     11| 11 ✔     46| 46 ✔     82| 82 ✔     91| 91 ✔     -1| -1 ✔    \n",
+      " 14| 14 ✔     99| 99 ✔     53| 53 ✔     12| 12 ✔    121|121 ✔     42| 42 ✔     84| 84 ✔     75| 75 ✔    \n",
+      " 68| 68 ✔      6|  6 ✔     68| 68 ✔     47| 47 ✔    127|127 ✔    116|116 ✔      3|  3 ✔     76| 76 ✔    \n",
+      "100|100 ✔     52| 52 ✔    104|104 ✔     78| 78 ✔     15| 15 ✔     20| 20 ✔     99| 99 ✔     58| 58 ✔    \n",
+      " 23| 23 ✔     79| 79 ✔     13| 13 ✔    117|117 ✔     85| 85 ✔     48| 48 ✔     49| 49 ✔     69| 69 ✔    \n",
+      "Symbol error rate e=0.03125\n",
+      "maximum bitrate r=326.953125 b/h\n"
+     ]
+    }
+   ],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(12,5))\n",
+    "run_ser_test(ax=ax)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 84,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "<ipython-input-84-f36c7b0ffb61>:13: 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": "b5eb2cd4f4224f75bc3dc73b6143d849",
+       "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": "stdout",
+     "output_type": "stream",
+     "text": [
+      "nbits=5\n",
+      "signal_amplitude=0.00029: ser=1.01000 ±0.012207515615390381, br=-5.62500\n",
+      "signal_amplitude=0.00020: ser=1.01469 ±0.013678792892649557, br=-8.26172\n",
+      "signal_amplitude=0.00052: ser=1.00406 ±0.018895270572288715, br=-2.28516\n",
+      "signal_amplitude=0.00043: ser=1.01000 ±0.018817586521655747, br=-5.62500\n",
+      "signal_amplitude=0.00024: ser=1.01156 ±0.014510233457804875, br=-6.50391\n",
+      "signal_amplitude=0.00036: ser=1.00687 ±0.014875787794264881, br=-3.86719\n",
+      "signal_amplitude=0.00032: ser=1.01156 ±0.01837117307087384, br=-6.50391\n",
+      "signal_amplitude=0.00022: ser=1.01000 ±0.013535254892317322, br=-5.62500\n",
+      "signal_amplitude=0.00057: ser=1.00281 ±0.012476540486048206, br=-1.58203\n",
+      "signal_amplitude=0.00039: ser=1.00812 ±0.014402148277253642, br=-4.57031\n",
+      "signal_amplitude=0.00047: ser=1.00438 ±0.012899854650343935, br=-2.46094\n",
+      "signal_amplitude=0.00027: ser=1.00937 ±0.014657549249448218, br=-5.27344\n",
+      "signal_amplitude=0.00063: ser=0.99938 ±0.019253652250936705, br=0.35156\n",
+      "signal_amplitude=0.00077: ser=0.99156 ±0.03231920868461974, br=4.74609\n",
+      "signal_amplitude=0.00093: ser=0.95156 ±0.06625442202223185, br=27.24609\n",
+      "signal_amplitude=0.00112: ser=0.76000 ±0.2099632594348354, br=135.00000\n",
+      "signal_amplitude=0.00136: ser=0.51375 ±0.30673813139223494, br=273.51562\n",
+      "signal_amplitude=0.00165: ser=0.39844 ±0.38814210912370745, br=338.37891\n",
+      "signal_amplitude=0.00070: ser=0.99281 ±0.023688242072809035, br=4.04297\n",
+      "signal_amplitude=0.00084: ser=0.96375 ±0.050769469787461836, br=20.39062\n",
+      "signal_amplitude=0.00102: ser=0.91063 ±0.10310321739645179, br=50.27344\n",
+      "signal_amplitude=0.00124: ser=0.72500 ±0.23567348639059932, br=154.68750\n",
+      "signal_amplitude=0.00150: ser=0.40969 ±0.3064419041596629, br=332.05078\n",
+      "signal_amplitude=0.00182: ser=0.32531 ±0.38085840544748384, br=379.51172\n",
+      "signal_amplitude=0.00200: ser=0.29000 ±0.3885339029608613, br=399.37500\n",
+      "nbits=6\n",
+      "signal_amplitude=0.00052: ser=1.00375 ±0.027432445434193427, br=-1.26562\n",
+      "signal_amplitude=0.00029: ser=1.01531 ±0.013528038013695853, br=-5.16797\n",
+      "signal_amplitude=0.00020: ser=1.02000 ±0.01698459780212649, br=-6.75000\n",
+      "signal_amplitude=0.00024: ser=1.01844 ±0.0197494066366562, br=-6.22266\n",
+      "signal_amplitude=0.00043: ser=1.01000 ±0.013535254892317322, br=-3.37500\n",
+      "signal_amplitude=0.00036: ser=1.01500 ±0.01860884366369926, br=-5.06250\n",
+      "signal_amplitude=0.00032: ser=1.00906 ±0.01443601182806387, br=-3.05859\n",
+      "signal_amplitude=0.00022: ser=1.01656 ±0.015200483133769137, br=-5.58984\n",
+      "signal_amplitude=0.00057: ser=0.98281 ±0.04926213365760764, br=5.80078\n",
+      "signal_amplitude=0.00027: ser=1.02000 ±0.015946688527716343, br=-6.75000\n",
+      "signal_amplitude=0.00047: ser=1.00687 ±0.02815276407388802, br=-2.32031\n",
+      "signal_amplitude=0.00039: ser=1.00906 ±0.016189792308735775, br=-3.05859\n",
+      "signal_amplitude=0.00077: ser=0.76906 ±0.23454244018940368, br=77.94141\n",
+      "signal_amplitude=0.00063: ser=0.94031 ±0.08557822627572974, br=20.14453\n",
+      "signal_amplitude=0.00112: ser=0.29750 ±0.347296478171029, br=237.09375\n",
+      "signal_amplitude=0.00093: ser=0.50125 ±0.3293776683952632, br=168.32812\n",
+      "signal_amplitude=0.00136: ser=0.37250 ±0.42536588111001566, br=211.78125\n",
+      "signal_amplitude=0.00165: ser=0.51000 ±0.46215950303980546, br=165.37500\n",
+      "signal_amplitude=0.00070: ser=0.90063 ±0.1848975645458858, br=33.53906\n",
+      "signal_amplitude=0.00084: ser=0.64687 ±0.26652421325275494, br=119.17969\n",
+      "signal_amplitude=0.00124: ser=0.38500 ±0.39889079606767064, br=207.56250\n",
+      "signal_amplitude=0.00102: ser=0.40875 ±0.3467111099315971, br=199.54688\n",
+      "signal_amplitude=0.00150: ser=0.40375 ±0.4435118198819508, br=201.23438\n",
+      "signal_amplitude=0.00182: ser=0.58531 ±0.46179168734397985, br=139.95703\n",
+      "signal_amplitude=0.00200: ser=0.61594 ±0.4584529436730666, br=129.62109\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.legend.Legend at 0x7fe2bdd64430>"
+      ]
+     },
+     "execution_count": 84,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "sample_duration=128\n",
+    "sample_reps = 25\n",
+    "sweep_points = 25\n",
+    "\n",
+    "default_params = dict(\n",
+    "        nbits=6,\n",
+    "        signal_amplitude=2.0e-3,\n",
+    "        decimation=10,\n",
+    "        threshold_factor=4.0,\n",
+    "        power_avg_width=2.5,\n",
+    "        max_lookahead=6.5)\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(12, 9))\n",
+    "\n",
+    "for nbits in [5, 6]: # FIXME make sim stable for higher bit counts\n",
+    "    print(f'nbits={nbits}')\n",
+    "    \n",
+    "    def calculate_ser(v):\n",
+    "        params = dict(default_params)\n",
+    "        params['signal_amplitude'] = v\n",
+    "        params['nbits'] = nbits\n",
+    "        sers, brs = [], []\n",
+    "        for i in range(sample_reps):\n",
+    "            ser, br = run_ser_test(**params, sample_duration=sample_duration, print=noprint, seed=np.random.randint(0xffffffff))\n",
+    "            sers.append(ser)\n",
+    "            brs.append(br)\n",
+    "        sers, brs = np.array(sers), np.array(brs)\n",
+    "        ser = np.mean(sers)\n",
+    "        print(f'signal_amplitude={v:<.5f}: ser={ser:<.5f} ±{np.std(sers):<.5f}, br={np.mean(brs):<.5f}')\n",
+    "        return ser\n",
+    "    \n",
+    "    vs = 0.2e-3 * 10 ** np.linspace(0, 1.0, sweep_points)\n",
+    "    with Pool(6) as p:\n",
+    "        data = p.map(calculate_ser, vs)\n",
+    "    \n",
+    "    ax.plot(vs, data, label=f'{nbits} bit')\n",
+    "ax.grid()\n",
+    "ax.set_xlabel('Amplitude in mHz')\n",
+    "ax.set_ylabel('Symbol error rate')\n",
+    "ax.legend()"
+   ]
+  }
+ ],
+ "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/dsss_experiments.ipynb b/lab-windows/dsss_experiments.ipynb
index be41b59..f2bbc0b 100644
--- a/lab-windows/dsss_experiments.ipynb
+++ b/lab-windows/dsss_experiments.ipynb
@@ -30,6 +30,15 @@
   },
   {
    "cell_type": "code",
+   "execution_count": 105,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sampling_rate = 10 # sp/s"
+   ]
+  },
+  {
+   "cell_type": "code",
    "execution_count": 3,
    "metadata": {},
    "outputs": [],
@@ -410,29 +419,29 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 96,
+   "execution_count": 145,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "(31,) (31,)\n",
-      "(31,) (31,)\n"
+      "(63,) (63,)\n",
+      "(63,) (63,)\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "<ipython-input-96-b3aae757ccad>:33: 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-145-babcf8a4e867>:33: 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, ((ax1, ax3), (ax2, ax4)) = plt.subplots(2, 2, figsize=(16, 9))\n"
      ]
     },
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "246c19a3c1424e7eb0b675ff060ea5b3",
+       "model_id": "10aa67d294304f2ba26c8e6d5555d5e6",
        "version_major": 2,
        "version_minor": 0
       },
@@ -446,10 +455,10 @@
     {
      "data": {
       "text/plain": [
-       "(0.002, 0.013899708)"
+       "(0.002, 0.014074279)"
       ]
      },
-     "execution_count": 96,
+     "execution_count": 145,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -457,15 +466,15 @@
    "source": [
     "decimation = 10\n",
     "signal_amplitude = 2.0e-3\n",
-    "nbits = 5\n",
+    "nbits = 6\n",
     "\n",
     "#test_data = np.random.randint(0, 2, 100)\n",
     "#test_data = np.array([0, 1, 0, 0, 1, 1, 1, 0])\n",
-    "#test_data = np.random.RandomState(seed=0).randint(0, 2 * (2**nbits), 64)\n",
+    "test_data = np.random.RandomState(seed=0).randint(0, 2 * (2**nbits), 64)\n",
     "#test_data = np.random.RandomState(seed=0).randint(0, 8, 64)\n",
     "#test_data = np.array(list(range(8)) * 8)\n",
     "#test_data = np.array([0, 1] * 32)\n",
-    "test_data = np.array(list(range(64)))\n",
+    "#test_data = np.array(list(range(64)))\n",
     "\n",
     "foo = np.repeat(modulate(test_data, nbits) * 2.0 - 1, decimation) * signal_amplitude\n",
     "noise = np.resize(mains_noise, len(foo))\n",
@@ -493,12 +502,20 @@
     "ax1.plot(foo + noise)\n",
     "ax1.plot(foo)\n",
     "ax1.set_title('raw')\n",
+    "ax1.grid(axis='y')\n",
     "\n",
     "ax2.plot(filtered)\n",
     "ax2.plot(foo)\n",
     "ax2.set_title('filtered')\n",
+    "ax2.grid(axis='y')\n",
+    "\n",
+    "for i in range(0, len(foo) + 1, decimation*(2**nbits - 1)):\n",
+    "    ax1.axvline(i, color='gray', alpha=0.5, lw=1)\n",
+    "    ax2.axvline(i, color='gray', alpha=0.5, lw=1)\n",
     "\n",
-    "ax3.plot(cor1.T)\n",
+    "for i, (color, trace) in enumerate(zip(plt.cm.winter(np.linspace(0, 1, cor1.shape[0])), cor1.T)):\n",
+    "    if i%3 == 0:\n",
+    "        ax3.plot(trace + 0.5 * i, alpha=1.0, color=color)\n",
     "ax3.set_title('corr raw')\n",
     "ax3.grid()\n",
     "\n",
@@ -678,21 +695,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 97,
+   "execution_count": 146,
    "metadata": {},
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "<ipython-input-97-2d2c2f814215>:11: 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-146-badd40342f73>:11: 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, (ax1, ax3) = plt.subplots(2, figsize=(12, 5))\n"
      ]
     },
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "b083c661b5b441d6b7fc45201faa0576",
+       "model_id": "509bf67d93b74741b48bca58529c4b9d",
        "version_major": 2,
        "version_minor": 0
       },
@@ -707,26 +724,24 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "cor_an (33, 20149)\n",
-      "cwt_res (33, 20149)\n",
-      "th (33, 20149)\n",
-      "[((33,), (33,)), ((33,), (33,)), ((33,), (33,)), ((33,), (33,)), ((33,), (33,))]\n",
-      "peaks: 180\n",
-      "avg_peak 1.058897833824206\n",
-      "skipped 3 symbols at 5889.0\n",
-      "skipped 2 symbols at 8369.0\n",
-      "skipped 2 symbols at 14568.5\n",
-      "skipped 2 symbols at 16739.0\n",
+      "cor_an (65, 40949)\n",
+      "cwt_res (65, 40949)\n",
+      "th (65, 40949)\n",
+      "[((65,), (65,)), ((65,), (65,)), ((65,), (65,)), ((65,), (65,)), ((65,), (65,))]\n",
+      "peaks: 982\n",
+      "avg_peak 1.6673786030736735\n",
+      "skipped 2 symbols at 30238.5\n",
       "decoding [ref|dec]:\n",
-      " -1| -1 ✔      1|  1 ✔      2|  2 ✔      3|  3 ✔      4|  4 ✔      5|  5 ✔      6|  6 ✔      7|  7 ✔    \n",
-      "  8|  8 ✔      9|  9 ✔     10| 10 ✔     11| 11 ✔     12| 12 ✔     13| 13 ✔     14| 14 ✔     15| 15 ✔    \n",
-      " 16| -1 ✘     17| -1 ✘     18| 18 ✔     19| 19 ✔     20| 20 ✔     21| 21 ✔     22| 22 ✔     23| 23 ✔    \n",
-      " 24| 24 ✔     25| -1 ✘     26| 26 ✔     27| 27 ✔     28| 28 ✔     29| 29 ✔     30| 30 ✔     31| 31 ✔    \n",
-      " 32| 32 ✔     33| 33 ✔     34| 34 ✔     35| 35 ✔     36| 36 ✔     37| 37 ✔     38| 38 ✔     39| 39 ✔    \n",
-      " 40| 40 ✔     41| 41 ✔     42| 42 ✔     43| 43 ✔     44| 44 ✔     45| -1 ✘     46| 46 ✔     47| 47 ✔    \n",
-      " 48| 48 ✔     49| 49 ✔     50| 50 ✔     51| 51 ✔     52| -1 ✘     53| 53 ✔     54| 54 ✔     55| 55 ✔    \n",
-      " 56| 56 ✔     57| 57 ✔     58| 58 ✔     59| 59 ✔     60| 60 ✔     61| 61 ✔     62| 62 ✔     63| 56 ✘    \n",
-      "Symbol error rate r=0.09375\n"
+      " 44| 44 ✔     47| 47 ✔    117|117 ✔     64| 64 ✔     67| 67 ✔    123|123 ✔     67| 67 ✔    103|103 ✔    \n",
+      "  9|  9 ✔     83| 83 ✔     21| 21 ✔    114|114 ✔     36| 36 ✔     87| 87 ✔     70| 70 ✔     88| 88 ✔    \n",
+      " 88| 88 ✔     12| 12 ✔     58| 58 ✔     65| 65 ✔    102|102 ✔     39| 39 ✔     87| 87 ✔     46| 46 ✔    \n",
+      " 88| 88 ✔     81| 81 ✔     37| 37 ✔     25| 25 ✔     77| 77 ✔     72| 72 ✔      9|  9 ✔     20| 20 ✔    \n",
+      "115|115 ✔     80| 80 ✔    115|115 ✔     69| 69 ✔    126|126 ✔     79| 79 ✔     47| 47 ✔     64| 64 ✔    \n",
+      " 82| 82 ✔     99| 99 ✔     88| 88 ✔     49| 49 ✔    115|115 ✔     29| 29 ✔     19| -1       19| 19 ✔    \n",
+      " 14| 14 ✔     39| 39 ✔     32| 32 ✔     65| 64 ✘      9|  9 ✔     57| 57 ✔    127|127 ✔     32| 32 ✔    \n",
+      " 31| 31 ✔     74| 74 ✔    116|116 ✔     23| 23 ✔     35| 35 ✔    126|126 ✔     75| 75 ✔    114| 26 ✘    \n",
+      "Symbol error rate e=0.046875\n",
+      "maximum bitrate r=321.6796875 b/h\n"
      ]
     }
    ],
@@ -869,12 +884,13 @@
     "print('decoding [ref|dec]:')\n",
     "failures = 0\n",
     "for i, (ref, found) in enumerate(itertools.zip_longest(test_data, decoded)):\n",
-    "    print(f'{ref or -1:>3d}|{found or -1:>3d} {\"✔\" if ref==found else \"✘\"}', end='    ')\n",
+    "    print(f'{ref or -1:>3d}|{found or -1:>3d} {\"✔\" if ref==found else \"✘\" if found else \" \"}', end='    ')\n",
     "    if ref != found:\n",
     "        failures += 1\n",
     "    if i%8 == 7:\n",
     "        print()\n",
-    "print(f'Symbol error rate r={failures/len(test_data)}')\n",
+    "print(f'Symbol error rate e={failures/len(test_data)}')\n",
+    "print(f'maximum bitrate r={sampling_rate / decimation / (2**nbits) * nbits * (1 - failures/len(test_data)) * 3600} b/h')\n",
     "#ax3.plot(th)"
    ]
   },