fw simulator: WIP

This commit is contained in:
jaseg 2020-04-17 17:59:08 +02:00
parent e505627ada
commit 87ae7dfcb3
24 changed files with 465 additions and 192 deletions

View file

@ -3,6 +3,8 @@
# Dependency directories
########################################################################################################################
$(info $(shell env))
CUBE_DIR ?= STM32CubeF4
CMSIS_DIR ?= cmsis
MSPDEBUG_DIR ?= mspdebug
@ -25,13 +27,12 @@ FMEAS_SAMPLING_RATE ?= $(shell echo $(FMEAS_ADC_SAMPLING_RATE) / \($(FMEAS_FFT_
DSSS_GOLD_CODE_NBITS ?= 5
DSSS_DECIMATION ?= 10
# TODO maybe auto adjust this based on detection rate?
DSSS_THESHOLD_FACTOR ?= 5.0f
DSSS_THRESHOLD_FACTOR ?= 5.0f
DSSS_WAVELET_WIDTH ?= 7.3
DSSS_WAVELET_LUT_SIZE ?= 69
DSSS_FILTER_FC ?= 3e-3
DSSS_FILTER_ORDER ?= 12
PAYLOAD_DATA_BIT ?= 64
TRANSMISSION_SYMBOLS ?= 32
PRESIG_STORE_SIZE ?= 3
@ -149,17 +150,16 @@ COMMON_CFLAGS += -DFMEAS_ADC_MAX=$(FMEAS_ADC_MAX)
COMMON_CFLAGS += -DFMEAS_ADC_SAMPLING_RATE=$(FMEAS_ADC_SAMPLING_RATE)
COMMON_CFLAGS += -DFMEAS_FFT_WINDOW_SIGMA=$(FMEAS_FFT_WINDOW_SIGMA)
COMMON_CFLAGS += -DDSSS_DECIMATION=$(DSSS_DECIMATION)
COMMON_CFLAGS += -DDSSS_THESHOLD_FACTOR=$(DSSS_THESHOLD_FACTOR)
COMMON_CFLAGS += -DDSSS_THRESHOLD_FACTOR=$(DSSS_THRESHOLD_FACTOR)
COMMON_CFLAGS += -DDSSS_WAVELET_WIDTH=$(DSSS_WAVELET_WIDTH)
COMMON_CFLAGS += -DDSSS_WAVELET_LUT_SIZE=$(DSSS_WAVELET_LUT_SIZE)
COMMON_CFLAGS += -DPAYLOAD_DATA_BIT=$(PAYLOAD_DATA_BIT)
COMMON_CFLAGS += -DTRANSMISSION_SYMBOLS=$(TRANSMISSION_SYMBOLS)
COMMON_CFLAGS += -DPRESIG_STORE_SIZE=$(PRESIG_STORE_SIZE)
# for musl
CFLAGS += -Dhidden=
SIM_CFLAGS += -lm -DSIMULATION -fsanitize=address
SIM_CFLAGS += -lm -DSIMULATION
SIM_CFLAGS += -Wall -Wextra -Wpedantic -Wshadow -Wimplicit-function-declaration -Wundef -Wno-unused-parameter
INT_CFLAGS += -Wall -Wextra -Wpedantic -Wshadow -Wimplicit-function-declaration -Wundef -Wno-unused-parameter
@ -203,7 +203,8 @@ binsize: $(BUILDDIR)/$(BINARY) $(BUILDDIR)/$(BINARY:.elf=-symbol-sizes.pdf)
@echo "▐▬▬▬▌ SyMbOL sIzE HiGhScORe LiSt ▐▬▬▬▌"
$(NM) --print-size --size-sort --radix=d $< | tail -n 20
src/dsss_demod.c: $(BUILDDIR)/generated/dsss_gold_code.h $(BUILDDIR)/generated/dsss_butter_filter.h
# $(BUILDDIR)/generated/dsss_butter_filter.h
src/dsss_demod.c: $(BUILDDIR)/generated/dsss_gold_code.h
$(BUILDDIR)/generated/dsss_gold_code.h: $(BUILDDIR)/generated/gold_code_$(DSSS_GOLD_CODE_NBITS).h
ln -srf $< $@

View file

@ -13,23 +13,23 @@
#include "simulation.h"
#include "generated/dsss_gold_code.h"
#include "generated/dsss_butter_filter.h"
// #include "generated/dsss_butter_filter.h"
/* Generated CWT wavelet LUT */
extern const float * const dsss_cwt_wavelet_table;
struct iir_biquad cwt_filter_bq[DSSS_FILTER_CLEN] = {DSSS_FILTER_COEFF};
//struct iir_biquad cwt_filter_bq[DSSS_FILTER_CLEN] = {DSSS_FILTER_COEFF};
void debug_print_vector(const char *name, size_t len, const float *data, size_t stride, bool index, bool debug);
static float gold_correlate_step(const size_t ncode, const float a[DSSS_CORRELATION_LENGTH], size_t offx, bool debug);
static float cwt_convolve_step(const float v[DSSS_WAVELET_LUT_SIZE], size_t offx);
static float run_iir(const float x, const int order, const struct iir_biquad q[order], struct iir_biquad_state st[order]);
static float run_biquad(float x, const struct iir_biquad *const q, struct iir_biquad_state *const restrict st);
//static float run_iir(const float x, const int order, const struct iir_biquad q[order], struct iir_biquad_state st[order]);
//static float run_biquad(float x, const struct iir_biquad *const q, struct iir_biquad_state *const restrict st);
static void matcher_init(struct matcher_state states[static DSSS_MATCHER_CACHE_SIZE]);
static void matcher_tick(struct matcher_state states[static DSSS_MATCHER_CACHE_SIZE],
uint64_t ts, int peak_ch, float peak_ampl);
static void group_received(struct dsss_demod_state *st);
static uint8_t decode_peak(int peak_ch, float peak_ampl);
static symbol_t decode_peak(int peak_ch, float peak_ampl);
#ifdef SIMULATION
void debug_print_vector(const char *name, size_t len, const float *data, size_t stride, bool index, bool debug) {
@ -61,7 +61,7 @@ void dsss_demod_init(struct dsss_demod_state *st) {
void dsss_demod_step(struct dsss_demod_state *st, float new_value, uint64_t ts) {
//const float hole_patching_threshold = 0.01 * DSSS_CORRELATION_LENGTH;
bool log = false;
bool log_groups = true;
bool log_groups = false;
st->signal[st->signal_wpos] = new_value;
st->signal_wpos = (st->signal_wpos + 1) % ARRAY_LENGTH(st->signal);
@ -99,7 +99,7 @@ void dsss_demod_step(struct dsss_demod_state *st, float new_value, uint64_t ts)
max_ts = ts;
}
if (fabsf(val) > DSSS_THESHOLD_FACTOR)
if (fabsf(val) > DSSS_THRESHOLD_FACTOR)
found = true;
}
if (log) DEBUG_PRINTN("%f %d ", max_val, found);
@ -134,12 +134,12 @@ void dsss_demod_step(struct dsss_demod_state *st, float new_value, uint64_t ts)
/* Map a sequence match to a data symbol. This maps the sequence's index number to the 2nd to n+2nd bit of the result,
* and maps the polarity of detection to the LSb. 5-bit example:
*
* [0, S, S, S, S, S, S, P] ; S ^= symbol index (0 - 2^n+1), P ^= symbol polarity
* [0, S, S, S, S, S, S, P] ; S ^= symbol index (0 - 2^n+1 so we need just about n+1 bit), P ^= symbol polarity
*
* Symbol polarity is preserved from transmitter to receiver. The symbol index is n+1 bit instead of n bit since we have
* 2^n+1 symbols to express, one too many for an n-bit index.
*/
uint8_t decode_peak(int peak_ch, float peak_ampl) {
symbol_t decode_peak(int peak_ch, float peak_ampl) {
return (peak_ch<<1) | (peak_ampl > 0);
}
@ -157,7 +157,7 @@ void matcher_tick(struct matcher_state states[static DSSS_MATCHER_CACHE_SIZE], u
const float score_depreciation = 0.1f; /* 0.0 -> no depreciation, 1.0 -> complete disregard */
const int current_phase = ts % DSSS_CORRELATION_LENGTH;
const int max_skips = TRANSMISSION_SYMBOLS/4*3;
bool debug = true;
bool debug = false;
bool header_printed = false;
for (size_t i=0; i<DSSS_MATCHER_CACHE_SIZE; i++) {
@ -172,7 +172,7 @@ void matcher_tick(struct matcher_state states[static DSSS_MATCHER_CACHE_SIZE], u
header_printed = true;
DEBUG_PRINTN("windows %zu\n", ts);
}
if (debug) DEBUG_PRINTN(" skip %d old=%f new=%f\n", i, states[i].candidate_score, score);
if (debug) DEBUG_PRINTN(" skip %zd old=%f new=%f\n", i, states[i].candidate_score, score);
/* We win, update candidate */
assert(i < DSSS_MATCHER_CACHE_SIZE);
states[i].candidate_score = score;
@ -195,7 +195,7 @@ void matcher_tick(struct matcher_state states[static DSSS_MATCHER_CACHE_SIZE], u
header_printed = true;
DEBUG_PRINTN("windows %zu\n", ts);
}
if (debug) DEBUG_PRINTN(" %d ", i);
if (debug) DEBUG_PRINTN(" %zd ", i);
/* Process window results */
assert(i < DSSS_MATCHER_CACHE_SIZE);
assert(0 <= states[i].data_pos && states[i].data_pos < TRANSMISSION_SYMBOLS);
@ -235,7 +235,7 @@ static float score_group(const struct group *g, int phase_delta) {
}
void group_received(struct dsss_demod_state *st) {
bool debug = true;
bool debug = false;
const int group_phase = st->group.max_ts % DSSS_CORRELATION_LENGTH;
/* This is the score of a decoding starting at this group (with no context) */
float base_score = score_group(&st->group, 0);
@ -305,6 +305,7 @@ void group_received(struct dsss_demod_state *st) {
}
}
#if 0
float run_iir(const float x, const int order, const struct iir_biquad q[order], struct iir_biquad_state st[order]) {
float intermediate = x;
for (int i=0; i<(order+1)/2; i++)
@ -320,6 +321,7 @@ float run_biquad(float x, const struct iir_biquad *const q, struct iir_biquad_st
st->reg[0] = intermediate;
return out;
}
#endif
float cwt_convolve_step(const float v[DSSS_WAVELET_LUT_SIZE], size_t offx) {
float sum = 0.0f;

View file

@ -10,8 +10,12 @@
/* FIXME: move to makefile */
#define DSSS_MATCHER_CACHE_SIZE 8
/* FIXME: move to more appropriate header */
#define PAYLOAD_DATA_BYTE ((PAYLOAD_DATA_BIT+7)/8)
#if DSSS_GOLD_CODE_NBITS < 8
typedef uint8_t symbol_t;
#else
typedef uint16_t symbol_t;
#endif
struct iir_biquad {
float a[2];
@ -43,12 +47,9 @@ struct matcher_state {
int last_skips;
int candidate_skips;
#if DSSS_GOLD_CODE_NBITS > 7
#error DSSS_GOLD_CODE_NBITS is too large for matcher_state.data data type (uint8_t)
#endif
uint8_t data[TRANSMISSION_SYMBOLS];
symbol_t data[TRANSMISSION_SYMBOLS];
int data_pos;
uint8_t candidate_data;
symbol_t candidate_data;
};
struct dsss_demod_state {
@ -66,7 +67,7 @@ struct dsss_demod_state {
};
extern void handle_dsss_received(uint8_t data[static TRANSMISSION_SYMBOLS]);
extern void handle_dsss_received(symbol_t data[static TRANSMISSION_SYMBOLS]);
void dsss_demod_init(struct dsss_demod_state *st);
void dsss_demod_step(struct dsss_demod_state *st, float new_value, uint64_t ts);

View file

@ -3,7 +3,7 @@
#ifdef SIMULATION
#include <stdio.h>
#define DEBUG_PRINTN(...) fprintf(stderr, __VA_ARGS__)
#define DEBUG_PRINTN(...) printf(__VA_ARGS__)
#define DEBUG_PRINTNF(fmt, ...) DEBUG_PRINTN("%s:%d: " fmt, __FILE__, __LINE__, ##__VA_ARGS__)
#define DEBUG_PRINT(fmt, ...) DEBUG_PRINTNF(fmt "\n", ##__VA_ARGS__)
#else

View file

@ -18,7 +18,7 @@ if __name__ == '__main__':
print(f'const float {varname}[{args.n}] = {{')
win = sig.ricker(args.n, args.w)
par = ' '.join(f'{f:>015.8g}f,' for f in win)
par = ' '.join(f'{f:>015.12e}f,' for f in win)
print(textwrap.fill(par,
initial_indent=' '*4, subsequent_indent=' '*4,
width=120,

View file

@ -12,7 +12,7 @@
#include "dsss_demod.h"
void handle_dsss_received(uint8_t data[static TRANSMISSION_SYMBOLS]) {
void handle_dsss_received(symbol_t data[static TRANSMISSION_SYMBOLS]) {
printf("data sequence received: [ ");
for (size_t i=0; i<TRANSMISSION_SYMBOLS; i++) {
printf("%+3d", ((data[i]&1) ? 1 : -1) * (data[i]>>1));

View file

@ -4,36 +4,204 @@ import os
from os import path
import subprocess
import json
from collections import namedtuple, defaultdict
from tqdm import tqdm
import uuid
import multiprocessing
import sqlite3
import time
from urllib.parse import urlparse
import functools
import tempfile
import itertools
import numpy as np
np.set_printoptions(linewidth=240)
from dsss_demod_test_waveform_gen import load_noise_meas_params, load_noise_synth_params,\
mains_noise_measured, mains_noise_synthetic, modulate as dsss_modulate
def build_test_binary(nbits, thf, decimation, symbols, cachedir):
build_id = str(uuid.uuid4())
builddir = path.join(cachedir, build_id)
os.mkdir(builddir)
cwd = path.join(path.dirname(__file__), '..')
env = os.environ.copy()
env['BUILDDIR'] = path.abspath(builddir)
env['DSSS_GOLD_CODE_NBITS'] = str(nbits)
env['DSSS_DECIMATION'] = str(decimation)
env['DSSS_THRESHOLD_FACTOR'] = str(thf)
env['DSSS_WAVELET_WIDTH'] = str(0.73 * decimation)
env['DSSS_WAVELET_LUT_SIZE'] = str(10 * decimation)
env['TRANSMISSION_SYMBOLS'] = str(symbols)
with open(path.join(builddir, 'make_stdout.txt'), 'w') as stdout,\
open(path.join(builddir, 'make_stderr.txt'), 'w') as stderr:
subprocess.run(['make', 'clean', os.path.abspath(path.join(builddir, 'tools/dsss_demod_test'))],
env=env, cwd=cwd, check=True, stdout=stdout, stderr=stderr)
return build_id
@functools.lru_cache()
def load_noise_gen(url):
schema, refpath = url.split('://')
if not path.isabs(refpath):
refpath = path.abspath(path.join(path.dirname(__file__), refpath))
if schema == 'meas':
return mains_noise_measured, load_noise_meas_params(refpath)
elif schema == 'synth':
return mains_noise_synthetic, load_noise_synth_params(refpath)
else:
raise ValueError('Invalid schema', schema)
def sequence_matcher(test_data, decoded, max_shift=3):
match_result = []
for shift in range(-max_shift, max_shift):
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
if ref != found:
failures += 1
match_result.append(failures)
failures = min(match_result)
return failures/len(test_data)
ResultParams = namedtuple('ResultParams', ['nbits', 'thf', 'decimation', 'symbols', 'seed', 'amplitude', 'background'])
def run_test(seed, amplitude_spec, background, nbits, decimation, symbols, thfs, lookup_binary, cachedir):
noise_gen, noise_params = load_noise_gen(background)
test_data = np.random.RandomState(seed=seed).randint(0, 2 * (2**nbits), symbols)
signal = np.repeat(dsss_modulate(test_data, nbits) * 2.0 - 1, decimation)
# We're re-using the seed here. This is not a problem.
noise = noise_gen(seed, len(signal), *noise_params)
amplitudes = amplitude_spec[0] * 10 ** np.linspace(0, amplitude_spec[1], amplitude_spec[2])
output = []
for amp in amplitudes:
with tempfile.NamedTemporaryFile(dir=cachedir) as f:
waveform = signal*amp + noise
f.write(waveform.astype('float').tobytes())
f.flush()
for thf in thfs:
cmdline = [lookup_binary(nbits, thf, decimation, symbols), f.name]
proc = subprocess.Popen(cmdline, stdout=subprocess.PIPE, text=True)
stdout, _stderr = proc.communicate()
if proc.returncode != 0:
raise SystemError(f'Subprocess signalled error: {proc.returncode=}')
lines = stdout.splitlines()
matched = [ l.partition('[')[2].partition(']')[0]
for l in lines if l.strip().startswith('data sequence received:') ]
matched = [ [ int(elem) for elem in l.split(',') ] for l in matched ]
ser = min(sequence_matcher(test_data, match) for match in matched) if matched else None
rpars = ResultParams(nbits, thf, decimation, symbols, seed, amp, background)
output.append((rpars, ser))
print(f'ran {rpars} {ser=} {" ".join(cmdline)}')
return output
def parallel_generator(db, table, columns, builder, param_list, desc, context={}, params_mapper=lambda *args: args):
with multiprocessing.Pool(multiprocessing.cpu_count()) as pool:
with db as conn:
jobs = []
for params in param_list:
found_res = conn.execute(
f'SELECT result FROM {table} WHERE ({",".join(columns)}) = ({",".join("?"*len(columns))})',
params_mapper(*params)).fetchone()
if found_res:
yield params, json.loads(*found_res)
else:
jobs.append((params, pool.apply_async(builder, params, context)))
pool.close()
print('Using', len(param_list) - len(jobs), 'cached jobs', flush=True)
with tqdm(total=len(jobs), desc=desc) as tq:
for params, res in jobs:
tq.update(1)
result = res.get()
with db as conn:
conn.execute(f'INSERT INTO {table} VALUES ({"?,"*len(params)}?,?)',
(*params_mapper(*params), json.dumps(result), timestamp()))
yield params, result
pool.join()
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser()
parser.add_argument(metavar='test_data_directory', dest='dir', help='Directory with test data .bin files')
default_binary = path.abspath(path.join(path.dirname(__file__), '../build/tools/dsss_demod_test'))
parser.add_argument(metavar='test_binary', dest='binary', nargs='?', default=default_binary)
parser.add_argument('-d', '--dump', help='Write raw measurements to JSON file')
parser.add_argument('-d', '--dump', help='Write results to JSON file')
parser.add_argument('-c', '--cachedir', default='dsss_test_cache', help='Directory to store build output and data in')
args = parser.parse_args()
bin_files = [ path.join(args.dir, d) for d in os.listdir(args.dir) if d.lower().endswith('.bin') ]
DecoderParams = namedtuple('DecoderParams', ['nbits', 'thf', 'decimation', 'symbols'])
dec_paramses = [ DecoderParams(nbits=nbits, thf=thf, decimation=decimation, symbols=20)
for nbits in [5, 6]
for thf in [4.5, 4.0, 5.0]
for decimation in [10, 5, 22] ]
# dec_paramses = [ DecoderParams(nbits=nbits, thf=thf, decimation=decimation, symbols=100)
# for nbits in [5, 6, 7, 8]
# 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 decimation in [1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 16, 22, 30, 40, 50] ]
savedata = {}
for p in bin_files:
output = subprocess.check_output([args.binary, p], stderr=subprocess.DEVNULL)
measurements = np.array([ float(value) for _offset, value in [ line.split() for line in output.splitlines() ] ])
savedata[p] = list(measurements)
build_cache_dir = path.join(args.cachedir, 'builds')
data_cache_dir = path.join(args.cachedir, 'data')
os.makedirs(build_cache_dir, exist_ok=True)
os.makedirs(data_cache_dir, exist_ok=True)
# Cut off first and last sample for mean and RMS calculations as these show boundary effects.
measurements = measurements[1:-1]
mean = np.mean(measurements)
rms = np.sqrt(np.mean(np.square(measurements - mean)))
build_db = sqlite3.connect(path.join(args.cachedir, 'build_db.sqlite3'))
build_db.execute('CREATE TABLE IF NOT EXISTS builds (nbits, thf, decimation, symbols, result, timestamp)')
timestamp = lambda: int(time.time()*1000)
print(f'{path.basename(p):<60}: mean={mean:<8.4f}Hz rms={rms*1000:.3f}mHz')
builds = dict(parallel_generator(build_db, table='builds', columns=['nbits', 'thf', 'decimation', 'symbols'],
builder=build_test_binary, param_list=dec_paramses, desc='Building decoders',
context=dict(cachedir=build_cache_dir)))
print('Done building decoders.')
GeneratorParams = namedtuple('GeneratorParams', ['seed', 'amplitude_spec', 'background'])
gen_params = [ GeneratorParams(rep, (5e-3, 1, 5), background)
#GeneratorParams(rep, (0.05e-3, 3.5, 50), background)
for rep in range(30)
for background in ['meas://fmeas_export_ocxo_2day.bin', 'synth://grid_freq_psd_spl_108pt.json'] ]
data_db = sqlite3.connect(path.join(args.cachedir, 'data_db.sqlite3'))
data_db.execute('CREATE TABLE IF NOT EXISTS waveforms'
'(seed, amplitude_spec, background, nbits, decimation, symbols, thresholds, result, timestamp)')
dec_param_groups = defaultdict(lambda: [])
for nbits, thf, decimation, symbols in dec_paramses:
dec_param_groups[(nbits, decimation, symbols)].append(thf)
waveform_params = [ (*gp, *dp, thfs) for gp in gen_params for dp, thfs in dec_param_groups.items() ]
print(f'Generated {len(waveform_params)} parameter sets')
def lookup_binary(*params):
return path.join(build_cache_dir, builds[tuple(params)], 'tools/dsss_demod_test')
def params_mapper(seed, amplitude_spec, background, nbits, decimation, symbols, thresholds):
amplitude_spec = ','.join(str(x) for x in amplitude_spec)
thresholds = ','.join(str(x) for x in thresholds)
return seed, amplitude_spec, background, nbits, decimation, symbols, thresholds
results = []
for _params, chunk in parallel_generator(data_db, 'waveforms',
['seed', 'amplitude_spec', 'background', 'nbits', 'decimation', 'symbols', 'thresholds'],
params_mapper=params_mapper,
builder=run_test,
param_list=waveform_params, desc='Generating waveforms',
context=dict(cachedir=data_cache_dir, lookup_binary=lookup_binary)):
results += chunk
if args.dump:
with open(args.dump, 'w') as f:
json.dump(savedata, f)
json.dump(results, f)

View file

@ -0,0 +1,71 @@
import functools
import numpy as np
import numbers
import math
from scipy import signal as sig
import scipy.fftpack
sampling_rate = 10 # sp/s
# 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)
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)) ])
def load_noise_meas_params(capture_file):
with open(capture_file, 'rb') as f:
meas_data = np.copy(np.frombuffer(f.read(), dtype='float32'))
meas_data -= np.mean(meas_data)
return (meas_data,)
def mains_noise_measured(seed, n, meas_data):
last_valid = len(meas_data) - n
st = np.random.RandomState(seed)
start = st.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 (np.linspace(*d['x_spec']), # spl_x
d['x_spec'][2], # spl_N
(d['t'], d['c'], d['k'])) # psd_spl
def mains_noise_synthetic(seed, n, psd_spl, spl_N, spl_x):
st = np.random.RandomState(seed)
noise = st.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]

Binary file not shown.

View file

@ -0,0 +1 @@
{"x_spec": [3.2595692805152726e-05, 5.0, 613575], "t": [3.2595692805152726e-05, 3.2595692805152726e-05, 3.2595692805152726e-05, 3.2595692805152726e-05, 0.0001423024947075771, 0.00015800362803968106, 0.00017543716661470822, 0.00019479425764873777, 0.0002162871388378975, 0.00024015146540428407, 0.00026664889389955537, 0.00029606995109590574, 0.00032873721941990017, 0.0003650088738553592, 0.0004052826090950758, 0.00045000000000000004, 0.000499651343175437, 0.0005547810327489297, 0.0006159935292916862, 0.0006839599873288199, 0.0007594256141046668, 0.0008432178402871724, 0.0009362553921977272, 0.0010395583650374223, 0.0011542594075560205, 0.001281616140796111, 0.0014230249470757708, 0.001580036280396809, 0.0017543716661470824, 0.0019479425764873776, 0.002162871388378975, 0.0024015146540428403, 0.002666488938995554, 0.002960699510959057, 0.0032873721941990056, 0.0036500887385535925, 0.004052826090950754, 0.0045000000000000005, 0.00499651343175437, 0.005547810327489296, 0.006159935292916869, 0.0068395998732882, 0.007594256141046669, 0.008432178402871724, 0.009362553921977271, 0.010395583650374221, 0.011542594075560205, 0.012816161407961109, 0.014230249470757707, 0.01580036280396809, 0.017543716661470823, 0.01947942576487376, 0.02162871388378975, 0.024015146540428405, 0.026664889389955565, 0.02960699510959057, 0.03287372194199005, 0.036500887385535925, 0.04052826090950754, 0.045, 0.0499651343175437, 0.05547810327489296, 0.06159935292916863, 0.06839599873288206, 0.07594256141046668, 0.08432178402871732, 0.09362553921977272, 0.10395583650374222, 0.11542594075560206, 0.12816161407961107, 0.14230249470757705, 0.15800362803968088, 0.1754371666147082, 0.1947942576487376, 0.21628713883789774, 0.24015146540428406, 0.26664889389955565, 0.2960699510959057, 0.32873721941990053, 0.36500887385535924, 0.40528260909507535, 0.45, 0.499651343175437, 0.5547810327489296, 0.6159935292916868, 0.6839599873288206, 0.7594256141046669, 0.8432178402871732, 0.9362553921977271, 1.0395583650374223, 1.1542594075560206, 1.2816161407961109, 1.4230249470757708, 1.5800362803968104, 1.7543716661470823, 1.9479425764873777, 2.162871388378975, 2.4015146540428405, 2.6664889389955535, 2.960699510959057, 3.287372194199002, 3.6500887385535927, 4.052826090950758, 4.5, 5.0, 5.0, 5.0, 5.0], "c": [0.7720161468716866, -0.5547528253056444, 0.30706059086000753, 0.19422577014134906, -1.1954636661840032, 0.9215976941641111, -0.6668136393976918, -1.341269161156733, -0.16311330594842666, -1.7639636752234251, -1.238385544822954, -0.32649555618555554, -0.03086589610280171, -2.358195657381619, -0.5759152419849985, 0.1892225800004134, -1.8122889670546236, -0.8109120798216202, -0.5500991736738969, -4.680192969256771, -2.8007700704649876, 0.16866469558571784, -1.1040811840849307, -3.0243574268705546, -4.018139927365795, -4.100581028618109, -0.556354762846191, -7.414377514669229, 1.36396325920194, -6.002559557058508, -2.2113451390305365, -4.578944771104116, -4.372644849632638, -3.945339124673235, -4.778747958903158, -2.370174137632325, -5.7372466088109295, -4.707506574819875, -4.834404729330929, -5.005244244061701, -5.82644896783577, -4.717966026411524, -6.146374820241562, -4.972788381244952, -5.854957092953355, -5.702174935205885, -6.222035857079607, -6.2128389666872, -6.212821706753751, -6.253599689326325, -6.681685577659057, -6.372364384360678, -6.771223202540934, -6.856809137231159, -6.986412256164045, -7.190466178818742, -7.577896455149433, -7.515731696006047, -7.598155006351761, -7.824526916149126, -8.141496591776512, -8.36794927682997, -8.80307396767114, -8.828816533544659, -9.357524260470413, -9.658130054343863, -10.005768472049466, -10.499801262514108, -11.028689820560558, -11.413688641742898, -11.906162042727946, -12.232342460719975, -12.438432746733596, -13.088338100203112, -12.308710772618745, -11.685074853925329, -11.397838681243094, -12.265219694936695, -13.600359694898529, -14.031425961884718, -12.236885080485473, -13.527508426900974, -13.698402018452601, -13.397911198962568, -14.144410560196603, -13.905769594095293, -14.410874830544122, -14.531727635304264, -14.59275291853806, -14.35404826562502, -14.58670053318149, -14.432515268864977, -14.363428024828353, -14.429222027493264, -14.73947634127499, -14.717315405960353, -14.678539669792505, -14.825278423641382, -14.80936417940876, -14.943375264882789, -14.680885181815674, -14.54841244844906, -14.634365225950589, -14.609444790868906, 0.0, 0.0, 0.0, 0.0], "k": 3}

File diff suppressed because one or more lines are too long

View file

@ -701,7 +701,7 @@
" ax.set_xlabel('Period T [s]')\n",
" ax.set_ylabel('Power Spectral Density [Hz^2/Hz]')\n",
"\n",
" for i, t in enumerate([2.0, 3.2, 30, 150, 220, 450, 600, 900]):\n",
" for i, t in enumerate([2.0, 3.2, 30, 150, 220, 450, 600, 900]):z\n",
" ax.axvline(1/t, color='red', alpha=0.5, zorder=-1)\n",
" ax.annotate(f'{t} s', xy=(1/t, 1e-7), xytext=(-10, 5), xycoords='data', textcoords='offset pixels', rotation=90)\n",
" #ax.text(1/60, 10,'60 s', ha='left')\n",

View file

@ -12,6 +12,7 @@ all: safety_reset.pdf
safety_reset.pdf: resources/grid_freq_estimation.pdf
safety_reset.pdf: resources/gps_clock_jitter_analysis.pdf
safety_reset.pdf: resources/dsss_experiments-ber.pdf
%.pdf: %.tex %.bib
pdflatex -shell-escape $<

View file

@ -1183,16 +1183,30 @@ indicates SER is related fairly monotonically to the signal-to-noise margins ins
\begin{figure}
\centering
\includegraphics[width=\textwidth]{../lab-windows/fig_out/dsss_gold_nbits_overview}
\includegraphics{../lab-windows/fig_out/dsss_gold_nbits_overview}
\caption{
Symbol Error Rate (SER) as a function of transmission amplitude. The line indicates the mean of several
measurements for each parameter set. The shaded areas indicate one standard deviation from the mean. Background
noise for each trial is a random segment of measured grid frequency. Background noise amplitude is the same for
all trials. Shown are four traces for four different DSSS sequence lengths. Using a 5-bit gold code, one DSSS
symbol measures 31 chips. 6 bit per symbol are 63 chips, 7 bit are 127 chips and 8 bit 255 chips. This
simulation uses a decimation of 10, which corresponds to an $1 \text{s}$ chip length at our $10 \text{Hz}$ grid
frequency sampling rate. At 5 bit per symbol, one symbol takes $31 \text{s}$ and one bit takes $6.2 \text{s}$
amortized. At 8 bit one symbol takes $255 \text{s} = 4 \text{min} 15 \text{s}$ and one bit takes $31.9 \text{s}$
amortized. Here, slower transmission speed buys coding gain. All else being the same this allows for a decrease
in transmission power.
}
\label{dsss_gold_nbits_overview}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\textwidth]{../lab-windows/fig_out/dsss_gold_nbits_sensitivity}
\includegraphics{../lab-windows/fig_out/dsss_gold_nbits_sensitivity}
\caption{
Amplitude at a SER of 0.5\ in mHz depending on symbol length. Here we can observe an increase of sensitivity
with increasing symbol length, but we can clearly see diminishing returns above 6 bit (63 chips). Considering
that each bit roughly doubles overall transmission time for a given data length it seems lower bit counts are
preferrable if the necessary transmitter power can be realized.
}
\label{dsss_gold_nbits_sensitivity}
\end{figure}
@ -1200,20 +1214,38 @@ indicates SER is related fairly monotonically to the signal-to-noise margins ins
\begin{figure}
\begin{subfigure}{\textwidth}
\centering
\includegraphics[width=\textwidth]{../lab-windows/fig_out/dsss_thf_amplitude_5678}
\includegraphics{../lab-windows/fig_out/dsss_thf_amplitude_5678}
\label{dsss_thf_amplitude_5678}
\caption{
\footnotesize SER vs.\ amplitude graph similar to fig.\ \ref{dsss_gold_nbits_overview} with dependence on
threshold factor color-coded. Each graph shows traces for a single DSSS symbol length.
}
\end{subfigure}
\end{figure}
\begin{figure}
\ContinuedFloat
\begin{subfigure}{\textwidth}
\centering
\includegraphics[width=\textwidth]{../lab-windows/fig_out/dsss_thf_sensitivity_5678}
\includegraphics{../lab-windows/fig_out/dsss_thf_sensitivity_5678}
\label{dsss_thf_sensitivity_5678}
\caption{
\footnotesize Graphs of amplitude at $SER=0.5$ for each symbol length as well as asymptotic SER for large
amplitudes. Areas shaded red indicate that $SER=0.5$ was not reached for any amplitude in the simulated
range. We can observe that smaller symbol lengths favor lower threshold factors, and that optimal threshold
factors for all symbol lengths are between $4.0$ and $5.0$.
}
\end{subfigure}
\caption{
}
Dependence of demodulator sensitivity on the threshold factor used for correlation peak detection in our
DSSS demodulator. This is an empirically-determined parameter specific to our demodulation algorithm. At low
threshold factors our classifier yields lots of spurious peaks that have to be thrown out by our maximum
likelihood estimator. These spurious peaks have a random time distribution and thus do not pose much of a
challenge to our MLE but at very low threshold factors the number of spurious peaks slows down decoding and
does still clog our MLE's internal size-limited candidate lists which leads to failed decodings. At very
high threshold factors decoding performance suffers greatly since many valid correlation peaks get
incorrectly ignored. The glitches at medium threshold factors in the 7- and 8-bit graphs are artifacts of
our prototype decoding algorithm that we have not fixed in the prototype implementation since we wanted to
focus on the final C version.}
\label{dsss_thf_sensitivity}
\end{figure}
@ -1223,16 +1255,31 @@ indicates SER is related fairly monotonically to the signal-to-noise margins ins
\includegraphics[width=\textwidth]{../lab-windows/fig_out/chip_duration_sensitivity_5}
\label{chip_duration_sensitivity_5}
\caption{
5 bit Gold code
}
\end{subfigure}
\end{figure}
\begin{figure}
\ContinuedFloat
\begin{subfigure}{\textwidth}
\centering
\includegraphics[width=\textwidth]{../lab-windows/fig_out/chip_duration_sensitivity_6}
\label{chip_duration_sensitivity_6}
\caption{
6 bit Gold code
}
\end{subfigure}
\caption{
Dependence of demodulator sensitivity on DSSS chip duration. Due to computational constraints this simulation is
limited to 5 bit and 6 bit DSSS sequences. There is a clearly visible sensitivity maximum at fairly short chip
lengths around $0.2 \text{s}$. Short chip durations shift the entire transmission band up in frequency. In fig.\
\ref{freq_meas_spectrum} we can see that noise energy is mostly concentrated at lower frequencies, so shifting
our signal up in frequency will reduce the amount of noise the decoder sees behind the correlator by shifting
the band of interest into a lower-noise spectral region. For a practical implementation chip duration is limited
by physical factors such as the maximum modulation slew rate ($\frac{\text{d}P}{\text{d}t}$), the maximum
Rate-Of-Change-Of-Frequency (ROCOF, $\frac{\text{d}f}{\text{d}t}$) the grid can tolerate and possible inertial
effects limiting response of frequency to load changes at certain load levels.
% FIXME are these inertial effects likely? Ask an expert.
}
\label{chip_duration_sensitivity}
\end{figure}
@ -1243,16 +1290,25 @@ indicates SER is related fairly monotonically to the signal-to-noise margins ins
\includegraphics[width=\textwidth]{../lab-windows/fig_out/chip_duration_sensitivity_cmp_meas_6}
\label{chip_duration_sensitivity_cmp_meas_6}
\caption{
Simulation using baseline frequency data from actual measurements.
}
\end{subfigure}
\end{figure}
\begin{figure}
\ContinuedFloat
\begin{subfigure}{\textwidth}
\centering
\includegraphics[width=\textwidth]{../lab-windows/fig_out/chip_duration_sensitivity_cmp_synth_6}
\label{chip_duration_sensitivity_cmp_synth_6}
\caption{
Simulation using synthetic frequency data.
}
\end{subfigure}
\caption{
Chip duration/sensitivity simulation results like in fig.\ \ref{chip_duration_sensitivity} compared between a
simulation using measured frequency data like previous graphs and one using artificially generated noise. There
is almost no visible difference indicating that we have found a good model of reality in our noise synthesizer,
but also that real grid frequency behaves like a frequency-shaped gaussian noise process.
}
\label{chip_duration_sensitivity_cmp}
\end{figure}
@ -1355,6 +1411,7 @@ correctly configure than it is to simply use separate hardware and secure the in
\includenotebook{Grid frequency estimation}{grid_freq_estimation}
\includenotebook{Frequency sensor clock stability analysis}{gps_clock_jitter_analysis}
\includenotebook{DSSS modulation experiments}{dsss_experiments-ber}
\chapter{Demonstrator schematics and code}