{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup\n",
    "## Import required packages"
   ]
  },
  {
   "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",
    "import traceback\n",
    "import glob\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib\n",
    "from matplotlib import ticker\n",
    "import numpy as np\n",
    "from scipy import signal as sig\n",
    "from scipy import fftpack as fftpack\n",
    "import ipywidgets\n",
    "from IPython.display import set_matplotlib_formats\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\n",
    "%matplotlib inline\n",
    "set_matplotlib_formats('png', 'pdf')\n",
    "font = {'family' : 'normal',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 6}\n",
    "matplotlib.rc('font', **font)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define mains frequency sampling rate\n",
    "\n",
    "This is the rate of mains frequency measurements, also called \"reporting rate\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampling_rate = 10 # sp/s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Library functions\n",
    "### Gold code generator"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "### Gold code modulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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": "markdown",
   "metadata": {},
   "source": [
    "### Gold code correlator\n",
    "\n",
    "This function, used by the prototype demodulation algorithm below, correlates a signal against all 2^n+1 Gold sequences. Given an input signal of length $k$ it produces an output matrix of dimensions $(2^n+1, k)$ with one column for each shift of the reference Gold sequences w.r.t. the input signal and one row per Gold sequence."
   ]
  },
  {
   "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",
    "    # Our input signal has large DC bias. Remove DC bias to reduce numerical errors during correlation.\n",
    "    sequence -= np.mean(sequence)\n",
    "    \n",
    "    return np.array([np.correlate(sequence, row, mode='full') for row in mask])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read recorded mains frequency data from exported capture file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean: 50.00341 len: 1946174\n"
     ]
    }
   ],
   "source": [
    "with open('data/fmeas_export_ocxo_2day.bin', 'rb') as f:\n",
    "    meas_data = np.copy(np.frombuffer(f.read(), dtype='float32'))\n",
    "    print('mean:', np.mean(meas_data), 'len:', len(meas_data))\n",
    "    meas_data -= np.mean(meas_data)\n",
    "    def mains_noise(n):\n",
    "        last_valid = len(meas_data) - n\n",
    "        start = np.random.randint(last_valid)\n",
    "        return meas_data[start:start+n]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test signal generator\n",
    "\n",
    "This generates deterministically random test data, modulates it using the Gold code modulator, scales it to a given target amplitude and adds noise from recorded data above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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 = mains_noise(len(signal))\n",
    "    \n",
    "    return test_data, signal + noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Signal exporters for hardware testing\n",
    "The following two functions generate test data to test the firmware implementation in software simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def do_export_clean():\n",
    "    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": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def do_export_noisy():\n",
    "    test_duration = 32\n",
    "    test_nbits = 5\n",
    "    test_signal_amplitude=2.0e-3\n",
    "    test_decimation=10\n",
    "    test_signal_amplitude = 200e-3\n",
    "    noise_level = 10e-3\n",
    "\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",
    "    test_data = np.array(range(test_duration))\n",
    "    signal = np.repeat(modulate(test_data, test_nbits, pad=False) * 2.0 - 1, test_decimation) * test_signal_amplitude\n",
    "    noise = colorednoise.powerlaw_psd_gaussian(1, len(signal)*10) * noise_level\n",
    "    noise[-int(1.5*len(signal)):][:len(signal)] += signal\n",
    "\n",
    "    with open(f'dsss_test_signals/dsss_test_noisy_padded.bin', 'wb') as f:\n",
    "        for e in noise:\n",
    "            f.write(struct.pack('<f', e))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The algorithm\n",
    "\n",
    "## First we define some components used in our algorithm.\n",
    "\n",
    "The following function is used to score a new correlation peak against previous peaks. The aim is to assign a high fitness the closer the peak lies to a multiple of one symbol period from the last peak. The first peak is the ideal case, subsequent peaks correspond to dropped symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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(-0.2, 5.5, 10000)\n",
    "    ax.plot(x, nonlinear_distance(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "findfont: Font family ['normal'] not found. Falling back to DejaVu Sans.\n"
     ]
    },
    {
     "data": {
      "application/pdf": "[encoded data, 10956 characters]\n",
      "image/png": "[encoded data, 16000 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_distance_func()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ricker wavelet computations for firmware implementation\n",
    "For our firmware implementation we need a ricker wavelet lookup table. To find out the size of this lookup table, we calculate the truncation error for a given size below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "noprint = lambda *args, **kwargs: None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Ricker wavelet, w=69 a=7.3')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "findfont: Font family ['normal'] not found. Falling back to DejaVu Sans.\n"
     ]
    },
    {
     "data": {
      "application/pdf": "[encoded data, 13980 characters]\n",
      "image/png": "[encoded data, 15168 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/pdf": "[encoded data, 9712 characters]\n",
      "image/png": "[encoded data, 11508 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/pdf": "[encoded data, 10904 characters]\n",
      "image/png": "[encoded data, 12060 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": "markdown",
   "metadata": {},
   "source": [
    "## Demodulation algorithm and testing function\n",
    "\n",
    "The following function contains our prototype demodulation algorithm implementation along with test code applying it to simulated input data. By repeatedly running this function while sweeping parameters we can create plots of our algorithm's performance under various conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "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",
    "    print('>    1 Generating test signal')\n",
    "    # Generate test data for this simulation run\n",
    "    test_data, signal = generate_test_signal(sample_duration, nbits, signal_amplitude, decimation, seed)\n",
    "    \n",
    "    # === Begin of our prototype demodulation algorithm. ===\n",
    "    # (1) Correlate the input signal against all 2^n+1 gold sequences using the correlation function defined above.\n",
    "    cor_an = correlate(signal, nbits=nbits, decimation=decimation)\n",
    "\n",
    "    # span to compute average power measurements for peak finding over, in samples\n",
    "    power_avg_width = int(power_avg_width * (2**nbits - 1) * decimation)\n",
    "\n",
    "    bit_period = (2**nbits) * decimation # duration of one DSSS symbol\n",
    "    peak_group_threshold = 0.05 * bit_period # Duration over which to group several detected peaks into one, in samples\n",
    "    hole_patching_threshold = 0.01 * bit_period # Duration over which to ignore temporary dropouts in signal during grouping\n",
    "    \n",
    "    print('>    2 CWT')\n",
    "    # (2) Calculate continuous wavelet transform of correlator output and a ricker wavelet. The width factor of the ricker wavelet at 0.73*decimation was\n",
    "    # determined empirically. This transform acts like a sharpening filter on our peaks and will increase SNR for the following thresholding step.\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",
    "    # (3) Threshold CWT'ed correlator outputs using the factors defined above. Classify any sample as part of a peak that is threshold_factor times\n",
    "    # larger than the average of the surrounding signal. \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",
    "    # Helper function for thresholding\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",
    "    print('>    3 Peak grouping')\n",
    "    # (4) Group samples above threshold value into spans\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 each span of samples above threshold, try to coalesce this span with adjacent spans if it is close enough\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",
    "    print('>    4 Maximum-likelihood decoding')\n",
    "    # Calculate average magnitude of all found peaks for normalization in next step\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",
    "    # (5) Perform Maximum likelihood estimation to group peaks into chains of peaks spaced at one-symbol duration intervals.\n",
    "    noprint = lambda *args, **kwargs: None\n",
    "    def mle_decode(peak_groups, print=print):\n",
    "        \"\"\" Maximum likelihood estimation decoding.\n",
    "        \n",
    "        This function tries to find sequences of peaks that are spaced at one-symbol intervals. It will yield candidate sequences along with their by fitness.\n",
    "        A sequence is evaluated better the higher its peaks, the closer they match one-symbol intervals from each other and the longer it is. \"\"\"\n",
    "        # For each peak, extract index inside capture (in samples), polarity and the index of the Gold sequence that matched to produce this peak\n",
    "        peak_groups = [ (pos, pol, idx) for _1, _2, pol, pos, idx in peak_groups ]\n",
    "        \n",
    "        # Initially populate candidate array with all peaks in first couple of symbol periods from signal start.\n",
    "        candidates = [ (abs(pol)/avg_peak, [(pos, pol, idx)]) for pos, pol, idx in peak_groups if pos < bit_period*2.5 ]\n",
    "\n",
    "        # Iterate while there are candidates remaining\n",
    "        while candidates:\n",
    "            chain_candidates = [] # candidates for next iteration\n",
    "            for chain_score, chain in candidates:\n",
    "                pos, ampl, _idx = chain[-1]\n",
    "                score_fun = lambda pos, npos, npol: pol_score_factor*abs(npol)/avg_peak + nonlinear_distance((npos-pos)/bit_period)\n",
    "                \n",
    "                # For this candidate, consider all peaks that might extend it to a longer chain up to some maximum timestamp in the future.\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",
    "                    # If we have hit the end of our signal or if we did not detect any more peaks, calculate this candidate's score and yield it to the caller.\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",
    "                    # If we have not yet hit the end of our signal, and we still have peaks left try to extend the current candidate with each of these peaks in turn.\n",
    "                    # Calculate the score of the resulting extended chains and if they are not too bad, append them to the candidates for the next iteration\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",
    "            \n",
    "            # For next iteration select top-n highest scoring candidates just computed\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",
    "    # Group peaks into chains and select the chain with the highest score\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",
    "    print('>    5 Matching')\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": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "751dcfcafb7d44e6aa6d77ae422e71c7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "0it [00:00, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scheduled 50 tasks. waiting...\n",
      "done\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f916f444c70>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 17604 characters]\n",
      "image/png": "[encoded data, 48932 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "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",
    "            traceback.print_exc()\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",
    "st = np.random.RandomState(0)\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, 5, 25, 64), (6, 4.0, 5, 25, 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",
    "            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*1e3, np.clip(sers, 0, 1), label=f'{nbits} bit')\n",
    "    ax.fill_between(vs*1e3, 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": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ea39fcd66b66478c8a0df52afdf9709f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "0it [00:00, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting 0.0001/0\n",
      "Starting 0.000316227766016838/0\n",
      ">    1 Generating test signalStarting 0.0014677992676220691/0\n",
      "Starting 0.00046415888336127784/0\n",
      "Starting 0.001/0\n",
      "Starting 0.0021544346900318834/0\n",
      "Starting 0.0006812920690579612/0\n",
      ">    1 Generating test signalStarting 0.0031622776601683794/0\n",
      "\n",
      ">    1 Generating test signal>    1 Generating test signal>    1 Generating test signal>    1 Generating test signal>    1 Generating test signal>    1 Generating test signal\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "Starting 0.00014677992676220695/0\n",
      ">    1 Generating test signal\n",
      "Starting 0.0002154434690031884/0\n",
      ">    1 Generating test signal\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-18-cdc407ab1b04>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     44\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     45\u001b[0m         \u001b[0mpool\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 46\u001b[0;31m         \u001b[0mpool\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     47\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     48\u001b[0m         \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf'scheduled {tq.total} tasks. waiting...'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/multiprocessing/pool.py\u001b[0m in \u001b[0;36mjoin\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    660\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_state\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mCLOSE\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTERMINATE\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    661\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"In unknown state\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 662\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_worker_handler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    663\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_task_handler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    664\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_result_handler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/threading.py\u001b[0m in \u001b[0;36mjoin\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m   1031\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1032\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mtimeout\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1033\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_wait_for_tstate_lock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1034\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1035\u001b[0m             \u001b[0;31m# the behavior of a negative timeout isn't documented, but\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/threading.py\u001b[0m in \u001b[0;36m_wait_for_tstate_lock\u001b[0;34m(self, block, timeout)\u001b[0m\n\u001b[1;32m   1047\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mlock\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# already determined that the C code is done\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1048\u001b[0m             \u001b[0;32massert\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_is_stopped\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1049\u001b[0;31m         \u001b[0;32melif\u001b[0m \u001b[0mlock\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mblock\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1050\u001b[0m             \u001b[0mlock\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrelease\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1051\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_stop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    },
    {
     "data": {
      "application/pdf": "[encoded data, 7304 characters]\n",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAC/CAYAAADjJVfXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAHsklEQVR4nO3bv4tdeRnH8c+zpskEScCMBov1R6G73W52yoWtspZKsImCFsIUImHZ7dxEsq7YiCKWwWItAinTLlhZz4AEK/8CYRQFSbGFeSwyMNmbmdyb3Zm5eXZfr+aS75zDffLl5j0nh3OruwPAXC+sewAAPh0hBxhOyAGGE3KA4YQcYDghBxhuacir6ltV9UFVfe+xtR9X1TtV9YsTnQ6ApZaGvLv/nuSDheVXuvu3SVJVF45/LABWdeZTnn/ot4mqajvJdpKcO3futZdeeulTvg3A58vu7u4/u3tzlWNr2Tc7q+pSkhtJzib5S5IPk3wnycUkG939/tPO39ra6p2dnVVmAWBfVe1299Yqxy69Iu/ufyT52cLynz7JYAAcP0+tAAwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMNyZZQdU1eUkV5NsJLnZ3Q+q6t0kHyV5ubt/csIzAvAUq1yRX0tyK8m9JFf21x4m+WqSB4edUFXbVbVTVTt7e3vHMCYAR1n11kovvFZ3v53kQVWdf+Lg7tvdvdXdW5ubm8cxJwBHWHprJcndPLoi30hyv6ouJdmoqhtJzif578mNB8AyS0Pe3btJdheWb5zMOAA8K0+tAAwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcEIOMJyQAwwn5ADDCTnAcGeWHVBVl5NcTbKR5GZ3P6iqN5O8muRf3f3HE54RgKdY5Yr8WpJbSe4lubK/9qMkD1c8H4ATtGqIe+H1K939myTfqKovLR5cVdtVtVNVO3t7e8cxJwBHWCXkd/Poivy7SS5U1aUkd6rqRpKzSf69eEJ33+7ure7e2tzcPM55AViw9B55d+8m2V1Y/uBEpgHgmbnHDTCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBLQ15Vl6vqV1X1u6o699j69ar6/YlOB8BSq1yRX0tyK8m9JFeSpKp+kOTPJzYVACs7s+JxvfD6epIvJ3m1qja7e+/xg6tqO8l2krz44ovHMScAR1gl5Hfz6Ip8I8n9qrrU3T9Nkqr6+mLEk6S7bye5nSRbW1u9+HMAjs/SkHf3bpLdI3721nEPBMCz8dQKwHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAwwk5wHBCDjCckAMMJ+QAw51ZdkBVXU5yNclGkpvd/aCqfpjka0kudff1E54RgKdY5Yr8WpJbSe4luZIk3X2nu3+d5IuHnVBV21W1U1U7e3t7xzQqAIdZ9dZKP/5aVS9U1XtJ/nDowd23u3uru7c2NzePYUwAjrJKyO/m0RX5d5NcqKpLeRTwi0neqKovnNx4ACyz9B55d+8m2V1Y/tnJjAPAs/LUCsBwQg4wnJADDCfkAMMJOcBwQg4wnJADDCfkAMMJOcBwQg4wnJADDCfkAMMJOcBwQg4wnJADDCfkAMMJOcBwQg4wnJADDCfkAMMJOcBwQg4wnJADDCfkAMMJOcBwQg4wnJADDHdm2QFVdTnJ1SQbSW5294Oq+nGSi0nOdfcvT3hGAJ5ilSvya0luJbmX5Mr+2ivd/dskqaoLJzEYAKtZekW+rxdeF9c/pqq2k2zv//GjqvrbJ5jts+hikn+ue4jnhL04YC8O2IsD3171wOo+tMUHB1S9luR7eXRr5X6SD5N8J482fKO7319y/k53b6060GeZvThgLw7YiwP24sCz7MXSK/Lu3k2yu7D8p08yGADHz1MrAMOdRshvn8J7TGEvDtiLA/bigL04sPJeLL1HDsDzza0VgOFWffxwJb48dOCIvfhhkq8ludTd19c64Ck6bC/2168n+WZ3v7XG8U7VEZ+LN5O8muRf3f3HtQ54io7Yi3eTfJTk5e7+yVoHPGVV9a0kP09yr7vv7a+t1M/jviL35aEDT+xFd9/p7l8n+eL6xlqLJ/aiqn6Q5M/rG2ltDvs38qMkD/P5+x/yYXvxMMlXkzxYz0jr091/T/LBwvJK/TyJD84zfXnoM+5je1FVL1TVe0n+sL6R1mbxc/F6kjeTvFpVm+sZaW0W9+Ir3f2bJN+oqi+taaZ1WdyL6u63kzyoqvNrmul59NR+HuutlSR38+g37EaS+1V1Kclfq+qdJOnu/xzz+z3PDtuLG0kqyRtVdb+7/7fG+U7TE3vR3T9Nkqr6enfvrXO4U3bY5+JOVd1IcjbJv9c422k7bC829vfifJL/rnG2U7f/9/9+krP7v8Q+zIr99NQKwHCft3tyAJ85Qg4wnJADDCfkAMMJOcBwQg4wnJADDPd/rk6zU3aQufoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "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=(6, 3))\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",
    "            print(f'Starting {v}/{i}\\n', end='', flush=True)\n",
    "            ser, br = run_ser_test(**params, sample_duration=duration, print=print, seed=seed)\n",
    "            sers.append(ser)\n",
    "            brs.append(br)\n",
    "            print(f'Finishing {v}/{i}\\n', end='', flush=True)\n",
    "        except Exception as e:\n",
    "            traceback.print_exc()\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",
    "st = np.random.RandomState(0)\n",
    "with tqdm(total = 0) as tq:\n",
    "    with multiprocessing.Pool(multiprocessing.cpu_count()) as pool:\n",
    "        for nbits, thf, reps, points, duration in [(10, 4.0, 5, 10, 10)]:\n",
    "            #print(f'nbits={nbits}')\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*1e3, np.clip(sers, 0, 1), label=f'{nbits} bit')\n",
    "    ax.fill_between(vs*1e3, 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": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/pdf": "[encoded data, 19436 characters]\n",
      "image/png": "[encoded data, 18492 characters]\n",
      "text/plain": [
       "<Figure size 216x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3, 2))\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.1 * 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()\n",
    "fig.tight_layout()\n",
    "fig.savefig('fig_out/dsss_gold_nbits_overview.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/pdf": "[encoded data, 57604 characters]\n",
      "image/png": "[encoded data, 100824 characters]\n",
      "text/plain": [
       "<Figure size 504x288 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "default_files = [\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",
    "        'data/dsss_experiments_res-par107-run115-0-2020-04-07-11-41-37.json',\n",
    "        'data/dsss_experiments_res-par107-run115-1-2020-04-07-13-23-42.json',\n",
    "        'data/dsss_experiments_res-par107-run115-2-2020-04-07-08-57-38.json',\n",
    "        'data/dsss_experiments_res-par107-run115-3-2020-04-07-15-48-04.json',\n",
    "    ]\n",
    "\n",
    "def load_results(*files):\n",
    "    results = []\n",
    "    for fn in files:\n",
    "        with open(fn, 'r') as f:\n",
    "            results += json.load(f)\n",
    "    return results\n",
    "\n",
    "def thf_dependence_plot(results, plot_nbits=6,\n",
    "                        ax=None, cbar_ax=None, intercept_ax=None,\n",
    "                        xlabel=True, x2label=False, ylabel=True, y2label=True, y2ticks=True, ic_ylim=[-0.5, 5],\n",
    "                        legend_loc=4, split_legend=False, calc_best_ampl=False):\n",
    "\n",
    "    thfs = [thf for (_nbits, thf, _reps, _points, _duration, _decimation), 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, decimation), series in results:\n",
    "        if nbits != plot_nbits:\n",
    "            continue\n",
    "        data = [ [ mean for mean, _std, _msg in reps if mean is not None ] for _amp, reps in series ]\n",
    "        amps = [ amp*1000 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",
    "        if ax:\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",
    "    if ax:\n",
    "        ax.grid()\n",
    "        ax.set_title(f'{plot_nbits}-bit Gold code')\n",
    "        if xlabel:\n",
    "            ax.set_xlabel('Amplitude [mHz]')\n",
    "        else:\n",
    "            ax.xaxis.set_ticklabels([])\n",
    "        if ylabel:\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='SER at large amplitudes')\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",
    "        if ax:\n",
    "            ax.set_xlim([min(x), max(x)])\n",
    "            l = ax.plot(x, y, label='Amplitude at SER=0.5', color='orange')\n",
    "        else:\n",
    "            l = None\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",
    "        if ax:\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",
    "        \n",
    "        best_ampl = (np.inf, np.nan)\n",
    "        for yval, stdval in zip(y, std):\n",
    "            if yval < best_ampl[0]:\n",
    "                best_ampl = [yval, stdval]\n",
    "        \n",
    "        return l, best_ampl\n",
    "\n",
    "    if intercept_ax:\n",
    "        if isinstance(intercept_ax, tuple):\n",
    "            intercept_ax, intercept_ax_twin = intercept_ax\n",
    "        else:\n",
    "            intercept_ax_twin = intercept_ax.twinx()\n",
    "            \n",
    "    if intercept_ax or calc_best_ampl:\n",
    "        l1, best_ampl = plot_intercepts(intercept_ax)\n",
    "    else:\n",
    "        best_ampl = None\n",
    "        \n",
    "    if intercept_ax:\n",
    "        l2 = plot_base_amp(intercept_ax_twin)\n",
    "        \n",
    "        intercept_ax.set_title(f'{plot_nbits}-bit Gold code')\n",
    "        if xlabel:\n",
    "            intercept_ax.set_xlabel('Threshold factor')\n",
    "        if x2label:\n",
    "            intercept_ax_twin.set_xlabel('Threshold factor')\n",
    "        if ylabel:\n",
    "            intercept_ax.set_ylabel('Amplitude [mHz]')\n",
    "        intercept_ax.set_ylim(ic_ylim)\n",
    "        intercept_ax_twin.set_ylim([-0.1, 1])\n",
    "        if y2label:\n",
    "            intercept_ax_twin.set_ylabel('Symbol Error Rate')\n",
    "        if not y2ticks:\n",
    "            intercept_ax_twin.set_yticklabels([])\n",
    "        if legend_loc is not None:\n",
    "            if split_legend:\n",
    "                intercept_ax.legend(l1, [l1[0].get_label()], loc=legend_loc)\n",
    "                intercept_ax_twin.legend(l2, [l2[0].get_label()], loc=legend_loc)\n",
    "            else:\n",
    "                intercept_ax.legend(l1 + l2, [l.get_label() for l in l1+l2], loc=legend_loc)\n",
    "\n",
    "    if cbar_ax:\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",
    "        \n",
    "    return best_ampl\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "\n",
    "def plot_gold_sensitivity(results, nbitses=[5,6,7,8]):\n",
    "    nbitses = np.array(nbitses)\n",
    "    ampls = np.array([ thf_dependence_plot(plot_nbits=dep_n, results=results, calc_best_ampl=True) for dep_n in nbitses ])\n",
    "    fig, ax = plt.subplots(figsize=(3, 2))\n",
    "    l = ax.plot(nbitses, ampls[:,0])\n",
    "    ax.fill_between(nbitses, ampls[:,0]-ampls[:,1], ampls[:,0]+ampls[:,1], color=l[0].get_color(), alpha=0.3)\n",
    "    ax.grid()\n",
    "    ax.set_xlabel('Gold code bits')\n",
    "    ax.set_ylabel('Amplitude at SER=0.5 [mHz]')\n",
    "    ax.set_ylim([0, 2])\n",
    "    ax.xaxis.set_major_locator(ticker.MultipleLocator(1.0))\n",
    "    fig.tight_layout()\n",
    "    return fig\n",
    "\n",
    "def plot_amplitude_ber(results, grid=(2, 3), nbitses=[5,6,7,8], figsize=(12, 9), xlim=None, xlog=False):\n",
    "    fig = plt.figure(figsize=figsize)\n",
    "    gs = plt.GridSpec(*grid, figure=fig, width_ratios=[1, 1, 0.05])\n",
    "    \n",
    "    cbar_ax = fig.add_subplot(gs[0, 2])\n",
    "    \n",
    "    axs = np.empty([2, 2], dtype=object)\n",
    "    for i, nbits in enumerate(nbitses):\n",
    "        row, col = i//2, i%2\n",
    "        \n",
    "        ax = axs[row, col] = fig.add_subplot(gs[row, col])\n",
    "        if xlog:\n",
    "            ax.set_xscale('log')\n",
    "        if xlim is not None:\n",
    "            ax.set_xlim(xlim)\n",
    "        if row == 1:\n",
    "            ax.get_shared_x_axes().join(axs[0, col])\n",
    "        if col == 1:\n",
    "            ax.get_shared_y_axes().join(axs[row, 0])\n",
    "            \n",
    "        xlabel = row==1 if len(nbitses) > 2 else True\n",
    "        thf_dependence_plot(plot_nbits=nbits, ax=ax, cbar_ax=cbar_ax if i==0 else None, xlabel=xlabel, ylabel=col==0, results=results)\n",
    "        \n",
    "    return fig\n",
    "\n",
    "def plot_thf_graph(results, nbitses=[5,6,7,8], ic_ylim=[-0.5, 5], figsize=(12, 9)):\n",
    "    fig, axs = plt.subplots(2, 2, figsize=figsize, sharex='col', sharey='row', gridspec_kw={'wspace': 0.1})\n",
    "    for nbits, ax, ax_below in zip(nbitses, axs.flatten(), [*axs.flatten()[2:], None, None]):\n",
    "        if len(nbitses) <= 2:\n",
    "            ax = ax, ax_below\n",
    "            ax_below.grid()\n",
    "            y2label = nbits in [5, 7]\n",
    "            legend_loc = 9\n",
    "            y2ticks = True\n",
    "        else:\n",
    "            y2ticks = y2label = nbits in [6, 8]\n",
    "            legend_loc = 1\n",
    "        thf_dependence_plot(plot_nbits=nbits, intercept_ax=ax,\n",
    "                            xlabel=nbits in [7, 8], ylabel=nbits in [5, 7], y2label=y2label, ic_ylim=ic_ylim,\n",
    "                            y2ticks=y2ticks, x2label=len(nbitses) <= 2,\n",
    "                            legend_loc=legend_loc if nbits == nbitses[-1] else None,\n",
    "                            split_legend =  len(nbitses) <= 2,\n",
    "                            results=results)\n",
    "    return fig\n",
    "\n",
    "#plot_gold_sensitivity(load_results(*default_files))\\\n",
    "#.savefig('fig_out/dsss_gold_nbits_sensitivity.pdf');\n",
    "\n",
    "#plot_amplitude_ber(load_results(*default_files), figsize=(7, 4))\\\n",
    "#.savefig('fig_out/dsss_thf_amplitude_5678.pdf');\n",
    "\n",
    "#plot_thf_graph(load_results(*default_files), figsize=(7, 4))\\\n",
    "#.savefig('fig_out/dsss_thf_sensitivity_5678.pdf')\n",
    "\n",
    "# Note: due to a mistake these \"par114\" files actually contain \"par115\" data.\n",
    "#plot_thf_graph(load_results('data/dsss_experiments_res-par114-run120-0-2020-04-08-20-03-56.json'),\n",
    "#               nbitses=[5,6], ic_ylim=[-4.99, 30],\n",
    "#              figsize=(7, 4))\\\n",
    "#.savefig('fig_out/dsss_thf_sensitivity_56.pdf');\n",
    "\n",
    "#plot_amplitude_ber(load_results('data/dsss_experiments_res-par114-run120-0-2020-04-08-20-03-56.json'), nbitses=[5,6],\n",
    "#                  xlog=True, xlim=[5e-2, 5e1],\n",
    "#                  figsize=(7, 4))\\\n",
    "#.savefig('fig_out/dsss_thf_amplitude_56.pdf');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/pdf": "[encoded data, 39504 characters]\n",
      "image/png": "[encoded data, 55244 characters]\n",
      "text/plain": [
       "<Figure size 504x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 29288 characters]\n",
      "image/png": "[encoded data, 47872 characters]\n",
      "text/plain": [
       "<Figure size 504x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 23780 characters]\n",
      "image/png": "[encoded data, 22056 characters]\n",
      "text/plain": [
       "<Figure size 504x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def load_results_fw_sim(*files, background=None, filter_decimation=None):\n",
    "    results = defaultdict(lambda: defaultdict(lambda: defaultdict(lambda: [])))\n",
    "    \n",
    "    for fn in files:\n",
    "        with open(fn, 'r') as f:\n",
    "            for (nbits, thf, decimation, symbols, seed, amp, background), result in json.load(f):\n",
    "                if filter_decimation is None or decimation == filter_decimation:\n",
    "                    results[background][(nbits, thf, symbols, decimation)][amp].append(result)\n",
    "                \n",
    "    if len(results) > 1:\n",
    "        if background is None:\n",
    "            raise ValueError('Results series contains series for multiple noise backgrounds. Please select one.')\n",
    "        \n",
    "        results = results[background]\n",
    "    else:\n",
    "        results = list(results.values())[0]\n",
    "    \n",
    "    out = []\n",
    "    for (nbits, thf, duration, decimation), series in results.items():\n",
    "        out_series = []\n",
    "        for amplitude, amplitude_series in sorted(series.items(), key=lambda x: x[0]):\n",
    "            reps = len(amplitude_series)\n",
    "            out_amplitude_series = [(ser if ser is not None else 1.0, None, None) for ser in amplitude_series]\n",
    "            out_series.append((amplitude, out_amplitude_series))\n",
    "        out.append(((nbits, thf, reps, len(series), duration, decimation), out_series))\n",
    "    return out\n",
    "\n",
    "plot_amplitude_ber(load_results(*default_files), nbitses=[5, 6], figsize=(7, 4))\\\n",
    ".savefig('fig_out/dsss_thf_amplitude_56_jupyter_impl.pdf');\n",
    "\n",
    "fw_sim_res = load_results_fw_sim(*glob.glob('data/fw_sim_ser_2/*.json'), filter_decimation=10)\n",
    "plot_amplitude_ber(results=fw_sim_res, nbitses=[5, 6], figsize=(7, 4))\\\n",
    ".savefig('fig_out/dsss_thf_amplitude_56_fw_impl.pdf');\n",
    "plot_thf_graph(results=fw_sim_res, nbitses=[5, 6], figsize=(7, 4))\\\n",
    ".savefig('fig_out/dsss_thf_sensitivity_56_fw_impl.pdf');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sorted([x[0] for x in fw_sim_res])\n",
    "#sorted({amp for _params, series in fw_sim_res for amp, reps in series}),\\\n",
    "#sorted({amp for _params, series in load_results(*default_files) for amp, reps in series})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "findfont: Font family ['normal'] not found. Falling back to DejaVu Sans.\n"
     ]
    },
    {
     "data": {
      "application/pdf": "[encoded data, 25532 characters]\n",
      "image/png": "[encoded data, 56304 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 25508 characters]\n",
      "image/png": "[encoded data, 55212 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 56912 characters]\n",
      "image/png": "[encoded data, 120140 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 55856 characters]\n",
      "image/png": "[encoded data, 116396 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "extra_amp_files = ['data/dsss_experiments_res-par114-run120-0-2020-04-08-20-03-56.json']\n",
    "synth_files = [\n",
    "    'data/dsss_experiments_res-par115-synth-run122-0-2020-04-11-20-07-33.json',\n",
    "    'data/dsss_experiments_res-par115-synth-run122-1-2020-04-11-20-39-19.json'\n",
    "]\n",
    "\n",
    "# Note: due to a mistake these \"par114\" files actually contain \"par115\" data.\n",
    "plot_thf_graph(extra_amp_files, nbitses=[5,6], ic_ylim=[-4.99, 30],);\n",
    "plot_thf_graph(synth_files, nbitses=[5,6], ic_ylim=[-4.99, 30]);\n",
    "\n",
    "plot_amplitude_ber(extra_amp_files, nbitses=[5,6], xlog=True, xlim=[5e-2, 5e1]);\n",
    "plot_amplitude_ber(synth_files, nbitses=[5,6], xlog=True, xlim=[5e-2, 5e1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/pdf": "[encoded data, 34268 characters]\n",
      "image/png": "[encoded data, 69180 characters]\n",
      "text/plain": [
       "<Figure size 504x360 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 32000 characters]\n",
      "image/png": "[encoded data, 66936 characters]\n",
      "text/plain": [
       "<Figure size 504x360 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "chip_duration_default_files = [\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",
    "#    'data/dsss_experiments_res-par107-run115-0-2020-04-07-11-41-37.json',\n",
    "#    'data/dsss_experiments_res-par107-run115-1-2020-04-07-13-23-42.json',\n",
    "#    'data/dsss_experiments_res-par107-run115-2-2020-04-07-08-57-38.json',\n",
    "#    'data/dsss_experiments_res-par107-run115-3-2020-04-07-15-48-04.json',\n",
    "    'data/dsss_experiments_res-par114-run119-0-2020-04-08-20-13-44.json'\n",
    "]\n",
    "\n",
    "def plot_chip_duration_sensitivity(only_nbits=5, files=chip_duration_default_files, only_thf=None, figsize=(12, 9)):\n",
    "    fig, ((ax, cbar_ax), (intercept_ax, empty)) = plt.subplots(2, 2, figsize=figsize, gridspec_kw={'width_ratios': [1, 0.05], 'hspace': 0.5})\n",
    "    empty.axis('off')\n",
    "    #fig.tight_layout()\n",
    "\n",
    "    results = []\n",
    "\n",
    "    for fn in files:\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 only_thf is not None and thf != only_thf:\n",
    "            continue\n",
    "        if nbits != only_nbits:\n",
    "            continue\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",
    "    ax.set_title(f'{only_nbits}-bit Gold code')\n",
    "\n",
    "    norm = matplotlib.colors.Normalize(vmin=np.log10(min(decimations)), vmax=np.log10(max(decimations)))\n",
    "    tick_decs = sorted(set(float(dec) for dec in decimations))\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=-40)\n",
    "\n",
    "    cb1t.set_yticklabels([f'{d/sampling_rate * 2**only_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**only_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",
    "    return fig\n",
    "\n",
    "plot_chip_duration_sensitivity(5, figsize=(7, 5))\\\n",
    ".savefig('fig_out/chip_duration_sensitivity_5.pdf');\n",
    "plot_chip_duration_sensitivity(6, figsize=(7, 5))\\\n",
    ".savefig('fig_out/chip_duration_sensitivity_6.pdf');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/pdf": "[encoded data, 34656 characters]\n",
      "image/png": "[encoded data, 123284 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 34112 characters]\n",
      "image/png": "[encoded data, 117040 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 34308 characters]\n",
      "image/png": "[encoded data, 119308 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 32572 characters]\n",
      "image/png": "[encoded data, 120908 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "new_files = [\n",
    "    'data/dsss_experiments_res-par111-run119-0-2020-04-09-04-02-53.json',\n",
    "    'data/dsss_experiments_res-par111-run119-1-2020-04-08-16-11-20.json',\n",
    "    'data/dsss_experiments_res-par111-run119-2-2020-04-08-18-07-22.json',\n",
    "    'data/dsss_experiments_res-par111-run119-3-2020-04-08-13-56-03.json',\n",
    "            ]\n",
    "plot_chip_duration_sensitivity(6, only_thf=3.5, files=new_files);\n",
    "plot_chip_duration_sensitivity(6, only_thf=4.0, files=new_files);\n",
    "plot_chip_duration_sensitivity(6, only_thf=4.5, files=new_files);\n",
    "plot_chip_duration_sensitivity(6, only_thf=5.0, files=new_files);\n",
    "#plot_chip_duration_sensitivity(6, files=new_files);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/pdf": "[encoded data, 34880 characters]\n",
      "image/png": "[encoded data, 125016 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 34760 characters]\n",
      "image/png": "[encoded data, 130120 characters]\n",
      "text/plain": [
       "<Figure size 864x648 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 32112 characters]\n",
      "image/png": "[encoded data, 67540 characters]\n",
      "text/plain": [
       "<Figure size 504x360 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/pdf": "[encoded data, 31784 characters]\n",
      "image/png": "[encoded data, 65240 characters]\n",
      "text/plain": [
       "<Figure size 504x360 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "synth_files = [\n",
    "    'data/dsss_experiments_res-par114-synth-run121-0-2020-04-11-17-50-31.json',\n",
    "    'data/dsss_experiments_res-par114-synth-run121-1-2020-04-12-03-46-19.json',\n",
    "    'data/dsss_experiments_res-par114-synth-run121-2-2020-04-11-18-44-15.json',\n",
    "    'data/dsss_experiments_res-par114-synth-run121-3-2020-04-11-15-25-53.json',\n",
    "            ]\n",
    "plot_chip_duration_sensitivity(5, only_thf=4.0, files=new_files);\n",
    "plot_chip_duration_sensitivity(5, only_thf=4.0, files=synth_files);\n",
    "plot_chip_duration_sensitivity(6, only_thf=5.0, files=new_files, figsize=(7, 5))\\\n",
    ".savefig('fig_out/chip_duration_sensitivity_cmp_meas_6.pdf');\n",
    "plot_chip_duration_sensitivity(6, only_thf=5.0, files=synth_files, figsize=(7, 5))\\\n",
    ".savefig('fig_out/chip_duration_sensitivity_cmp_synth_6.pdf');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ma-thesis-env",
   "language": "python",
   "name": "ma-thesis-env"
  },
  "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
