{ "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": [ ":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": [ "" ] }, "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 }