summaryrefslogtreecommitdiff
path: root/tools
diff options
context:
space:
mode:
Diffstat (limited to 'tools')
-rw-r--r--tools/dsss_experiments-ber.py468
1 files changed, 468 insertions, 0 deletions
diff --git a/tools/dsss_experiments-ber.py b/tools/dsss_experiments-ber.py
new file mode 100644
index 0000000..02a2c50
--- /dev/null
+++ b/tools/dsss_experiments-ber.py
@@ -0,0 +1,468 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[17]:
+
+
+import numpy as np
+import numbers
+import math
+from scipy import signal as sig
+import scipy.fftpack
+import random
+import struct
+import random
+import ipywidgets
+import itertools
+import datetime
+import multiprocessing
+import traceback
+from collections import defaultdict
+import json
+
+from tqdm import tqdm
+import colorednoise
+
+np.set_printoptions(linewidth=240)
+
+
+
+sampling_rate = 10 # sp/s
+
+
+# In[4]:
+
+
+# From https://github.com/mubeta06/python/blob/master/signal_processing/sp/gold.py
+preferred_pairs = {5:[[2],[1,2,3]], 6:[[5],[1,4,5]], 7:[[4],[4,5,6]],
+ 8:[[1,2,3,6,7],[1,2,7]], 9:[[5],[3,5,6]],
+ 10:[[2,5,9],[3,4,6,8,9]], 11:[[9],[3,6,9]]}
+
+def gen_gold(seq1, seq2):
+ gold = [seq1, seq2]
+ for shift in range(len(seq1)):
+ gold.append(seq1 ^ np.roll(seq2, -shift))
+ return gold
+
+def gold(n):
+ n = int(n)
+ if not n in preferred_pairs:
+ raise KeyError('preferred pairs for %s bits unknown' % str(n))
+ t0, t1 = preferred_pairs[n]
+ (seq0, _st0), (seq1, _st1) = sig.max_len_seq(n, taps=t0), sig.max_len_seq(n, taps=t1)
+ return gen_gold(seq0, seq1)
+
+
+# In[5]:
+
+
+def modulate(data, nbits=5):
+ # 0, 1 -> -1, 1
+ mask = np.array(gold(nbits))*2 - 1
+
+ sel = mask[data>>1]
+ data_lsb_centered = ((data&1)*2 - 1)
+
+ signal = (np.multiply(sel, np.tile(data_lsb_centered, (2**nbits-1, 1)).T).flatten() + 1) // 2
+ return np.hstack([ np.zeros(len(mask)), signal, np.zeros(len(mask)) ])
+
+
+# In[6]:
+
+
+def correlate(sequence, nbits=5, decimation=1, mask_filter=lambda x: x):
+ mask = np.tile(np.array(gold(nbits))[:,:,np.newaxis]*2 - 1, (1, 1, decimation)).reshape((2**nbits + 1, (2**nbits-1) * decimation))
+
+ sequence -= np.mean(sequence)
+
+ return np.array([np.correlate(sequence, row, mode='full') for row in mask])
+
+
+# In[7]:
+
+
+def load_noise_meas_params(capture_file):
+ with open(capture_file, 'rb') as f:
+ meas_data = np.copy(np.frombuffer(f.read(), dtype='float32'))
+ print('mean:', np.mean(meas_data), 'len:', len(meas_data))
+ meas_data -= np.mean(meas_data)
+ return {'meas_data': meas_data}
+
+def mains_noise_measured(n, meas_data):
+ last_valid = len(meas_data) - n
+ start = np.random.randint(last_valid)
+ return meas_data[start:start+n]
+
+def load_noise_synth_params(specfile):
+ with open(specfile) as f:
+ d = json.load(f)
+ return {'spl_x': np.linspace(*d['x_spec']),
+ 'spl_N': d['x_spec'][2],
+ 'psd_spl': (d['t'], d['c'], d['k']) }
+
+def mains_noise_synthetic(n, psd_spl, spl_N, spl_x):
+ noise = np.random.normal(size=spl_N) * 2
+ spec = scipy.fftpack.fft(noise) **2
+
+ spec *= np.exp(scipy.interpolate.splev(spl_x, psd_spl))
+
+ spec **= 1/2
+
+ renoise = scipy.fftpack.ifft(spec)
+ return renoise[10000:][:n]
+
+# In[8]:
+
+
+def generate_test_signal(duration, nbits=6, signal_amplitude=2.0e-3, decimation=10, seed=0):
+ test_data = np.random.RandomState(seed=seed).randint(0, 2 * (2**nbits), duration)
+
+ signal = np.repeat(modulate(test_data, nbits) * 2.0 - 1, decimation) * signal_amplitude
+ noise = mains_noise(len(signal))
+
+ return test_data, signal + noise
+
+
+# In[9]:
+
+
+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
+
+# In[11]:
+
+
+noprint = lambda *args, **kwargs: None
+
+
+# In[12]:
+
+
+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):
+
+ test_data, signal = generate_test_signal(sample_duration, nbits, signal_amplitude, decimation, seed)
+ cor_an = correlate(signal, nbits=nbits, decimation=decimation)
+
+ power_avg_width = int(power_avg_width * (2**nbits - 1) * decimation)
+
+ bit_period = (2**nbits) * decimation
+ peak_group_threshold = 0.05 * bit_period
+ hole_patching_threshold = 0.01 * bit_period
+
+ cwt_res = np.array([ sig.cwt(row, sig.ricker, [0.73 * decimation]).flatten() for row in cor_an ])
+ if ax:
+ ax.grid()
+ ax.plot(cwt_res.T)
+
+ th = np.array([ np.convolve(np.abs(row), np.ones((power_avg_width,))/power_avg_width, mode='same') for row in cwt_res ])
+
+ def compare_th(elem):
+ idx, (th, val) = elem
+ #print('compare_th:', th.shape, val.shape)
+ return np.any(np.abs(val) > th*threshold_factor)
+
+ peaks = [ list(group) for val, group in itertools.groupby(enumerate(zip(th.T, cwt_res.T)), compare_th) if val ]
+ peaks_processed = []
+ peak_group = []
+ for group in peaks:
+ pos = np.mean([idx for idx, _val in group])
+ #pol = np.mean([max(val.min(), val.max(), key=abs) for _idx, (_th, val) in group])
+ pol = max([max(val.min(), val.max(), key=abs) for _idx, (_th, val) in group], key=abs)
+ pol_idx = np.argmax(np.bincount([ np.argmax(np.abs(val)) for _idx, (_th, val) in group ]))
+ peaks_processed.append((pos, pol, pol_idx))
+ #print(f'group', pos, pol, pol_idx)
+ #for pol, (_idx, (_th, val)) in zip([max(val.min(), val.max(), key=abs) for _idx, (_th, val) in group], group):
+ # print(' ', pol, val)
+ #if ax:
+ # ax.axvline(pos, color='cyan', alpha=0.3)
+ msg = f'peak at {pos} = {pol} idx {pol_idx}: '
+
+ if peak_group:
+ msg += f'continuing previous group: {peak_group[-1]},'
+ group_start, last_pos, last_pol, peak_pos, last_pol_idx = peak_group[-1]
+
+ if abs(pol) > abs(last_pol):
+ msg += 'larger, '
+ if ax:
+ ax.axvline(pos, color='magenta', alpha=0.5)
+ peak_group[-1] = (group_start, pos, pol, pos, pol_idx)
+
+ else:
+ msg += 'smaller, '
+ if ax:
+ ax.axvline(pos, color='blue', alpha=0.5)
+ peak_group[-1] = (group_start, pos, last_pol, peak_pos, last_pol_idx)
+ else:
+ last_pos = None
+
+ if not peak_group or pos - last_pos > peak_group_threshold:
+ msg += 'terminating, '
+ if peak_group:
+ msg += f'previous group: {peak_group[-1]},'
+ peak_pos = peak_group[-1][3]
+ if ax:
+ ax.axvline(peak_pos, color='red', alpha=0.6)
+ #ax3.text(peak_pos-20, 2.0, f'{0 if pol < 0 else 1}', horizontalalignment='right', verticalalignment='center', color='black')
+
+ msg += f'new group: {(pos, pos, pol, pos, pol_idx)} '
+ peak_group.append((pos, pos, pol, pos, pol_idx))
+ if ax:
+ ax.axvline(pos, color='cyan', alpha=0.5)
+
+ if debug_range:
+ low, high = debug_range
+ if low < pos < high:
+ print(msg)
+ print(group)
+
+ avg_peak = np.mean(np.abs(np.array([last_pol for _1, _2, last_pol, _3, _4 in peak_group])))
+ print('avg_peak', avg_peak)
+
+ noprint = lambda *args, **kwargs: None
+ def mle_decode(peak_groups, print=print):
+ peak_groups = [ (pos, pol, idx) for _1, _2, pol, pos, idx in peak_groups ]
+ candidates = [ (abs(pol)/avg_peak, [(pos, pol, idx)]) for pos, pol, idx in peak_groups if pos < bit_period*2.5 ]
+
+ while candidates:
+ chain_candidates = []
+ for chain_score, chain in candidates:
+ pos, ampl, _idx = chain[-1]
+ score_fun = lambda pos, npos, npol: pol_score_factor*abs(npol)/avg_peak + nonlinear_distance((npos-pos)/bit_period)
+ 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)
+
+ print(f' candidates for {pos}, {ampl}:')
+ for score, npos, npol, nidx in next_candidates:
+ print(f' {score:.4f} {npos:.2f} {npol:.2f} {nidx:.2f}')
+
+ nch, cor_len = cor_an.shape
+ if cor_len - pos < 1.5*bit_period or not next_candidates:
+ score = sum(score_fun(opos, npos, npol) for (opos, _opol, _oidx), (npos, npol, _nidx) in zip(chain[:-1], chain[1:])) / len(chain)
+ yield score, chain
+
+ else:
+ print('extending')
+ for score, npos, npol, nidx in next_candidates[:3]:
+ if score > 0.5:
+ new_chain_score = chain_score * 0.9 + score * 0.1
+ chain_candidates.append((new_chain_score, chain + [(npos, npol, nidx)]))
+ print('chain candidates:')
+ for score, chain in sorted(chain_candidates, reverse=True):
+ print(' ', [(score, [(f'{pos:.2f}', f'{pol:.2f}') for pos, pol, _idx in chain])])
+ candidates = [ (chain_score, chain) for chain_score, chain in sorted(chain_candidates, reverse=True)[:10] ]
+
+ res = sorted(mle_decode(peak_group, print=noprint), reverse=True)
+ #for i, (score, chain) in enumerate(res):
+ # print(f'Chain {i}@{score:.4f}: {chain}')
+ (_score, chain), *_ = res
+
+ def viz(chain, peaks):
+ last_pos = None
+ for pos, pol, nidx in chain:
+ if last_pos:
+ delta = int(round((pos - last_pos) / bit_period))
+ if delta > 1:
+ print(f'skipped {delta-1} symbols at {pos}/{last_pos}')
+
+ # Hole patching routine
+ for i in range(1, delta):
+ est_pos = last_pos + (pos - last_pos) / delta * i
+
+ icandidates = [ (ipos, ipol, iidx) for ipos, ipol, iidx in peaks if abs(est_pos - ipos) < hole_patching_threshold ]
+ if not icandidates:
+ yield None
+ continue
+
+ ipos, ipol, iidx = max(icandidates, key = lambda e: abs(e[1]))
+
+ decoded = iidx*2 + (0 if ipol < 0 else 1)
+ print(f'interpolating, last_pos={last_pos}, delta={delta}, pos={pos}, est={est_pos} dec={decoded}')
+ yield decoded
+
+ decoded = nidx*2 + (0 if pol < 0 else 1)
+ yield decoded
+ if ax:
+ ax.axvline(pos, color='blue', alpha=0.5)
+ ax.text(pos-20, 0.0, f'{decoded}', horizontalalignment='right', verticalalignment='center', color='black')
+
+ last_pos = pos
+
+ decoded = list(viz(chain, peaks_processed))
+ print('decoding [ref|dec]:')
+ match_result = []
+ for shift in range(-ser_maxshift, ser_maxshift):
+ msg = f'=== shift = {shift} ===\n'
+ failures = -shift if shift < 0 else 0 # we're skipping the first $shift symbols
+ a = test_data if shift > 0 else test_data[-shift:]
+ b = decoded if shift < 0 else decoded[shift:]
+ for i, (ref, found) in enumerate(itertools.zip_longest(a, b)):
+ if ref is None: # end of signal
+ break
+ 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 " "} '
+ if ref != found:
+ failures += 1
+ if i%8 == 7:
+ msg += '\n'
+ match_result.append((failures, msg))
+ failures, msg = min(match_result, key=lambda e: e[0])
+ print(msg)
+ ser = failures/len(test_data)
+ print(f'Symbol error rate e={ser}: {failures}/{len(test_data)}')
+ br = sampling_rate / decimation / (2**nbits) * nbits * (1 - ser) * 3600
+ print(f'maximum bitrate r={br} b/h')
+ return ser, br
+
+
+# In[13]:
+
+
+default_params = dict(
+ power_avg_width=2.5,
+ max_lookahead=6.5)
+
+def calculate_ser(reps, v, seed, nbits, thf, duration, decimation):
+ params = dict(default_params)
+ params['signal_amplitude'] = v
+ params['decimation'] = decimation
+ params['nbits'] = nbits
+ params['threshold_factor'] = thf
+
+ out = []
+ for _rep in range(reps):
+ sers, brs = [], []
+ try:
+ ser, br = run_ser_test(**params, sample_duration=duration, print=noprint, seed=seed)
+ out.append((ser, br, None))
+ except Exception as e:
+ #traceback.print_exc()
+ out.append((None, None, f'got {e} seed {seed:08x} params {params}'))
+ return out
+
+import argparse
+parser = argparse.ArgumentParser()
+parser.add_argument('-i', '--index', type=int, default=0)
+parser.add_argument('-b', '--batches', type=int, default=1)
+parser.add_argument('-z', '--chunksize', type=int, default=50)
+parser.add_argument('-r', '--run-id', type=str, default='default')
+group = parser.add_mutually_exclusive_group()
+group.add_argument('-s', '--synthetic-noise', nargs='?', const='grid_freq_psd_spl_108pt.json', type=str)
+group.add_argument('-m', '--measured-noise', nargs='?', default='fmeas_export_ocxo_2day.bin', const='fmeas_export_ocxo_2day.bin', type=str)
+args = parser.parse_args()
+
+# scheduled = [
+# (5, 1.5, 50, 50, 128),
+# (6, 2.0, 50, 50, 128),
+# (6, 2.5, 50, 50, 128),
+# (6, 3.0, 50, 50, 128),
+# (6, 3.5, 50, 50, 128),
+# (6, 4.0, 50, 50, 128),
+# (6, 4.5, 50, 50, 128),
+# (6, 5.0, 50, 50, 128),
+# (6, 5.5, 50, 50, 128),
+# (6, 6.0, 50, 50, 128),
+# (6, 6.5, 50, 50, 128),
+# (6, 7.0, 50, 50, 128),
+# (6, 7.5, 50, 50, 128),
+# (6, 8.0, 50, 50, 128),
+# (6, 8.5, 50, 50, 128),
+# (6, 9.0, 50, 50, 128),
+# (6, 9.5, 50, 50, 128),
+# (6, 10.0, 50, 50, 128),
+# ]
+#scheduled = [ (5, 4.0, 50, 70, 100, dec) for dec in list(range(19)) + [int(round(20 * 10**dec)) for dec in np.linspace(0, 0.7, 20)] ] #int(round(1 * 10**dec))) for dec in np.linspace(0, 2, 30) ]
+
+#param_name = 'par101'
+#scheduled = [ (nbits, thf, 100, 30, duration, 10) for thf in [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0] for nbits, duration in [(5, 128), (6, 128), (7, 64), (8, 64)] ]# int(round(3 * 10**dec)) for dec in np.linspace(0, 2, 30) ]
+
+#param_name = 'par102'
+#scheduled = [ (nbits, thf, 1000, 50, duration, int(round(3 * 10**dec))) for thf in [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0] for nbits, duration in [(5, 128), (6, 128), (7, 64), (8, 64)] for dec in np.linspace(0, 2, 30) ]
+
+#param_name = 'par103'
+#scheduled = [ (nbits, thf, 1000, 50, duration, 10) for thf in [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0] for nbits, duration in [(5, 128), (6, 128), (7, 64), (8, 64)] ]
+
+#param_name = 'par105'
+#scheduled = [ (nbits, thf, 1000, 30, duration, int(round(3 * 10**dec))) for thf in [2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0] for nbits, duration in [(5, 128), (6, 128)] for dec in np.linspace(0, 2, 30) ]
+
+#param_name = 'dbg106'
+#scheduled = [ (nbits, thf, 20, 15, duration, 10) for thf in [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0] for nbits, duration in [(5, 128), (6, 128)] ]
+
+#param_name = 'par107'
+#scheduled = [ (nbits, thf, 100, 30, duration, 10) for thf in [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0] for nbits, duration in [(5, 128), (6, 128), (7, 64), (8, 64)] ]
+
+#param_name = 'par108'
+#scheduled = [ (nbits, thf, 100, 30, duration, int(round(3 * 10**dec))) for thf in [2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0] for nbits, duration in [(5, 128), (6, 128)] for dec in np.linspace(0, 2, 30) ]
+
+#param_name = 'par110'
+#scheduled = [ (nbits, thf, 100, 30, duration, int(round(1 * 10**dec))) for thf in [3.5, 4.0, 4.5, 5.0] for nbits, duration in [(5, 128), (6, 128)] for dec in np.linspace(0.15, 2, 15) ]
+
+#param_name = 'par111'
+#scheduled = [ (nbits, thf, 100, 30, duration, dec) for thf in [3.5, 4.0, 4.5, 5.0] for nbits, duration in [(5, 128), (6, 128)] for dec in [1, 2, 3, 4, 5, 6, 9, 12, 16, 22, 30, 40, 50] ]
+
+# Faster chip duration graph
+#param_name = 'par112'
+#scheduled = [ (nbits, thf, 50, 30, duration, dec) for thf in [4.0] for nbits, duration in [(5, 100)] for dec in [1, 2, 3, 4, 5, 6, 9, 12, 16, 22, 30, 40, 50] ]
+
+#param_name = 'par113'
+#scheduled = [ (nbits, thf, 100, 30, duration, dec) for thf in [4.0, 4.5, 5.0, 5.5] for nbits, duration in [(5, 128), (6, 128)] for dec in [1, 2, 3, 4, 5, 6, 9, 12, 16, 22, 30, 40, 50] ]
+
+#param_name = 'par114'
+#scheduled = [ (nbits, thf, 100, 30, duration, dec) for nbits, thf, duration in [(5, 4.0, 128), (6, 5.0, 128)] for dec in [1, 2, 3, 4, 5, 6, 9, 12, 16, 22, 30, 40, 50] ]
+
+# Extra amplitude for 5 and 6 bit graphs NOTE: accidental duplicate param_name, earlier run is above set, later run this one.
+#param_name = 'par114'
+param_name = 'par115'
+scheduled = [ (nbits, thf, 100, 0.05e-3 * 10 ** np.linspace(0, 3.5, 50), duration, 10) for thf in [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0] for nbits, duration in [(5, 128), (6, 128)] ]
+
+random.Random(0).shuffle(scheduled)
+batch = scheduled[args.index::args.batches]
+print('Parameter set name', param_name)
+
+if args.synthetic_noise is not None:
+ noise_params = load_noise_synth_params(args.synthetic_noise)
+ mains_noise = lambda n: mains_noise_synthetic(n, **noise_params)
+ noise_source = 'synth'
+else:
+ noise_params = load_noise_meas_params(args.measured_noise)
+ mains_noise = lambda n: mains_noise_measured(n, **noise_params)
+ noise_source = 'meas'
+
+def chunked(n, size):
+ for start in range(0, n, size):
+ end = min(reps, start + size)
+ yield start, end
+
+results = {}
+with tqdm(total = 0) as tq:
+ with multiprocessing.Pool(multiprocessing.cpu_count()) as pool:
+ #for nbits, thf, reps, points, duration in [(5, 4.0, 50, 25, 128), (6, 4.0, 50, 25, 128), (7, 5.0, 25, 25, 128), (8, 6.0, 25, 25, 64)]:
+ #for nbits, thf, reps, points, duration in [(6, 4.5, 5, 25, 128), (6, 5.0, 5, 25, 128)]:
+ for nbits, thf, reps, points, duration, decimation in batch:
+ st = np.random.RandomState(0)
+
+ if isinstance(points, numbers.Number):
+ vs = 0.05e-3 * 10 ** np.linspace(0, 2.0, points)
+ points_id = points
+ else:
+ vs = points
+ points_id = len(vs)
+ tq.write('scheduled {}reps x {}pt {}b @ th={} t={}sym dec={}sp'.format(reps, points_id, nbits, thf, duration, decimation))
+ results[(nbits, thf, reps, points_id, duration, decimation)] = { v:
+ [ pool.apply_async(calculate_ser, (end - start, v, st.randint(0xffffffff), nbits, thf, duration, decimation), callback=lambda _res: tq.update(end - start)) for start, end in chunked(reps, args.chunksize) ]
+ for v in vs }
+ tq.total += len(vs) * reps
+ tq.refresh()
+
+ tq.write(f'scheduled {tq.total} tasks. waiting...')
+ pool.close()
+ pool.join()
+
+ results = [ (key, [ (v, [ x for res in chunks for x in res.get() ]) for v, chunks in series.items() ]) for key, series in results.items() ]
+ tq.write('done')
+
+ #sers, brs = np.array(sers), np.array(brs)
+ #ser, std = np.mean(sers), np.std(sers)
+ #print(f'signal_amplitude={v:<.5f}: ser={ser:<.5f} ±{std:<.5f}, br={np.mean(brs):<.5f}')
+ #return ser, std
+
+with open(f'dsss_experiments_res-{param_name}-{noise_source}-{args.run_id}-{args.index}-{datetime.datetime.now():%Y-%m-%d-%H-%M-%S}.json', 'w') as f:
+ json.dump(results, f)
+ print(f'Wrote results to {f.name}')
+
generated by cgit (git 2.34.1) at 2024-11-15 03:25:07 +0000