Add BER curve for DSSS experiments

This commit is contained in:
jaseg 2020-02-17 17:43:57 +00:00
parent 61accdf087
commit 9a833efed9
2 changed files with 556 additions and 35 deletions

View file

@ -0,0 +1,505 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [],
"source": [
"from matplotlib import pyplot as plt\n",
"import numpy as np\n",
"from scipy import signal as sig\n",
"import struct\n",
"import random\n",
"import ipywidgets\n",
"import itertools\n",
"from multiprocessing import Pool\n",
"\n",
"import colorednoise\n",
"\n",
"np.set_printoptions(linewidth=240)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib widget"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"sampling_rate = 10 # sp/s"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [],
"source": [
"# From https://github.com/mubeta06/python/blob/master/signal_processing/sp/gold.py\n",
"preferred_pairs = {5:[[2],[1,2,3]], 6:[[5],[1,4,5]], 7:[[4],[4,5,6]],\n",
" 8:[[1,2,3,6,7],[1,2,7]], 9:[[5],[3,5,6]], \n",
" 10:[[2,5,9],[3,4,6,8,9]], 11:[[9],[3,6,9]]}\n",
"\n",
"def gen_gold(seq1, seq2):\n",
" gold = [seq1, seq2]\n",
" for shift in range(len(seq1)):\n",
" gold.append(seq1 ^ np.roll(seq2, -shift))\n",
" return gold\n",
"\n",
"def gold(n):\n",
" n = int(n)\n",
" if not n in preferred_pairs:\n",
" raise KeyError('preferred pairs for %s bits unknown' % str(n))\n",
" t0, t1 = preferred_pairs[n]\n",
" (seq0, _st0), (seq1, _st1) = sig.max_len_seq(n, taps=t0), sig.max_len_seq(n, taps=t1)\n",
" return gen_gold(seq0, seq1)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def modulate(data, nbits=5):\n",
" # 0, 1 -> -1, 1\n",
" mask = np.array(gold(nbits))*2 - 1\n",
" \n",
" sel = mask[data>>1]\n",
" data_lsb_centered = ((data&1)*2 - 1)\n",
"\n",
" return (np.multiply(sel, np.tile(data_lsb_centered, (2**nbits-1, 1)).T).flatten() + 1) // 2"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def correlate(sequence, nbits=5, decimation=1, mask_filter=lambda x: x):\n",
" mask = np.tile(np.array(gold(nbits))[:,:,np.newaxis]*2 - 1, (1, 1, decimation)).reshape((2**nbits + 1, (2**nbits-1) * decimation))\n",
"\n",
" sequence -= np.mean(sequence)\n",
" \n",
" return np.array([np.correlate(sequence, row, mode='full') for row in mask])"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"mean: 49.98625\n"
]
}
],
"source": [
"with open('/mnt/c/Users/jaseg/shared/raw_freq.bin', 'rb') as f:\n",
" mains_noise = np.copy(np.frombuffer(f.read(), dtype='float32'))\n",
" print('mean:', np.mean(mains_noise))\n",
" mains_noise -= np.mean(mains_noise)"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [],
"source": [
"def generate_test_signal(duration, nbits=6, signal_amplitude=2.0e-3, decimation=10, seed=0):\n",
" test_data = np.random.RandomState(seed=seed).randint(0, 2 * (2**nbits), duration)\n",
" \n",
" signal = np.repeat(modulate(test_data, nbits) * 2.0 - 1, decimation) * signal_amplitude\n",
" noise = np.resize(mains_noise, len(signal))\n",
" \n",
" return test_data, signal + noise"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"nonlinear_distance = lambda x: 100**(2*np.abs(0.5-x%1)) / (np.abs(x)+3)**2\n",
"\n",
"def plot_distance_func():\n",
" fig, ax = plt.subplots()\n",
" x = np.linspace(-1.5, 5.5, 10000)\n",
" ax.plot(x, nonlinear_distance(x))"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"noprint = lambda *args, **kwargs: None"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [],
"source": [
"def run_ser_test(sample_duration=128, nbits=6, signal_amplitude=2.0e-3, decimation=10, threshold_factor=4.0, power_avg_width=2.5, max_lookahead=6.5, seed=0, ax=None, print=print):\n",
"\n",
" test_data, signal = generate_test_signal(sample_duration, nbits, signal_amplitude, decimation, seed)\n",
" cor_an = correlate(signal, nbits=nbits, decimation=decimation)\n",
"\n",
" power_avg_width = int(power_avg_width * (2**nbits - 1) * decimation)\n",
"\n",
" bit_period = (2**nbits) * decimation\n",
" peak_group_threshold = 0.1 * bit_period\n",
" \n",
" cwt_res = np.array([ sig.cwt(row, sig.ricker, [0.73 * decimation]).flatten() for row in cor_an ])\n",
" if ax:\n",
" ax.grid()\n",
" ax.plot(cwt_res.T)\n",
" \n",
" th = np.array([ np.convolve(np.abs(row), np.ones((power_avg_width,))/power_avg_width, mode='same') for row in cwt_res ])\n",
"\n",
" def compare_th(elem):\n",
" idx, (th, val) = elem\n",
" #print('compare_th:', th.shape, val.shape)\n",
" return np.any(np.abs(val) > th*threshold_factor)\n",
"\n",
" peaks = [ list(group) for val, group in itertools.groupby(enumerate(zip(th.T, cwt_res.T)), compare_th) if val ]\n",
" peak_group = []\n",
" for group in peaks:\n",
" pos = np.mean([idx for idx, _val in group])\n",
" pol = np.mean([max(val.min(), val.max(), key=abs) for _idx, (_th, val) in group])\n",
" pol_idx = np.argmax(np.bincount([ np.argmax(np.abs(val)) for _idx, (_th, val) in group ]))\n",
" #print(f'group', pos, pol, pol_idx)\n",
" #for pol, (_idx, (_th, val)) in zip([max(val.min(), val.max(), key=abs) for _idx, (_th, val) in group], group):\n",
" # print(' ', pol, val)\n",
" if ax:\n",
" ax.axvline(pos, color='cyan', alpha=0.3)\n",
"\n",
" if not peak_group or pos - peak_group[-1][1] > peak_group_threshold:\n",
" if peak_group:\n",
" peak_pos = peak_group[-1][3]\n",
" if ax:\n",
" ax.axvline(peak_pos, color='red', alpha=0.6)\n",
" #ax3.text(peak_pos-20, 2.0, f'{0 if pol < 0 else 1}', horizontalalignment='right', verticalalignment='center', color='black')\n",
"\n",
" peak_group.append((pos, pos, pol, pos, pol_idx))\n",
" #ax3.axvline(pos, color='cyan', alpha=0.5)\n",
"\n",
" else:\n",
" group_start, last_pos, last_pol, peak_pos, last_pol_idx = peak_group[-1]\n",
"\n",
" if abs(pol) > abs(last_pol):\n",
" #ax3.axvline(pos, color='magenta', alpha=0.5)\n",
" peak_group[-1] = (group_start, pos, pol, pos, pol_idx)\n",
" else:\n",
" #ax3.axvline(pos, color='blue', alpha=0.5)\n",
" peak_group[-1] = (group_start, pos, last_pol, peak_pos, last_pol_idx)\n",
"\n",
" avg_peak = np.mean(np.abs(np.array([last_pol for _1, _2, last_pol, _3, _4 in peak_group])))\n",
" print('avg_peak', avg_peak)\n",
"\n",
" noprint = lambda *args, **kwargs: None\n",
" def mle_decode(peak_groups, print=print):\n",
" peak_groups = [ (pos, pol, idx) for _1, _2, pol, pos, idx in peak_groups ]\n",
" candidates = [ (0, [(pos, pol, idx)]) for pos, pol, idx in peak_groups if pos < bit_period*2.5 ]\n",
"\n",
" while candidates:\n",
" chain_candidates = []\n",
" for chain_score, chain in candidates:\n",
" pos, ampl, _idx = chain[-1]\n",
" score_fun = lambda pos, npos, npol: abs(npol)/avg_peak + nonlinear_distance((npos-pos)/bit_period)\n",
" next_candidates = sorted([ (score_fun(pos, npos, npol), npos, npol, nidx) for npos, npol, nidx in peak_groups if pos < npos < pos + bit_period*max_lookahead ], reverse=True)\n",
"\n",
" print(f' candidates for {pos}, {ampl}:')\n",
" for score, npos, npol, nidx in next_candidates:\n",
" print(f' {score:.4f} {npos:.2f} {npol:.2f} {nidx:.2f}')\n",
"\n",
" nch, cor_len = cor_an.shape\n",
" if cor_len - pos < 1.5*bit_period or not next_candidates:\n",
" score = sum(score_fun(opos, npos, npol) for (opos, _opol, _oidx), (npos, npol, _nidx) in zip(chain[:-1], chain[1:])) / len(chain)\n",
" yield score, chain\n",
"\n",
" else:\n",
" print('extending')\n",
" for score, npos, npol, nidx in next_candidates[:3]:\n",
" if score > 0.5:\n",
" new_chain_score = chain_score * 0.9 + score * 0.1\n",
" chain_candidates.append((new_chain_score, chain + [(npos, npol, nidx)]))\n",
" print('chain candidates:')\n",
" for score, chain in sorted(chain_candidates, reverse=True):\n",
" print(' ', [(score, [(f'{pos:.2f}', f'{pol:.2f}') for pos, pol, _idx in chain])])\n",
" candidates = [ (chain_score, chain) for chain_score, chain in sorted(chain_candidates, reverse=True)[:10] ]\n",
"\n",
" res = sorted(mle_decode(peak_group, print=noprint), reverse=True)\n",
" #for i, (score, chain) in enumerate(res):\n",
" # print(f'Chain {i}@{score:.4f}: {chain}')\n",
" (_score, chain), *_ = res\n",
"\n",
" def viz(chain):\n",
" last_pos = None\n",
" for pos, pol, nidx in chain:\n",
" if last_pos:\n",
" delta = int(round((pos - last_pos) / bit_period))\n",
" if delta > 1:\n",
" print(f'skipped {delta} symbols at {pos}')\n",
" for i in range(delta-1):\n",
" yield None\n",
" decoded = nidx*2 + (0 if pol < 0 else 1)\n",
" yield decoded\n",
" if ax:\n",
" ax.axvline(pos, color='blue', alpha=0.5)\n",
" ax.text(pos-20, 0.0, f'{decoded}', horizontalalignment='right', verticalalignment='center', color='black')\n",
"\n",
" last_pos = pos\n",
"\n",
" decoded = list(viz(chain))\n",
" print('decoding [ref|dec]:')\n",
" failures = 0\n",
" for i, (ref, found) in enumerate(itertools.zip_longest(test_data, decoded)):\n",
" print(f'{ref or -1:>3d}|{found or -1:>3d} {\"✔\" if ref==found else \"✘\" if found else \" \"}', end=' ')\n",
" if ref != found:\n",
" failures += 1\n",
" if i%8 == 7:\n",
" print()\n",
" ser = failures/len(test_data)\n",
" print(f'Symbol error rate e={ser}')\n",
" br = sampling_rate / decimation / (2**nbits) * nbits * (1 - ser) * 3600\n",
" print(f'maximum bitrate r={br} b/h')\n",
" return ser, br"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "9b7546a233fb4b6cb6e8f809250ba768",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"(63,) (63,)\n",
"(63,) (63,)\n",
"avg_peak 1.6845488102985742\n",
"skipped 2 symbols at 10079.0\n",
"skipped 2 symbols at 30870.0\n",
"skipped 2 symbols at 42209.5\n",
"decoding [ref|dec]:\n",
" 44| 44 ✔ 47| 47 ✔ 117|117 ✔ 64| 64 ✔ 67| 67 ✔ 123|123 ✔ 67| 67 ✔ 103|103 ✔ \n",
" 9| 9 ✔ 83| 83 ✔ 21| 21 ✔ 114|114 ✔ 36| 36 ✔ 87| 87 ✔ 70| -1 88| 88 ✔ \n",
" 88| 88 ✔ 12| 12 ✔ 58| 58 ✔ 65| 65 ✔ 102|102 ✔ 39| 39 ✔ 87| 87 ✔ 46| 46 ✔ \n",
" 88| 88 ✔ 81| 81 ✔ 37| 37 ✔ 25| 25 ✔ 77| 77 ✔ 72| 72 ✔ 9| 9 ✔ 20| 20 ✔ \n",
"115|115 ✔ 80| 80 ✔ 115|115 ✔ 69| 69 ✔ 126|126 ✔ 79| 79 ✔ 47| 47 ✔ 64| 64 ✔ \n",
" 82| 82 ✔ 99| 99 ✔ 88| 88 ✔ 49| 49 ✔ 115|115 ✔ 29| 29 ✔ 19| 19 ✔ 19| -1 \n",
" 14| 14 ✔ 39| 39 ✔ 32| 32 ✔ 65| 65 ✔ 9| 9 ✔ 57| 57 ✔ 127|127 ✔ 32| 32 ✔ \n",
" 31| 31 ✔ 74| 74 ✔ 116|116 ✔ 23| 23 ✔ 35| 35 ✔ 126|126 ✔ 75| 75 ✔ 114|114 ✔ \n",
" 55| 55 ✔ 28| -1 34| 34 ✔ -1| -1 ✔ -1| -1 ✔ 36| 36 ✔ 53| 53 ✔ 5| 5 ✔ \n",
" 38| 38 ✔ 104|104 ✔ 116|116 ✔ 17| 17 ✔ 79| 79 ✔ 4| 4 ✔ 105|105 ✔ 42| 42 ✔ \n",
" 58| 58 ✔ 31| 31 ✔ 120|120 ✔ 1| 1 ✔ 65| 65 ✔ 103|103 ✔ 41| 41 ✔ 57| 57 ✔ \n",
" 35| 35 ✔ 102|103 ✘ 119|119 ✔ 11| 11 ✔ 46| 46 ✔ 82| 82 ✔ 91| 91 ✔ -1| -1 ✔ \n",
" 14| 14 ✔ 99| 99 ✔ 53| 53 ✔ 12| 12 ✔ 121|121 ✔ 42| 42 ✔ 84| 84 ✔ 75| 75 ✔ \n",
" 68| 68 ✔ 6| 6 ✔ 68| 68 ✔ 47| 47 ✔ 127|127 ✔ 116|116 ✔ 3| 3 ✔ 76| 76 ✔ \n",
"100|100 ✔ 52| 52 ✔ 104|104 ✔ 78| 78 ✔ 15| 15 ✔ 20| 20 ✔ 99| 99 ✔ 58| 58 ✔ \n",
" 23| 23 ✔ 79| 79 ✔ 13| 13 ✔ 117|117 ✔ 85| 85 ✔ 48| 48 ✔ 49| 49 ✔ 69| 69 ✔ \n",
"Symbol error rate e=0.03125\n",
"maximum bitrate r=326.953125 b/h\n"
]
}
],
"source": [
"fig, ax = plt.subplots(figsize=(12,5))\n",
"run_ser_test(ax=ax)"
]
},
{
"cell_type": "code",
"execution_count": 84,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-84-f36c7b0ffb61>:13: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
" fig, ax = plt.subplots(figsize=(12, 9))\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "b5eb2cd4f4224f75bc3dc73b6143d849",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"nbits=5\n",
"signal_amplitude=0.00029: ser=1.01000 ±0.012207515615390381, br=-5.62500\n",
"signal_amplitude=0.00020: ser=1.01469 ±0.013678792892649557, br=-8.26172\n",
"signal_amplitude=0.00052: ser=1.00406 ±0.018895270572288715, br=-2.28516\n",
"signal_amplitude=0.00043: ser=1.01000 ±0.018817586521655747, br=-5.62500\n",
"signal_amplitude=0.00024: ser=1.01156 ±0.014510233457804875, br=-6.50391\n",
"signal_amplitude=0.00036: ser=1.00687 ±0.014875787794264881, br=-3.86719\n",
"signal_amplitude=0.00032: ser=1.01156 ±0.01837117307087384, br=-6.50391\n",
"signal_amplitude=0.00022: ser=1.01000 ±0.013535254892317322, br=-5.62500\n",
"signal_amplitude=0.00057: ser=1.00281 ±0.012476540486048206, br=-1.58203\n",
"signal_amplitude=0.00039: ser=1.00812 ±0.014402148277253642, br=-4.57031\n",
"signal_amplitude=0.00047: ser=1.00438 ±0.012899854650343935, br=-2.46094\n",
"signal_amplitude=0.00027: ser=1.00937 ±0.014657549249448218, br=-5.27344\n",
"signal_amplitude=0.00063: ser=0.99938 ±0.019253652250936705, br=0.35156\n",
"signal_amplitude=0.00077: ser=0.99156 ±0.03231920868461974, br=4.74609\n",
"signal_amplitude=0.00093: ser=0.95156 ±0.06625442202223185, br=27.24609\n",
"signal_amplitude=0.00112: ser=0.76000 ±0.2099632594348354, br=135.00000\n",
"signal_amplitude=0.00136: ser=0.51375 ±0.30673813139223494, br=273.51562\n",
"signal_amplitude=0.00165: ser=0.39844 ±0.38814210912370745, br=338.37891\n",
"signal_amplitude=0.00070: ser=0.99281 ±0.023688242072809035, br=4.04297\n",
"signal_amplitude=0.00084: ser=0.96375 ±0.050769469787461836, br=20.39062\n",
"signal_amplitude=0.00102: ser=0.91063 ±0.10310321739645179, br=50.27344\n",
"signal_amplitude=0.00124: ser=0.72500 ±0.23567348639059932, br=154.68750\n",
"signal_amplitude=0.00150: ser=0.40969 ±0.3064419041596629, br=332.05078\n",
"signal_amplitude=0.00182: ser=0.32531 ±0.38085840544748384, br=379.51172\n",
"signal_amplitude=0.00200: ser=0.29000 ±0.3885339029608613, br=399.37500\n",
"nbits=6\n",
"signal_amplitude=0.00052: ser=1.00375 ±0.027432445434193427, br=-1.26562\n",
"signal_amplitude=0.00029: ser=1.01531 ±0.013528038013695853, br=-5.16797\n",
"signal_amplitude=0.00020: ser=1.02000 ±0.01698459780212649, br=-6.75000\n",
"signal_amplitude=0.00024: ser=1.01844 ±0.0197494066366562, br=-6.22266\n",
"signal_amplitude=0.00043: ser=1.01000 ±0.013535254892317322, br=-3.37500\n",
"signal_amplitude=0.00036: ser=1.01500 ±0.01860884366369926, br=-5.06250\n",
"signal_amplitude=0.00032: ser=1.00906 ±0.01443601182806387, br=-3.05859\n",
"signal_amplitude=0.00022: ser=1.01656 ±0.015200483133769137, br=-5.58984\n",
"signal_amplitude=0.00057: ser=0.98281 ±0.04926213365760764, br=5.80078\n",
"signal_amplitude=0.00027: ser=1.02000 ±0.015946688527716343, br=-6.75000\n",
"signal_amplitude=0.00047: ser=1.00687 ±0.02815276407388802, br=-2.32031\n",
"signal_amplitude=0.00039: ser=1.00906 ±0.016189792308735775, br=-3.05859\n",
"signal_amplitude=0.00077: ser=0.76906 ±0.23454244018940368, br=77.94141\n",
"signal_amplitude=0.00063: ser=0.94031 ±0.08557822627572974, br=20.14453\n",
"signal_amplitude=0.00112: ser=0.29750 ±0.347296478171029, br=237.09375\n",
"signal_amplitude=0.00093: ser=0.50125 ±0.3293776683952632, br=168.32812\n",
"signal_amplitude=0.00136: ser=0.37250 ±0.42536588111001566, br=211.78125\n",
"signal_amplitude=0.00165: ser=0.51000 ±0.46215950303980546, br=165.37500\n",
"signal_amplitude=0.00070: ser=0.90063 ±0.1848975645458858, br=33.53906\n",
"signal_amplitude=0.00084: ser=0.64687 ±0.26652421325275494, br=119.17969\n",
"signal_amplitude=0.00124: ser=0.38500 ±0.39889079606767064, br=207.56250\n",
"signal_amplitude=0.00102: ser=0.40875 ±0.3467111099315971, br=199.54688\n",
"signal_amplitude=0.00150: ser=0.40375 ±0.4435118198819508, br=201.23438\n",
"signal_amplitude=0.00182: ser=0.58531 ±0.46179168734397985, br=139.95703\n",
"signal_amplitude=0.00200: ser=0.61594 ±0.4584529436730666, br=129.62109\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7fe2bdd64430>"
]
},
"execution_count": 84,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sample_duration=128\n",
"sample_reps = 25\n",
"sweep_points = 25\n",
"\n",
"default_params = dict(\n",
" nbits=6,\n",
" signal_amplitude=2.0e-3,\n",
" decimation=10,\n",
" threshold_factor=4.0,\n",
" power_avg_width=2.5,\n",
" max_lookahead=6.5)\n",
"\n",
"fig, ax = plt.subplots(figsize=(12, 9))\n",
"\n",
"for nbits in [5, 6]: # FIXME make sim stable for higher bit counts\n",
" print(f'nbits={nbits}')\n",
" \n",
" def calculate_ser(v):\n",
" params = dict(default_params)\n",
" params['signal_amplitude'] = v\n",
" params['nbits'] = nbits\n",
" sers, brs = [], []\n",
" for i in range(sample_reps):\n",
" ser, br = run_ser_test(**params, sample_duration=sample_duration, print=noprint, seed=np.random.randint(0xffffffff))\n",
" sers.append(ser)\n",
" brs.append(br)\n",
" sers, brs = np.array(sers), np.array(brs)\n",
" ser = np.mean(sers)\n",
" print(f'signal_amplitude={v:<.5f}: ser={ser:<.5f} ±{np.std(sers):<.5f}, br={np.mean(brs):<.5f}')\n",
" return ser\n",
" \n",
" vs = 0.2e-3 * 10 ** np.linspace(0, 1.0, sweep_points)\n",
" with Pool(6) as p:\n",
" data = p.map(calculate_ser, vs)\n",
" \n",
" ax.plot(vs, data, label=f'{nbits} bit')\n",
"ax.grid()\n",
"ax.set_xlabel('Amplitude in mHz')\n",
"ax.set_ylabel('Symbol error rate')\n",
"ax.legend()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "labenv",
"language": "python",
"name": "labenv"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}

View file

@ -28,6 +28,15 @@
"%matplotlib widget"
]
},
{
"cell_type": "code",
"execution_count": 105,
"metadata": {},
"outputs": [],
"source": [
"sampling_rate = 10 # sp/s"
]
},
{
"cell_type": "code",
"execution_count": 3,
@ -410,29 +419,29 @@
},
{
"cell_type": "code",
"execution_count": 96,
"execution_count": 145,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(31,) (31,)\n",
"(31,) (31,)\n"
"(63,) (63,)\n",
"(63,) (63,)\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-96-b3aae757ccad>:33: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
"<ipython-input-145-babcf8a4e867>:33: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
" fig, ((ax1, ax3), (ax2, ax4)) = plt.subplots(2, 2, figsize=(16, 9))\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "246c19a3c1424e7eb0b675ff060ea5b3",
"model_id": "10aa67d294304f2ba26c8e6d5555d5e6",
"version_major": 2,
"version_minor": 0
},
@ -446,10 +455,10 @@
{
"data": {
"text/plain": [
"(0.002, 0.013899708)"
"(0.002, 0.014074279)"
]
},
"execution_count": 96,
"execution_count": 145,
"metadata": {},
"output_type": "execute_result"
}
@ -457,15 +466,15 @@
"source": [
"decimation = 10\n",
"signal_amplitude = 2.0e-3\n",
"nbits = 5\n",
"nbits = 6\n",
"\n",
"#test_data = np.random.randint(0, 2, 100)\n",
"#test_data = np.array([0, 1, 0, 0, 1, 1, 1, 0])\n",
"#test_data = np.random.RandomState(seed=0).randint(0, 2 * (2**nbits), 64)\n",
"test_data = np.random.RandomState(seed=0).randint(0, 2 * (2**nbits), 64)\n",
"#test_data = np.random.RandomState(seed=0).randint(0, 8, 64)\n",
"#test_data = np.array(list(range(8)) * 8)\n",
"#test_data = np.array([0, 1] * 32)\n",
"test_data = np.array(list(range(64)))\n",
"#test_data = np.array(list(range(64)))\n",
"\n",
"foo = np.repeat(modulate(test_data, nbits) * 2.0 - 1, decimation) * signal_amplitude\n",
"noise = np.resize(mains_noise, len(foo))\n",
@ -493,12 +502,20 @@
"ax1.plot(foo + noise)\n",
"ax1.plot(foo)\n",
"ax1.set_title('raw')\n",
"ax1.grid(axis='y')\n",
"\n",
"ax2.plot(filtered)\n",
"ax2.plot(foo)\n",
"ax2.set_title('filtered')\n",
"ax2.grid(axis='y')\n",
"\n",
"ax3.plot(cor1.T)\n",
"for i in range(0, len(foo) + 1, decimation*(2**nbits - 1)):\n",
" ax1.axvline(i, color='gray', alpha=0.5, lw=1)\n",
" ax2.axvline(i, color='gray', alpha=0.5, lw=1)\n",
"\n",
"for i, (color, trace) in enumerate(zip(plt.cm.winter(np.linspace(0, 1, cor1.shape[0])), cor1.T)):\n",
" if i%3 == 0:\n",
" ax3.plot(trace + 0.5 * i, alpha=1.0, color=color)\n",
"ax3.set_title('corr raw')\n",
"ax3.grid()\n",
"\n",
@ -678,21 +695,21 @@
},
{
"cell_type": "code",
"execution_count": 97,
"execution_count": 146,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-97-2d2c2f814215>:11: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
"<ipython-input-146-badd40342f73>:11: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
" fig, (ax1, ax3) = plt.subplots(2, figsize=(12, 5))\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "b083c661b5b441d6b7fc45201faa0576",
"model_id": "509bf67d93b74741b48bca58529c4b9d",
"version_major": 2,
"version_minor": 0
},
@ -707,26 +724,24 @@
"name": "stdout",
"output_type": "stream",
"text": [
"cor_an (33, 20149)\n",
"cwt_res (33, 20149)\n",
"th (33, 20149)\n",
"[((33,), (33,)), ((33,), (33,)), ((33,), (33,)), ((33,), (33,)), ((33,), (33,))]\n",
"peaks: 180\n",
"avg_peak 1.058897833824206\n",
"skipped 3 symbols at 5889.0\n",
"skipped 2 symbols at 8369.0\n",
"skipped 2 symbols at 14568.5\n",
"skipped 2 symbols at 16739.0\n",
"cor_an (65, 40949)\n",
"cwt_res (65, 40949)\n",
"th (65, 40949)\n",
"[((65,), (65,)), ((65,), (65,)), ((65,), (65,)), ((65,), (65,)), ((65,), (65,))]\n",
"peaks: 982\n",
"avg_peak 1.6673786030736735\n",
"skipped 2 symbols at 30238.5\n",
"decoding [ref|dec]:\n",
" -1| -1 ✔ 1| 1 ✔ 2| 2 ✔ 3| 3 ✔ 4| 4 ✔ 5| 5 ✔ 6| 6 ✔ 7| 7 ✔ \n",
" 8| 8 ✔ 9| 9 ✔ 10| 10 ✔ 11| 11 ✔ 12| 12 ✔ 13| 13 ✔ 14| 14 ✔ 15| 15 ✔ \n",
" 16| -1 ✘ 17| -1 ✘ 18| 18 ✔ 19| 19 ✔ 20| 20 ✔ 21| 21 ✔ 22| 22 ✔ 23| 23 ✔ \n",
" 24| 24 ✔ 25| -1 ✘ 26| 26 ✔ 27| 27 ✔ 28| 28 ✔ 29| 29 ✔ 30| 30 ✔ 31| 31 ✔ \n",
" 32| 32 ✔ 33| 33 ✔ 34| 34 ✔ 35| 35 ✔ 36| 36 ✔ 37| 37 ✔ 38| 38 ✔ 39| 39 ✔ \n",
" 40| 40 ✔ 41| 41 ✔ 42| 42 ✔ 43| 43 ✔ 44| 44 ✔ 45| -1 ✘ 46| 46 ✔ 47| 47 ✔ \n",
" 48| 48 ✔ 49| 49 ✔ 50| 50 ✔ 51| 51 ✔ 52| -1 ✘ 53| 53 ✔ 54| 54 ✔ 55| 55 ✔ \n",
" 56| 56 ✔ 57| 57 ✔ 58| 58 ✔ 59| 59 ✔ 60| 60 ✔ 61| 61 ✔ 62| 62 ✔ 63| 56 ✘ \n",
"Symbol error rate r=0.09375\n"
" 44| 44 ✔ 47| 47 ✔ 117|117 ✔ 64| 64 ✔ 67| 67 ✔ 123|123 ✔ 67| 67 ✔ 103|103 ✔ \n",
" 9| 9 ✔ 83| 83 ✔ 21| 21 ✔ 114|114 ✔ 36| 36 ✔ 87| 87 ✔ 70| 70 ✔ 88| 88 ✔ \n",
" 88| 88 ✔ 12| 12 ✔ 58| 58 ✔ 65| 65 ✔ 102|102 ✔ 39| 39 ✔ 87| 87 ✔ 46| 46 ✔ \n",
" 88| 88 ✔ 81| 81 ✔ 37| 37 ✔ 25| 25 ✔ 77| 77 ✔ 72| 72 ✔ 9| 9 ✔ 20| 20 ✔ \n",
"115|115 ✔ 80| 80 ✔ 115|115 ✔ 69| 69 ✔ 126|126 ✔ 79| 79 ✔ 47| 47 ✔ 64| 64 ✔ \n",
" 82| 82 ✔ 99| 99 ✔ 88| 88 ✔ 49| 49 ✔ 115|115 ✔ 29| 29 ✔ 19| -1 19| 19 ✔ \n",
" 14| 14 ✔ 39| 39 ✔ 32| 32 ✔ 65| 64 ✘ 9| 9 ✔ 57| 57 ✔ 127|127 ✔ 32| 32 ✔ \n",
" 31| 31 ✔ 74| 74 ✔ 116|116 ✔ 23| 23 ✔ 35| 35 ✔ 126|126 ✔ 75| 75 ✔ 114| 26 ✘ \n",
"Symbol error rate e=0.046875\n",
"maximum bitrate r=321.6796875 b/h\n"
]
}
],
@ -869,12 +884,13 @@
"print('decoding [ref|dec]:')\n",
"failures = 0\n",
"for i, (ref, found) in enumerate(itertools.zip_longest(test_data, decoded)):\n",
" print(f'{ref or -1:>3d}|{found or -1:>3d} {\"✔\" if ref==found else \"✘\"}', end=' ')\n",
" print(f'{ref or -1:>3d}|{found or -1:>3d} {\"✔\" if ref==found else \"✘\" if found else \" \"}', end=' ')\n",
" if ref != found:\n",
" failures += 1\n",
" if i%8 == 7:\n",
" print()\n",
"print(f'Symbol error rate r={failures/len(test_data)}')\n",
"print(f'Symbol error rate e={failures/len(test_data)}')\n",
"print(f'maximum bitrate r={sampling_rate / decimation / (2**nbits) * nbits * (1 - failures/len(test_data)) * 3600} b/h')\n",
"#ax3.plot(th)"
]
},