{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 92,
   "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",
    "from collections import defaultdict\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": 100,
   "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, ser_maxshift=3):\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 = defaultdict(lambda: 0)\n",
    "    for shift in range(-ser_maxshift, ser_maxshift):\n",
    "        print(f'=== shift = {shift} ===')\n",
    "        a = test_data if shift > 0 else test_data[-shift:]\n",
    "        b = decoded if shift < 0 else decoded[shift:]\n",
    "        for i, (ref, found) in enumerate(itertools.zip_longest(a, b)):\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[shift] += 1\n",
    "            if i%8 == 7:\n",
    "                print()\n",
    "        if i%8 != 7:\n",
    "            print()\n",
    "    ser = min(failures.values())/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": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-111-2589572e37aa>: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": "794478ea08254c9da3608685434de9d6",
       "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=0.99141 ±0.01434, br=4.83398\n",
      "signal_amplitude=0.00020: ser=0.99344 ±0.01365, br=3.69141\n",
      "signal_amplitude=0.00052: ser=0.98953 ±0.01621, br=5.88867\n",
      "signal_amplitude=0.00024: ser=0.98984 ±0.01426, br=5.71289\n",
      "signal_amplitude=0.00043: ser=0.98906 ±0.01531, br=6.15234\n",
      "signal_amplitude=0.00036: ser=0.98859 ±0.01451, br=6.41602\n",
      "signal_amplitude=0.00032: ser=0.98531 ±0.01621, br=8.26172\n",
      "signal_amplitude=0.00022: ser=0.99125 ±0.01377, br=4.92188\n",
      "signal_amplitude=0.00027: ser=0.99172 ±0.01339, br=4.65820\n",
      "signal_amplitude=0.00047: ser=0.98938 ±0.01821, br=5.97656\n",
      "signal_amplitude=0.00057: ser=0.98453 ±0.02223, br=8.70117\n",
      "signal_amplitude=0.00039: ser=0.99328 ±0.01499, br=3.77930\n",
      "signal_amplitude=0.00063: ser=0.97531 ±0.02277, br=13.88672\n",
      "signal_amplitude=0.00077: ser=0.95766 ±0.03962, br=23.81836\n",
      "signal_amplitude=0.00093: ser=0.89844 ±0.07690, br=57.12891\n",
      "signal_amplitude=0.00112: ser=0.69250 ±0.15978, br=172.96875\n",
      "signal_amplitude=0.00136: ser=0.34797 ±0.13550, br=366.76758\n",
      "signal_amplitude=0.00165: ser=0.15875 ±0.07678, br=473.20312\n",
      "signal_amplitude=0.00070: ser=0.97562 ±0.02444, br=13.71094\n",
      "signal_amplitude=0.00084: ser=0.93437 ±0.05440, br=36.91406\n",
      "signal_amplitude=0.00102: ser=0.83734 ±0.10122, br=91.49414\n",
      "signal_amplitude=0.00124: ser=0.56047 ±0.18288, br=247.23633\n",
      "signal_amplitude=0.00150: ser=0.24266 ±0.09904, br=426.00586\n",
      "signal_amplitude=0.00182: ser=0.10359 ±0.03178, br=504.22852\n",
      "signal_amplitude=0.00200: ser=0.06453 ±0.03562, br=526.20117\n",
      "nbits=6\n",
      "signal_amplitude=0.00024: ser=1.00156 ±0.01855, br=-0.52734\n",
      "signal_amplitude=0.00020: ser=1.00172 ±0.01365, br=-0.58008\n",
      "signal_amplitude=0.00052: ser=0.95688 ±0.04998, br=14.55469\n",
      "signal_amplitude=0.00029: ser=1.00078 ±0.01391, br=-0.26367\n",
      "signal_amplitude=0.00043: ser=0.98406 ±0.02215, br=5.37891\n",
      "signal_amplitude=0.00036: ser=0.99844 ±0.01733, br=0.52734\n",
      "signal_amplitude=0.00027: ser=0.99469 ±0.01339, br=1.79297\n",
      "signal_amplitude=0.00032: ser=0.99766 ±0.01541, br=0.79102\n",
      "signal_amplitude=0.00047: ser=0.97813 ±0.03233, br=7.38281\n",
      "signal_amplitude=0.00022: ser=1.00016 ±0.01612, br=-0.05273\n",
      "signal_amplitude=0.00057: ser=0.93609 ±0.07689, br=21.56836\n",
      "signal_amplitude=0.00039: ser=0.99625 ±0.01883, br=1.26562\n",
      "signal_amplitude=0.00063: ser=0.86875 ±0.08992, br=44.29688\n",
      "signal_amplitude=0.00077: ser=0.58109 ±0.17566, br=141.38086\n",
      "signal_amplitude=0.00093: ser=0.28594 ±0.13348, br=240.99609\n",
      "signal_amplitude=0.00112: ser=0.12547 ±0.06352, br=295.15430\n",
      "signal_amplitude=0.00165: ser=0.06500 ±0.01642, br=315.56250\n",
      "signal_amplitude=0.00136: ser=0.08516 ±0.06513, br=308.75977\n",
      "signal_amplitude=0.00070: ser=0.75609 ±0.15173, br=82.31836\n",
      "signal_amplitude=0.00084: ser=0.46125 ±0.18383, br=181.82812\n",
      "signal_amplitude=0.00102: ser=0.19703 ±0.10488, br=271.00195\n",
      "signal_amplitude=0.00124: ser=0.10562 ±0.06833, br=301.85156\n",
      "signal_amplitude=0.00182: ser=0.06563 ±0.02007, br=315.35156\n",
      "signal_amplitude=0.00150: ser=0.06984 ±0.01818, br=313.92773\n",
      "signal_amplitude=0.00200: ser=0.06500 ±0.02196, br=315.56250\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fe3073ed0a0>"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_duration=128\n",
    "sample_reps = 50\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",
    "            seed = np.random.randint(0xffffffff)\n",
    "            ser, br = run_ser_test(**params, sample_duration=sample_duration, print=noprint, seed=seed)\n",
    "            sers.append(ser)\n",
    "            brs.append(br)\n",
    "            #print(f'nbits={nbits} ampl={v:>.5f} seed={seed:08x} > ser={ser:.5f}')\n",
    "        sers, brs = np.array(sers), np.array(brs)\n",
    "        ser, std = np.mean(sers), np.std(sers)\n",
    "        print(f'signal_amplitude={v:<.5f}: ser={ser:<.5f} ±{std:<.5f}, br={np.mean(brs):<.5f}')\n",
    "        return ser, std\n",
    "    \n",
    "    vs = 0.2e-3 * 10 ** np.linspace(0, 1.0, sweep_points)\n",
    "    with Pool(6) as p:\n",
    "        data = np.array(p.map(calculate_ser, vs))\n",
    "        sers, stds = data[:,0], data[:,1]\n",
    "    \n",
    "    l, = ax.plot(vs, np.clip(sers, 0, 1), label=f'{nbits} bit')\n",
    "    ax.fill_between(vs, np.clip(sers + stds, 0, 1), np.clip(sers - stds, 0, 1), facecolor=l.get_color(), alpha=0.3)\n",
    "ax.grid()\n",
    "ax.set_xlabel('Amplitude in mHz')\n",
    "ax.set_ylabel('Symbol error rate')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-101-d477d71c27f9>: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": "cfaf2edf2f214ca6a3e97994a14a80f5",
       "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": [
      "avg_peak 1.6752232386509318\n",
      "skipped 2 symbols at 15749.0\n",
      "skipped 3 symbols at 42209.0\n",
      "skipped 2 symbols at 57959.5\n",
      "skipped 2 symbols at 61108.5\n",
      "skipped 2 symbols at 77489.5\n",
      "decoding [ref|dec]:\n",
      "=== shift = -3 ===\n",
      "102| 69 ✘      2|124 ✘      3|102 ✘     78|  2 ✘     29|  3 ✘    122| 78 ✘     73| 29 ✘     98|123 ✘    \n",
      " 34| 73 ✘     -1| 98 ✘     97| 34 ✘      7| -1       97| 97 ✔     86|  7 ✘    120| 97 ✘     95| 86 ✘    \n",
      " 90|120 ✘     49| 95 ✘     89| 90 ✘     83| 49 ✘     19| 89 ✘     84| 83 ✘    117| -1       92| 84 ✘    \n",
      "119|117 ✘     16| 92 ✘     45|119 ✘     23| 16 ✘     16| 45 ✘    111| 23 ✘      9| 16 ✘     89|111 ✘    \n",
      " 18|  9 ✘     36| 89 ✘      2| 18 ✘    115| 36 ✘     40|  2 ✘    100|115 ✘    105| 40 ✘     93|100 ✘    \n",
      " 85|105 ✘    107| 93 ✘     90| 85 ✘     62|107 ✘    116| 90 ✘     42| 62 ✘    123|116 ✘     40| 42 ✘    \n",
      " -1|123 ✘     77| 40 ✘     40| -1       57| 77 ✘    110| 40 ✘     29| 57 ✘     94|110 ✘      1| 29 ✘    \n",
      " 29| 94 ✘     71|  1 ✘    119| 29 ✘     15| 71 ✘    115|119 ✘    120| 15 ✘     70|115 ✘     50| -1      \n",
      " 71| -1       50| 50 ✔     61| 71 ✘     38| 50 ✘      4| 61 ✘      3| 38 ✘    124|  4 ✘     95|  3 ✘    \n",
      " 27|124 ✘     48| 95 ✘    116| 27 ✘      3| 64 ✘     63|116 ✘     19|  3 ✘     79| 63 ✘      2| 19 ✘    \n",
      " 43| 79 ✘     92|  2 ✘      8| 43 ✘     65| 92 ✘     35|  8 ✘     30| 65 ✘     73| 35 ✘     73| 30 ✘    \n",
      " 38| 73 ✘     58| -1       49| 38 ✘     45| 58 ✘     58| 49 ✘     46| 45 ✘    116| -1      101| 46 ✘    \n",
      "  5|116 ✘     78|101 ✘    126|  5 ✘    105| 78 ✘    108|126 ✘     59| 76 ✘     46|108 ✘     27| 59 ✘    \n",
      " 14| 46 ✘     57| 27 ✘     81| 14 ✘      3| 57 ✘      9| 81 ✘    126|  3 ✘     18|  9 ✘     76|126 ✘    \n",
      "101| 55 ✘    124| 76 ✘      4|101 ✘      3|124 ✘    102|  4 ✘     79|  3 ✘    121|102 ✘    103| 79 ✘    \n",
      " 92| -1       30|103 ✘      4| 92 ✘    103| 30 ✘     59|  4 ✘     -1|103 ✘     -1| 48 ✘    \n",
      "=== shift = -2 ===\n",
      "124| 69 ✘    102|124 ✘      2|102 ✘      3|  2 ✘     78|  3 ✘     29| 78 ✘    122| 29 ✘     73|123 ✘    \n",
      " 98| 73 ✘     34| 98 ✘     -1| 34 ✘     97| -1        7| 97 ✘     97|  7 ✘     86| 97 ✘    120| 86 ✘    \n",
      " 95|120 ✘     90| 95 ✘     49| 90 ✘     89| 49 ✘     83| 89 ✘     19| 83 ✘     84| -1      117| 84 ✘    \n",
      " 92|117 ✘    119| 92 ✘     16|119 ✘     45| 16 ✘     23| 45 ✘     16| 23 ✘    111| 16 ✘      9|111 ✘    \n",
      " 89|  9 ✘     18| 89 ✘     36| 18 ✘      2| 36 ✘    115|  2 ✘     40|115 ✘    100| 40 ✘    105|100 ✘    \n",
      " 93|105 ✘     85| 93 ✘    107| 85 ✘     90|107 ✘     62| 90 ✘    116| 62 ✘     42|116 ✘    123| 42 ✘    \n",
      " 40|123 ✘     -1| 40 ✘     77| -1       40| 77 ✘     57| 40 ✘    110| 57 ✘     29|110 ✘     94| 29 ✘    \n",
      "  1| 94 ✘     29|  1 ✘     71| 29 ✘    119| 71 ✘     15|119 ✘    115| 15 ✘    120|115 ✘     70| -1      \n",
      " 50| -1       71| 50 ✘     50| 71 ✘     61| 50 ✘     38| 61 ✘      4| 38 ✘      3|  4 ✘    124|  3 ✘    \n",
      " 95|124 ✘     27| 95 ✘     48| 27 ✘    116| 64 ✘      3|116 ✘     63|  3 ✘     19| 63 ✘     79| 19 ✘    \n",
      "  2| 79 ✘     43|  2 ✘     92| 43 ✘      8| 92 ✘     65|  8 ✘     35| 65 ✘     30| 35 ✘     73| 30 ✘    \n",
      " 73| 73 ✔     38| -1       58| 38 ✘     49| 58 ✘     45| 49 ✘     58| 45 ✘     46| -1      116| 46 ✘    \n",
      "101|116 ✘      5|101 ✘     78|  5 ✘    126| 78 ✘    105|126 ✘    108| 76 ✘     59|108 ✘     46| 59 ✘    \n",
      " 27| 46 ✘     14| 27 ✘     57| 14 ✘     81| 57 ✘      3| 81 ✘      9|  3 ✘    126|  9 ✘     18|126 ✘    \n",
      " 76| 55 ✘    101| 76 ✘    124|101 ✘      4|124 ✘      3|  4 ✘    102|  3 ✘     79|102 ✘    121| 79 ✘    \n",
      "103| -1       92|103 ✘     30| 92 ✘      4| 30 ✘    103|  4 ✘     59|103 ✘     -1| 48 ✘    \n",
      "=== shift = -1 ===\n",
      " 69| 69 ✔    124|124 ✔    102|102 ✔      2|  2 ✔      3|  3 ✔     78| 78 ✔     29| 29 ✔    122|123 ✘    \n",
      " 73| 73 ✔     98| 98 ✔     34| 34 ✔     -1| -1 ✔     97| 97 ✔      7|  7 ✔     97| 97 ✔     86| 86 ✔    \n",
      "120|120 ✔     95| 95 ✔     90| 90 ✔     49| 49 ✔     89| 89 ✔     83| 83 ✔     19| -1       84| 84 ✔    \n",
      "117|117 ✔     92| 92 ✔    119|119 ✔     16| 16 ✔     45| 45 ✔     23| 23 ✔     16| 16 ✔    111|111 ✔    \n",
      "  9|  9 ✔     89| 89 ✔     18| 18 ✔     36| 36 ✔      2|  2 ✔    115|115 ✔     40| 40 ✔    100|100 ✔    \n",
      "105|105 ✔     93| 93 ✔     85| 85 ✔    107|107 ✔     90| 90 ✔     62| 62 ✔    116|116 ✔     42| 42 ✔    \n",
      "123|123 ✔     40| 40 ✔     -1| -1 ✔     77| 77 ✔     40| 40 ✔     57| 57 ✔    110|110 ✔     29| 29 ✔    \n",
      " 94| 94 ✔      1|  1 ✔     29| 29 ✔     71| 71 ✔    119|119 ✔     15| 15 ✔    115|115 ✔    120| -1      \n",
      " 70| -1       50| 50 ✔     71| 71 ✔     50| 50 ✔     61| 61 ✔     38| 38 ✔      4|  4 ✔      3|  3 ✔    \n",
      "124|124 ✔     95| 95 ✔     27| 27 ✔     48| 64 ✘    116|116 ✔      3|  3 ✔     63| 63 ✔     19| 19 ✔    \n",
      " 79| 79 ✔      2|  2 ✔     43| 43 ✔     92| 92 ✔      8|  8 ✔     65| 65 ✔     35| 35 ✔     30| 30 ✔    \n",
      " 73| 73 ✔     73| -1       38| 38 ✔     58| 58 ✔     49| 49 ✔     45| 45 ✔     58| -1       46| 46 ✔    \n",
      "116|116 ✔    101|101 ✔      5|  5 ✔     78| 78 ✔    126|126 ✔    105| 76 ✘    108|108 ✔     59| 59 ✔    \n",
      " 46| 46 ✔     27| 27 ✔     14| 14 ✔     57| 57 ✔     81| 81 ✔      3|  3 ✔      9|  9 ✔    126|126 ✔    \n",
      " 18| 55 ✘     76| 76 ✔    101|101 ✔    124|124 ✔      4|  4 ✔      3|  3 ✔    102|102 ✔     79| 79 ✔    \n",
      "121| -1      103|103 ✔     92| 92 ✔     30| 30 ✔      4|  4 ✔    103|103 ✔     59| 48 ✘    \n",
      "=== shift = 0 ===\n",
      " 10| 69 ✘     69|124 ✘    124|102 ✘    102|  2 ✘      2|  3 ✘      3| 78 ✘     78| 29 ✘     29|123 ✘    \n",
      "122| 73 ✘     73| 98 ✘     98| 34 ✘     34| -1       -1| 97 ✘     97|  7 ✘      7| 97 ✘     97| 86 ✘    \n",
      " 86|120 ✘    120| 95 ✘     95| 90 ✘     90| 49 ✘     49| 89 ✘     89| 83 ✘     83| -1       19| 84 ✘    \n",
      " 84|117 ✘    117| 92 ✘     92|119 ✘    119| 16 ✘     16| 45 ✘     45| 23 ✘     23| 16 ✘     16|111 ✘    \n",
      "111|  9 ✘      9| 89 ✘     89| 18 ✘     18| 36 ✘     36|  2 ✘      2|115 ✘    115| 40 ✘     40|100 ✘    \n",
      "100|105 ✘    105| 93 ✘     93| 85 ✘     85|107 ✘    107| 90 ✘     90| 62 ✘     62|116 ✘    116| 42 ✘    \n",
      " 42|123 ✘    123| 40 ✘     40| -1       -1| 77 ✘     77| 40 ✘     40| 57 ✘     57|110 ✘    110| 29 ✘    \n",
      " 29| 94 ✘     94|  1 ✘      1| 29 ✘     29| 71 ✘     71|119 ✘    119| 15 ✘     15|115 ✘    115| -1      \n",
      "120| -1       70| 50 ✘     50| 71 ✘     71| 50 ✘     50| 61 ✘     61| 38 ✘     38|  4 ✘      4|  3 ✘    \n",
      "  3|124 ✘    124| 95 ✘     95| 27 ✘     27| 64 ✘     48|116 ✘    116|  3 ✘      3| 63 ✘     63| 19 ✘    \n",
      " 19| 79 ✘     79|  2 ✘      2| 43 ✘     43| 92 ✘     92|  8 ✘      8| 65 ✘     65| 35 ✘     35| 30 ✘    \n",
      " 30| 73 ✘     73| -1       73| 38 ✘     38| 58 ✘     58| 49 ✘     49| 45 ✘     45| -1       58| 46 ✘    \n",
      " 46|116 ✘    116|101 ✘    101|  5 ✘      5| 78 ✘     78|126 ✘    126| 76 ✘    105|108 ✘    108| 59 ✘    \n",
      " 59| 46 ✘     46| 27 ✘     27| 14 ✘     14| 57 ✘     57| 81 ✘     81|  3 ✘      3|  9 ✘      9|126 ✘    \n",
      "126| 55 ✘     18| 76 ✘     76|101 ✘    101|124 ✘    124|  4 ✘      4|  3 ✘      3|102 ✘    102| 79 ✘    \n",
      " 79| -1      121|103 ✘    103| 92 ✘     92| 30 ✘     30|  4 ✘      4|103 ✘    103| 48 ✘     59| -1      \n",
      "=== shift = 1 ===\n",
      " 10|124 ✘     69|102 ✘    124|  2 ✘    102|  3 ✘      2| 78 ✘      3| 29 ✘     78|123 ✘     29| 73 ✘    \n",
      "122| 98 ✘     73| 34 ✘     98| -1       34| 97 ✘     -1|  7 ✘     97| 97 ✔      7| 86 ✘     97|120 ✘    \n",
      " 86| 95 ✘    120| 90 ✘     95| 49 ✘     90| 89 ✘     49| 83 ✘     89| -1       83| 84 ✘     19|117 ✘    \n",
      " 84| 92 ✘    117|119 ✘     92| 16 ✘    119| 45 ✘     16| 23 ✘     45| 16 ✘     23|111 ✘     16|  9 ✘    \n",
      "111| 89 ✘      9| 18 ✘     89| 36 ✘     18|  2 ✘     36|115 ✘      2| 40 ✘    115|100 ✘     40|105 ✘    \n",
      "100| 93 ✘    105| 85 ✘     93|107 ✘     85| 90 ✘    107| 62 ✘     90|116 ✘     62| 42 ✘    116|123 ✘    \n",
      " 42| 40 ✘    123| -1       40| 77 ✘     -1| 40 ✘     77| 57 ✘     40|110 ✘     57| 29 ✘    110| 94 ✘    \n",
      " 29|  1 ✘     94| 29 ✘      1| 71 ✘     29|119 ✘     71| 15 ✘    119|115 ✘     15| -1      115| -1      \n",
      "120| 50 ✘     70| 71 ✘     50| 50 ✔     71| 61 ✘     50| 38 ✘     61|  4 ✘     38|  3 ✘      4|124 ✘    \n",
      "  3| 95 ✘    124| 27 ✘     95| 64 ✘     27|116 ✘     48|  3 ✘    116| 63 ✘      3| 19 ✘     63| 79 ✘    \n",
      " 19|  2 ✘     79| 43 ✘      2| 92 ✘     43|  8 ✘     92| 65 ✘      8| 35 ✘     65| 30 ✘     35| 73 ✘    \n",
      " 30| -1       73| 38 ✘     73| 58 ✘     38| 49 ✘     58| 45 ✘     49| -1       45| 46 ✘     58|116 ✘    \n",
      " 46|101 ✘    116|  5 ✘    101| 78 ✘      5|126 ✘     78| 76 ✘    126|108 ✘    105| 59 ✘    108| 46 ✘    \n",
      " 59| 27 ✘     46| 14 ✘     27| 57 ✘     14| 81 ✘     57|  3 ✘     81|  9 ✘      3|126 ✘      9| 55 ✘    \n",
      "126| 76 ✘     18|101 ✘     76|124 ✘    101|  4 ✘    124|  3 ✘      4|102 ✘      3| 79 ✘    102| -1      \n",
      " 79|103 ✘    121| 92 ✘    103| 30 ✘     92|  4 ✘     30|103 ✘      4| 48 ✘    103| -1       59| -1      \n",
      "=== shift = 2 ===\n",
      " 10|102 ✘     69|  2 ✘    124|  3 ✘    102| 78 ✘      2| 29 ✘      3|123 ✘     78| 73 ✘     29| 98 ✘    \n",
      "122| 34 ✘     73| -1       98| 97 ✘     34|  7 ✘     -1| 97 ✘     97| 86 ✘      7|120 ✘     97| 95 ✘    \n",
      " 86| 90 ✘    120| 49 ✘     95| 89 ✘     90| 83 ✘     49| -1       89| 84 ✘     83|117 ✘     19| 92 ✘    \n",
      " 84|119 ✘    117| 16 ✘     92| 45 ✘    119| 23 ✘     16| 16 ✔     45|111 ✘     23|  9 ✘     16| 89 ✘    \n",
      "111| 18 ✘      9| 36 ✘     89|  2 ✘     18|115 ✘     36| 40 ✘      2|100 ✘    115|105 ✘     40| 93 ✘    \n",
      "100| 85 ✘    105|107 ✘     93| 90 ✘     85| 62 ✘    107|116 ✘     90| 42 ✘     62|123 ✘    116| 40 ✘    \n",
      " 42| -1      123| 77 ✘     40| 40 ✔     -1| 57 ✘     77|110 ✘     40| 29 ✘     57| 94 ✘    110|  1 ✘    \n",
      " 29| 29 ✔     94| 71 ✘      1|119 ✘     29| 15 ✘     71|115 ✘    119| -1       15| -1      115| 50 ✘    \n",
      "120| 71 ✘     70| 50 ✘     50| 61 ✘     71| 38 ✘     50|  4 ✘     61|  3 ✘     38|124 ✘      4| 95 ✘    \n",
      "  3| 27 ✘    124| 64 ✘     95|116 ✘     27|  3 ✘     48| 63 ✘    116| 19 ✘      3| 79 ✘     63|  2 ✘    \n",
      " 19| 43 ✘     79| 92 ✘      2|  8 ✘     43| 65 ✘     92| 35 ✘      8| 30 ✘     65| 73 ✘     35| -1      \n",
      " 30| 38 ✘     73| 58 ✘     73| 49 ✘     38| 45 ✘     58| -1       49| 46 ✘     45|116 ✘     58|101 ✘    \n",
      " 46|  5 ✘    116| 78 ✘    101|126 ✘      5| 76 ✘     78|108 ✘    126| 59 ✘    105| 46 ✘    108| 27 ✘    \n",
      " 59| 14 ✘     46| 57 ✘     27| 81 ✘     14|  3 ✘     57|  9 ✘     81|126 ✘      3| 55 ✘      9| 76 ✘    \n",
      "126|101 ✘     18|124 ✘     76|  4 ✘    101|  3 ✘    124|102 ✘      4| 79 ✘      3| -1      102|103 ✘    \n",
      " 79| 92 ✘    121| 30 ✘    103|  4 ✘     92|103 ✘     30| 48 ✘      4| -1      103| -1       59| -1      \n",
      "Symbol error rate e=0.0859375\n",
      "maximum bitrate r=308.49609375 b/h\n",
      "nbits=6 ampl=0.00387 seed=cbb3b8cf > ser=0.08594\n"
     ]
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 9))\n",
    "\n",
    "params = dict(default_params)\n",
    "params['signal_amplitude'] = 2.0e-3\n",
    "params['nbits'] = 6\n",
    "seed = 0xcbb3b8cf\n",
    "ser, br = run_ser_test(**params, sample_duration=sample_duration, print=print, seed=seed, ax=ax)\n",
    "print(f'nbits={nbits} ampl={v:>.5f} seed={seed:08x} > ser={ser:.5f}')"
   ]
  }
 ],
 "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
}