path: root/lab-windows/dsss_experiments-ber.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import struct\n",
    "import random\n",
    "import itertools\n",
    "import datetime\n",
    "import multiprocessing\n",
    "from collections import defaultdict\n",
    "import json\n",
    "\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "from scipy import signal as sig\n",
    "import ipywidgets\n",
    "\n",
    "from tqdm.notebook import tqdm\n",
    "import colorednoise\n",
    "\n",
    "np.set_printoptions(linewidth=240)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampling_rate = 10 # sp/s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def modulate(data, nbits=5, pad=True):\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",
    "    signal = (np.multiply(sel, np.tile(data_lsb_centered, (2**nbits-1, 1)).T).flatten() + 1) // 2\n",
    "    if pad:\n",
    "        return np.hstack([ np.zeros(len(mask)), signal, np.zeros(len(mask)) ])\n",
    "    else:\n",
    "        return signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean: 49.98625 len: 166118\n"
     ]
    }
   ],
   "source": [
    "with open('data/raw_freq.bin', 'rb') as f:\n",
    "    mains_noise = np.copy(np.frombuffer(f.read(), dtype='float32'))\n",
    "    print('mean:', np.mean(mains_noise), 'len:', len(mains_noise))\n",
    "    mains_noise -= np.mean(mains_noise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "    \n",
    "    return test_data, signal + noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_duration = 200\n",
    "test_nbits = 5\n",
    "test_signal_amplitude=2.0e-3\n",
    "test_decimation=10\n",
    "\n",
    "for test_signal_amplitude in [2.0e-3, 20e-3, 200e-3, 2]:\n",
    "    test_data = np.random.RandomState(seed=0).randint(0, 2 * (2**test_nbits), test_duration)\n",
    "    #test_data = np.array([0, 1, 2, 3] * 50)\n",
    "    signal = np.repeat(modulate(test_data, test_nbits, pad=False) * 2.0 - 1, test_decimation) * test_signal_amplitude\n",
    "    with open(f'dsss_test_signals/dsss_test_noiseless_{test_signal_amplitude*1000:.0f}mHz.bin', 'wb') as f:\n",
    "        for e in signal:\n",
    "            f.write(struct.pack('<f', e))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "nonlinear_distance = lambda x: 100**(2*np.abs(0.5-x%1)) / (np.abs(x)+3)**2 * (np.clip(np.abs(x), 0, 0.5) * 2)**5\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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "373401133cfe408aa15738e48c58dfaa",
       "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": [
    "plot_distance_func()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "noprint = lambda *args, **kwargs: None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-44-38927d86fcc7>: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()\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2015d1f2951340d49397c6c289096b01",
       "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": [
       "Text(0.5, 1.0, 'Ricker wavelet, w=69 a=7.3')"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "w = 69\n",
    "a = 7.3\n",
    "ax.plot(range(-w//2+1, w//2+1), sig.ricker(w, a))\n",
    "ax.grid()\n",
    "ax.axvline(0, color='orange')\n",
    "ax.set_title(f'Ricker wavelet, w={w} a={a}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-56-8e404350b85f>: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()\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3cd0c25b914b41df835005ef6fb51b34",
       "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()\n",
    "r = list(range(60, 120))\n",
    "ax.plot(r, [sum(sig.ricker(w, a)) for w in r])\n",
    "ax.set_yscale('log')\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-67-2097444ed7d2>: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()\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c5d1587fa0c944378aed31941ab66ff6",
       "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()\n",
    "sw = 256\n",
    "w = sig.ricker(sw, a)\n",
    "r = list(range(1, sw//2 - 10))\n",
    "d = [-sum(w[:i]) - sum(w[-i:]) for i in r]\n",
    "ax.plot([sw-2*x for x in r], d)\n",
    "ax.set_yscale('log')\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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, pol_score_factor=1.0, seed=0, ax=None, print=print, ser_maxshift=3, debug_range=None):\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.05 * bit_period\n",
    "    hole_patching_threshold = 0.01 * 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",
    "    peaks_processed = []\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 = max([max(val.min(), val.max(), key=abs) for _idx, (_th, val) in group], key=abs)\n",
    "        pol_idx = np.argmax(np.bincount([ np.argmax(np.abs(val)) for _idx, (_th, val) in group ]))\n",
    "        peaks_processed.append((pos, pol, pol_idx))\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",
    "        msg = f'peak at {pos} = {pol} idx {pol_idx}: '\n",
    "\n",
    "        if peak_group:\n",
    "            msg += f'continuing previous group: {peak_group[-1]},'\n",
    "            group_start, last_pos, last_pol, peak_pos, last_pol_idx = peak_group[-1]\n",
    "\n",
    "            if abs(pol) > abs(last_pol):\n",
    "                msg += 'larger, '\n",
    "                if ax:\n",
    "                    ax.axvline(pos, color='magenta', alpha=0.5)\n",
    "                peak_group[-1] = (group_start, pos, pol, pos, pol_idx)\n",
    "                \n",
    "            else:\n",
    "                msg += 'smaller, '\n",
    "                if ax:\n",
    "                    ax.axvline(pos, color='blue', alpha=0.5)\n",
    "                peak_group[-1] = (group_start, pos, last_pol, peak_pos, last_pol_idx)\n",
    "        else:\n",
    "            last_pos = None\n",
    "        \n",
    "        if not peak_group or pos - last_pos > peak_group_threshold:\n",
    "            msg += 'terminating, '\n",
    "            if peak_group:\n",
    "                msg += f'previous group: {peak_group[-1]},'\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",
    "            msg += f'new group: {(pos, pos, pol, pos, pol_idx)} '\n",
    "            peak_group.append((pos, pos, pol, pos, pol_idx))\n",
    "            if ax:\n",
    "                ax.axvline(pos, color='cyan', alpha=0.5)\n",
    "        \n",
    "        if debug_range:\n",
    "            low, high = debug_range\n",
    "            if low < pos < high:\n",
    "                print(msg)\n",
    "                print(group)\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 = [ (abs(pol)/avg_peak, [(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: pol_score_eebfactor*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, peaks):\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-1} symbols at {pos}/{last_pos}')\n",
    "                    \n",
    "                    # Hole patching routine\n",
    "                    for i in range(1, delta):\n",
    "                        est_pos = last_pos + (pos - last_pos) / delta * i\n",
    "\n",
    "                        icandidates = [ (ipos, ipol, iidx) for ipos, ipol, iidx in peaks if abs(est_pos - ipos) < hole_patching_threshold ]\n",
    "                        if not icandidates:\n",
    "                            yield None\n",
    "                            continue\n",
    "\n",
    "                        ipos, ipol, iidx = max(icandidates, key = lambda e: abs(e[1]))\n",
    "\n",
    "                        decoded = iidx*2 + (0 if ipol < 0 else 1)\n",
    "                        print(f'interpolating, last_pos={last_pos}, delta={delta}, pos={pos}, est={est_pos} dec={decoded}')\n",
    "                        yield decoded\n",
    "            \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, peaks_processed))\n",
    "    print('decoding [ref|dec]:')\n",
    "    match_result = []\n",
    "    for shift in range(-ser_maxshift, ser_maxshift):\n",
    "        msg = f'=== shift = {shift} ===\\n'\n",
    "        failures = -shift if shift < 0 else 0 # we're skipping the first $shift symbols\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",
    "            if ref is None: # end of signal\n",
    "                break\n",
    "            msg += f'{ref if ref is not None else -1:>3d}|{found if found is not None else -1:>3d} {\"✔\" if ref==found else \"✘\" if found else \" \"}   '\n",
    "            if ref != found:\n",
    "                failures += 1\n",
    "            if i%8 == 7:\n",
    "                msg += '\\n'\n",
    "        match_result.append((failures, msg))\n",
    "    failures, msg = min(match_result, key=lambda e: e[0])\n",
    "    print(msg)\n",
    "    ser = failures/len(test_data)\n",
    "    print(f'Symbol error rate e={ser}: {failures}/{len(test_data)}')\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": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "672fae23dc99473a967c46b4b90d09d7",
       "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 0x7fc38cb743d0>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 9))\n",
    "\n",
    "duration = 64\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",
    "#    print(f'seed={seed:08x} > ser={ser:.5f}')\n",
    "\n",
    "#ax.set_xlim([40000, 43000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fig, ax = plt.subplots(figsize=(12,5))\n",
    "#run_ser_test(ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "985b1935796343d19c27e0b0ae31363c",
       "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": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "418a9dd281534ffa8e0aa5f250424127",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=1.0, bar_style='info', max=1.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nbits=5\n",
      "nbits=6\n",
      "signal_amplitude=0.00015: ser=0.98958 ±0.00737, br=5.85938\n",
      "signal_amplitude=0.00032: ser=0.97396 ±0.00737, br=14.64844\n",
      "signal_amplitude=0.00046: ser=0.97396 ±0.00737, br=14.64844\n",
      "signal_amplitude=0.00010: ser=0.97396 ±0.00737, br=14.64844\n",
      "signal_amplitude=0.00068: ser=0.94792 ±0.01949, br=29.29688\n",
      "signal_amplitude=0.00022: ser=0.96875 ±0.01276, br=17.57812\n",
      "signal_amplitude=0.00100: ser=0.78125 ±0.02552, br=123.04688\n",
      "signal_amplitude=0.00147: ser=0.28646 ±0.03683, br=401.36719\n",
      "signal_amplitude=0.00215: ser=0.06250 ±0.03375, br=527.34375\n",
      "signal_amplitude=0.00316: ser=0.00521 ±0.00737, br=559.57031\n",
      "signal_amplitude=0.00015: ser=0.97917 ±0.01473, br=7.03125\n",
      "signal_amplitude=0.00010: ser=0.98958 ±0.00737, br=3.51562\n",
      "signal_amplitude=0.00022: ser=0.98438 ±0.00000, br=5.27344\n",
      "signal_amplitude=0.00032: ser=0.98438 ±0.00000, br=5.27344\n",
      "signal_amplitude=0.00046: ser=0.97917 ±0.01949, br=7.03125\n",
      "signal_amplitude=0.00068: ser=0.68229 ±0.09051, br=107.22656\n",
      "signal_amplitude=0.00100: ser=0.15104 ±0.05156, br=286.52344\n",
      "signal_amplitude=0.00147: ser=0.01562 ±0.00000, br=332.22656\n",
      "signal_amplitude=0.00215: ser=0.01562 ±0.00000, br=332.22656\n",
      "signal_amplitude=0.00316: ser=0.00000 ±0.00000, br=337.50000\n",
      "scheduled 20 tasks. waiting...\n",
      "done\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc382e9b7f0>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "default_params = dict(\n",
    "        decimation=10,\n",
    "        power_avg_width=2.5,\n",
    "        max_lookahead=6.5)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12, 9))\n",
    "\n",
    "def calculate_ser(v, seed, nbits, thf, reps, duration):\n",
    "    st = np.random.RandomState(seed)\n",
    "    params = dict(default_params)\n",
    "    params['signal_amplitude'] = v\n",
    "    params['nbits'] = nbits\n",
    "    params['threshold_factor'] = thf\n",
    "    sers, brs = [], []\n",
    "    for i in range(reps):\n",
    "        seed = st.randint(0xffffffff)\n",
    "        try:\n",
    "            ser, br = run_ser_test(**params, sample_duration=duration, print=noprint, seed=seed)\n",
    "            sers.append(ser)\n",
    "            brs.append(br)\n",
    "        except Exception as e:\n",
    "            print('got', e, 'seed', seed, 'params', params)\n",
    "            #sers.append(1.0)\n",
    "            #brs.append(0.0)\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",
    "results = {}\n",
    "with tqdm(total = 0) as tq:\n",
    "    with multiprocessing.Pool(multiprocessing.cpu_count()//2) as pool:\n",
    "        for nbits, thf, reps, points, duration in [(5, 4.0, 3, 10, 64), (6, 4.0, 3, 10, 64)]: #[(5, 4.0, 50, 25, 128), (6, 4.0, 25, 25, 64), (7, 5.0, 10, 10, 64), (8, 6.0, 5, 10, 32)]:\n",
    "            print(f'nbits={nbits}')\n",
    "\n",
    "            st = np.random.RandomState(0)\n",
    "\n",
    "            vs = 0.1e-3 * 10 ** np.linspace(0, 1.5, points)\n",
    "            results[nbits] = [ pool.apply_async(calculate_ser, (v, st.randint(0xffffffff), nbits, thf, reps, duration), callback=lambda _res: tq.update(1)) for v in vs ]\n",
    "            tq.total += len(vs)\n",
    "            tq.refresh()\n",
    "        \n",
    "        pool.close()\n",
    "        pool.join()\n",
    "\n",
    "        print(f'scheduled {tq.total} tasks. waiting...')\n",
    "        results = { nbits: [ res.get() for res in series ] for nbits, series in results.items() }\n",
    "        print('done')\n",
    "\n",
    "with open(f'dsss_experiments_res-{datetime.datetime.now():%Y-%m-%d %H:%M:%S}.json', 'w') as f:\n",
    "    json.dump(results, f)\n",
    "    \n",
    "for nbits, res in results.items():\n",
    "    data = np.array(res)\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": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8b9593f48de94ec399fd8e61dc11d1de",
       "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.legend.Legend at 0x7fc3c4c54b20>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 9))\n",
    "\n",
    "#    sers, brs = np.array(sers), np.array(brs)\n",
    "#    ser, std = np.mean(sers), np.std(sers)\n",
    "#    results = { nbits: [ res.get() for res in series ] for nbits, series in results.items() }\n",
    "\n",
    "with open(f'data/dsss_experiments_res-2020-02-19-19-30-05.json', 'r') as f:\n",
    "    results = json.load(f)\n",
    "\n",
    "for nbits, series in results.items():\n",
    "    series = [ [ mean for mean, _std, _msg in reps if mean is not None ] for reps in series ]\n",
    "    sers = np.array([ np.mean(values) for values in series ])\n",
    "    stds = np.array([ np.std(values)  for values in series ])\n",
    "\n",
    "    # FIXME HACK HACK HACK\n",
    "    vs = 0.1e-3 * 10 ** np.linspace(0, 1.5, 25)\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": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "15dc2e1a7cdd460e8cbde0352ea18270",
       "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/fromnumeric.py:3334: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:216: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  ret = _var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:185: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean = um.true_divide(\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:209: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3334: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:216: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  ret = _var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:185: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean = um.true_divide(\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:209: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3334: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:216: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  ret = _var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:185: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean = um.true_divide(\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:209: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3334: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:216: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  ret = _var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:185: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean = um.true_divide(\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:209: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3334: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:216: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  ret = _var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:185: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean = um.true_divide(\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:209: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3334: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:216: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  ret = _var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:185: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean = um.true_divide(\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:209: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3334: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:216: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  ret = _var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:185: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean = um.true_divide(\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:209: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3334: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:216: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  ret = _var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:185: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean = um.true_divide(\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:209: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3334: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:216: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  ret = _var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:185: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean = um.true_divide(\n",
      "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_methods.py:209: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\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",
    "    'data/dsss_experiments_res-2020-02-20-12-18-35.json',\n",
    "    'data/dsss_experiments_res-2020-02-20-12-26-07.json',\n",
    "    'data/dsss_experiments_res-2020-02-20-12-29-02.json'\n",
    "]:\n",
    "    with open(fn, 'r') as f:\n",
    "        results += json.load(f)\n",
    "\n",
    "thfs = [thf for (_nbits, thf, _reps, _points, _duration), series in results]\n",
    "cmap = matplotlib.cm.viridis\n",
    "cm_func = lambda x: cmap((x - min(thfs)) / (max(thfs) - min(thfs)))\n",
    "\n",
    "thf_sers = {}\n",
    "for (nbits, thf, reps, points, duration), 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",
    "    thf_sers[thf] = list(zip(amps, sers, stds))\n",
    "    \n",
    "    l, = ax.plot(amps, np.clip(sers, 0, 1), label=f'thf={thf}', color=cm_func(thf))\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",
    "def plot_base_amp(ax):\n",
    "    base_sers = {}\n",
    "    for thf, sers in thf_sers.items():\n",
    "        base = np.mean([ser for amp, ser, std in sorted(sers)[-2:]])\n",
    "        base_std = np.sqrt(np.mean([std**2 for amp, ser, std in sorted(sers)[-2:]]))\n",
    "        base_sers[thf] = (base, base_std)\n",
    "\n",
    "    x = sorted(base_sers.keys())\n",
    "    y = np.array([ base_sers[thf][0] for thf in x ])\n",
    "    std = np.array([ base_sers[thf][1] for thf in x ])\n",
    "    l = ax.plot(x, y, label='Base amplitude')\n",
    "    ax.fill_between(x, y-std, y+std, color=l[0].get_color(), alpha=0.3)\n",
    "    return l\n",
    "\n",
    "def plot_intercepts(ax, SER_TH = 0.5):\n",
    "    intercepts = {}\n",
    "    for thf, sers in thf_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[thf] = (icp, ic_std)\n",
    "                break\n",
    "            last_amp, last_ser = amp, ser\n",
    "        else:\n",
    "            intercepts[thf] = None, None\n",
    "    \n",
    "    ser_valid = [thf for thf, (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[thf] for thf 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",
    "    \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.grid()\n",
    "    return l\n",
    "\n",
    "l1 = plot_intercepts(intercept_ax)\n",
    "l2 = plot_base_amp(intercept_ax.twinx())\n",
    "intercept_ax.legend(l1 + l2, [l.get_label() for l in l1+l2], loc=4)\n",
    "\n",
    "norm = matplotlib.colors.Normalize(vmin=min(thfs), vmax=max(thfs))\n",
    "cb1 = matplotlib.colorbar.ColorbarBase(cbar_ax, cmap=cmap, norm=norm, orientation='vertical', label=\"Threshold factor\")\n",
    "#fig.colorbar(ticks=thfs, label='Threshold factor', ax=ax)\n",
    "#ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fig.savefig('dsss_prototype_symbol_error_rate_5-8_bit.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e381ce02f19d4b85bbbaa9ea99e65623",
       "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-19-50856f807a2a>:20: 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], 'hspace': 0.4})\n",
    "empty.axis('off')\n",
    "#fig.tight_layout()\n",
    "\n",
    "results = []\n",
    "\n",
    "for fn in [\n",
    "#    'data/dsss_experiments_res-2020-02-20-14-10-13.json',\n",
    "#    'data/dsss_experiments_res-2020-02-20-13-21-57.json',\n",
    "#    'data/dsss_experiments_res-2020-02-20-13-23-47.json',\n",
    "    'data/dsss_experiments_res-2020-02-20-19-51-21.json',\n",
    "    'data/dsss_experiments_res-2020-02-20-20-43-32.json',\n",
    "    'data/dsss_experiments_res-2020-02-20-21-36-42.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 if decimation > 0]\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",
    "    if not decimation > 0:\n",
    "        continue\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",
    "    amps = [ amp*1000 for amp in amps ]\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 [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",
    "tick_decs = sorted(decimations)\n",
    "tick_decs = tick_decs[:4] + tick_decs[5::5]\n",
    "yticks = [np.log10(d) for d in tick_decs]\n",
    "cb1 = matplotlib.colorbar.ColorbarBase(cbar_ax, cmap=cmap, norm=norm, orientation='vertical', ticks=yticks)\n",
    "cb1t = cbar_ax.twinx()\n",
    "cb1t.set_ylim(cbar_ax.get_ylim())\n",
    "cb1t.set_yticks(yticks)\n",
    "\n",
    "cbar_ax.set_yticklabels([f'{d/sampling_rate:.1f}' for d in tick_decs])\n",
    "cbar_ax.set_ylabel(\"chip duration [s]\", labelpad=-70)\n",
    "\n",
    "cb1t.set_yticklabels([f'{d/sampling_rate * 2**nbits:.1f}' for d in tick_decs])\n",
    "cb1t.set_ylabel(\"symbol duration [s]\")\n",
    "\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",
    "    ax.set_xlim([min(x), max(x)])\n",
    "    y = [ v*1000 if v is not None else v for v in y ]\n",
    "    l = ax.plot(x, y, label='Amplitude at SER=0.5 [mHz]', color='orange')\n",
    "    #ax.legend(loc=3)\n",
    "    ax.set_ylabel('Amplitude at SER=0.5 [mHz]')\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",
    "    xticks = [1, 2, 5, 10, 20, 50]\n",
    "    ax.set_xticks(xticks)\n",
    "    ax.set_xticklabels([ f'{x/sampling_rate:.1f}' for x in xticks ])\n",
    "    ax.set_xlim([1, 60])\n",
    "    ax.set_xlabel('chip duration [s]')\n",
    "    \n",
    "    axt = ax.twiny()\n",
    "    axt.set_xlim(ax.get_xlim())\n",
    "    axt.set_xscale('log')\n",
    "    axt.set_xticks(xticks)\n",
    "    axt.set_xticklabels([ f'{x/sampling_rate * 2**nbits:.1f}' for x in xticks ])\n",
    "    axt.set_xlabel('symbol duration [s]')\n",
    "    \n",
    "    return l\n",
    "\n",
    "l1 = plot_intercepts(intercept_ax)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "41.6"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "13 * 2**5  / 10"
   ]
  }
 ],
 "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
}