Add simulation script

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}')
+