1043 lines
47 KiB
Text
1043 lines
47 KiB
Text
{
|
|
"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": 5,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def modulate(data, nbits=5):\n",
|
|
" # 0, 1 -> -1, 1\n",
|
|
" mask = np.array(gold(nbits))*2 - 1\n",
|
|
" \n",
|
|
" sel = mask[data>>1]\n",
|
|
" data_lsb_centered = ((data&1)*2 - 1)\n",
|
|
"\n",
|
|
" signal = (np.multiply(sel, np.tile(data_lsb_centered, (2**nbits-1, 1)).T).flatten() + 1) // 2\n",
|
|
" return np.hstack([ np.zeros(len(mask)), signal, np.zeros(len(mask)) ])"
|
|
]
|
|
},
|
|
{
|
|
"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": 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 = np.resize(mains_noise, len(signal))\n",
|
|
" \n",
|
|
" return test_data, signal + noise"
|
|
]
|
|
},
|
|
{
|
|
"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": "c5fba7c7cec24f09a4f23f7c7d87eb90",
|
|
"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": 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_factor*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": 14,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "5106e62521e246a0b0766cce3f9e556f",
|
|
"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 0x7fc200d359d0>]"
|
|
]
|
|
},
|
|
"execution_count": 14,
|
|
"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": 15,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#fig, ax = plt.subplots(figsize=(12,5))\n",
|
|
"#run_ser_test(ax=ax)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 16,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "15e67925b1bf47b5bbde289b353ed201",
|
|
"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": "6432dbccc7db43d7a227384bca772716",
|
|
"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.00046: ser=0.97396 ±0.00737, br=14.64844\n",
|
|
"signal_amplitude=0.00015: ser=0.98958 ±0.00737, br=5.85938\n",
|
|
"signal_amplitude=0.00022: ser=0.96875 ±0.01276, br=17.57812\n",
|
|
"signal_amplitude=0.00010: ser=0.97396 ±0.00737, br=14.64844\n",
|
|
"signal_amplitude=0.00032: 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.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.00068: ser=0.68229 ±0.09051, br=107.22656\n",
|
|
"signal_amplitude=0.00046: ser=0.97917 ±0.01949, br=7.03125\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 0x7fc23a223550>"
|
|
]
|
|
},
|
|
"execution_count": 16,
|
|
"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": 17,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "bda9f8070b0142928936a7752f318e97",
|
|
"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 0x7fc20046b730>"
|
|
]
|
|
},
|
|
"execution_count": 17,
|
|
"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": 18,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "33b1f777dda149eda3a2f031087372c9",
|
|
"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": 19,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#fig.savefig('dsss_prototype_symbol_error_rate_5-8_bit.svg')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 25,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "406505452bb84597a6c39bc12267dcc6",
|
|
"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-25-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": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"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
|
|
}
|