diff options
-rw-r--r-- | controller/fw/Makefile | 4 | ||||
-rw-r--r-- | controller/fw/src/dsss_demod.c | 182 | ||||
-rw-r--r-- | controller/fw/src/dsss_demod.h | 23 | ||||
-rw-r--r-- | controller/fw/tools/dsss_demod_test.c | 12 | ||||
-rw-r--r-- | lab-windows/dsss_experiments-ber.ipynb | 126 |
5 files changed, 288 insertions, 59 deletions
diff --git a/controller/fw/Makefile b/controller/fw/Makefile index c549e47..fac9169 100644 --- a/controller/fw/Makefile +++ b/controller/fw/Makefile @@ -57,9 +57,9 @@ DSSS_WAVELET_WIDTH ?= 7.3 DSSS_WAVELET_LUT_SIZE ?= 69 DSSS_FILTER_FC ?= 3e-3 DSSS_FILTER_ORDER ?= 12 -DSSS_GROUP_CACHE_SIZE ?= 12 PAYLOAD_DATA_BIT ?= 64 +TRANSMISSION_SYMBOLS ?= 32 CC := $(PREFIX)gcc CXX := $(PREFIX)g++ @@ -105,8 +105,8 @@ COMMON_CFLAGS += -DDSSS_DECIMATION=$(DSSS_DECIMATION) COMMON_CFLAGS += -DDSSS_THESHOLD_FACTOR=$(DSSS_THESHOLD_FACTOR) COMMON_CFLAGS += -DDSSS_WAVELET_WIDTH=$(DSSS_WAVELET_WIDTH) COMMON_CFLAGS += -DDSSS_WAVELET_LUT_SIZE=$(DSSS_WAVELET_LUT_SIZE) -COMMON_CFLAGS += -DDSSS_GROUP_CACHE_SIZE=$(DSSS_GROUP_CACHE_SIZE) COMMON_CFLAGS += -DPAYLOAD_DATA_BIT=$(PAYLOAD_DATA_BIT) +COMMON_CFLAGS += -DTRANSMISSION_SYMBOLS=$(TRANSMISSION_SYMBOLS) # for musl CFLAGS += -Dhidden= diff --git a/controller/fw/src/dsss_demod.c b/controller/fw/src/dsss_demod.c index 1cad6f8..32d6104 100644 --- a/controller/fw/src/dsss_demod.c +++ b/controller/fw/src/dsss_demod.c @@ -2,6 +2,8 @@ #include <unistd.h> #include <stdbool.h> #include <math.h> +#include <stdlib.h> +#include <assert.h> #include <arm_math.h> @@ -18,11 +20,15 @@ extern const float * const dsss_cwt_wavelet_table; struct iir_biquad cwt_filter_bq[DSSS_FILTER_CLEN] = {DSSS_FILTER_COEFF}; -float gold_correlate_step(const size_t ncode, const float a[DSSS_CORRELATION_LENGTH], size_t offx, bool debug); -float cwt_convolve_step(const float v[DSSS_WAVELET_LUT_SIZE], size_t offx); -float run_iir(const float x, const int order, const struct iir_biquad q[order], struct iir_biquad_state st[order]); -float run_biquad(float x, const struct iir_biquad *const q, struct iir_biquad_state *const restrict st); 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 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, uint64_t ts); #ifdef SIMULATION void debug_print_vector(const char *name, size_t len, const float *data, size_t stride, bool index, bool debug) { @@ -45,11 +51,19 @@ void debug_print_vector(const char *name, size_t len, const float *data, size_t void debug_print_vector(const char *name, size_t len, const float *data, size_t stride, bool index, bool debug) {} #endif +void dsss_demod_init(struct dsss_demod_state *st) { + memset(st, 0, sizeof(*st)); + matcher_init(st->matcher_cache); +} + #ifdef SIMULATION void dsss_demod_step(struct dsss_demod_state *st, float new_value, uint64_t ts, int record_channel) { + bool debug = false; + /* bool debug = (record_channel == -1) && (ts > 1000) && (ts % DSSS_CORRELATION_LENGTH == DSSS_CORRELATION_LENGTH-1); + */ if (debug) DEBUG_PRINT("Iteration %zd: signal=%f", ts, new_value); #else @@ -93,16 +107,19 @@ void dsss_demod_step(struct dsss_demod_state *st, float new_value) { for (size_t i=0; i<DSSS_GOLD_CODE_COUNT; i++) { float val = cwt[i] / avg[i]; - if (fabs(val) > DSSS_THESHOLD_FACTOR) - found = true; - if (fabs(val) > fabs(max_val)) { max_val = val; max_ch = i; max_ts = ts; + + if (fabs(val) > DSSS_THESHOLD_FACTOR) + found = true; } } + /* FIXME: skipped sample handling here */ + matcher_tick(st->matcher_cache, ts, max_ch, max_val); + if (found) { /* Continue ongoing group */ st->group.len++; @@ -120,6 +137,7 @@ void dsss_demod_step(struct dsss_demod_state *st, float new_value) { if (record_channel == -1) DEBUG_PRINT("GROUP FOUND: %8d len=%3d max=%f ch=%d offx=%d", ts, st->group.len, st->group.max, st->group.max_ch, st->group.max_ts); + group_received(st, ts); /* reset grouping state */ st->group.len = 0; @@ -128,49 +146,139 @@ void dsss_demod_step(struct dsss_demod_state *st, float new_value) { st->group.max = 0.0f; } -float score_group(const struct group *g, uint64_t ts) { - return fabs(g->max); /* Possibly at time penalty 1/(ts-max_ts) later */ +/* 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 + * + * 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) { + return (peak_ch<<1) | (peak_ampl > 0); } -ssize_t group_cache_insertion_index(const struct group *g, const struct group *cache, size_t cache_size, uint64_t ts) { - float min_score = INFINITY; - ssize_t min_idx = -1; - for (size_t i=0; i<cache_size; i++) { - /* If we find an empty or expired entry, use that */ - if (cache[i].max_ts == 0 || ts - cache[i].max_ts > group_cache_expiration) - return i; +void matcher_init(struct matcher_state states[static DSSS_MATCHER_CACHE_SIZE]) { + for (size_t i=0; i<DSSS_MATCHER_CACHE_SIZE; i++) + states[i].last_phase = -1; /* mark as inactive */ +} - /* Otherwise check weakest entry */ - float score = score_group(&cache[i]); - if (score < min_score) { - min_idx = i; - min_score = score; +/* TODO make these constants configurable from Makefile */ +const int group_phase_tolerance = (int)(DSSS_CORRELATION_LENGTH * 0.10); + +void matcher_tick(struct matcher_state states[static DSSS_MATCHER_CACHE_SIZE], uint64_t ts, int peak_ch, float peak_ampl) { + /* TODO make these constants configurable from Makefile */ + const float skip_sampling_depreciation = 0.2f; /* 0.0 -> no depreciation, 1.0 -> complete disregard */ + const float score_depreciation = 0.1f; /* 0.0 -> no depreciation, 1.0 -> complete disregard */ + const uint64_t current_phase = ts % DSSS_CORRELATION_LENGTH; + + for (size_t i=0; i<DSSS_MATCHER_CACHE_SIZE; i++) { + if (states[i].last_phase == -1) + continue; /* Inactive entry */ + + if (current_phase == states[i].last_phase) { + /* Skip sampling */ + float score = fabs(peak_ampl) * (1.0f - skip_sampling_depreciation); + if (score > states[i].candidate_score) { + /* We win, update candidate */ + states[i].candidate_score = score; + states[i].candidate_phase = current_phase; + states[i].candidate_data = decode_peak(peak_ch, peak_ampl); + } + } + + /* Note of caution on group_phase_tolerance: Group detection has some latency since a group is only considered + * "detected" after signal levels have fallen back below the detection threshold. This means we only get to + * process a group a couple ticks after its peak. We have to make sure the window is still open at this point. + * This means we have to match against group_phase_tolerance should a little bit loosely. + */ + if (abs(states[i].last_phase - current_phase) == group_phase_tolerance + DSSS_DECIMATION) { + /* Process window results */ + states[i].data[ states[i].data_pos ] = states[i].candidate_data; + states[i].data_pos = states[i].data_pos + 1; + states[i].last_score = score_depreciation * states[i].last_score + + (1.0f - score_depreciation) * states[i].candidate_score; + states[i].candidate_score = 0.0f; + + if (states[i].data_pos == TRANSMISSION_SYMBOLS) { + /* Frame received completely */ + DEBUG_PRINT("match on index %d phase %d score %.5f", i, states[i].last_phase, states[i].last_score); + handle_dsss_received(states[i].data); + states[i].last_phase = -1; /* invalidate entry */ + } } } +} + +static float gaussian(float a, float b, float c, float x) { + float n = x-b; + return a*expf(-n*n / (2.0f* c*c)); +} - /* Return weakest group if weaker than candidate */ - if (min_score < score_group(g)) - return min_idx; - return -1; +static float score_group(const struct group *g, int phase_delta) { + /* TODO make these constants configurable from Makefile */ + const float distance_func_phase_tolerance = 10.0f; + return fabsf(g->max) * gaussian(1.0f, 0.0f, distance_func_phase_tolerance, phase_delta); } void group_received(struct dsss_demod_state *st, uint64_t ts) { - /* TODO make these constants configurable from Makefile */ - const uint64_t group_cache_expiration = DSSS_CORRELATION_LENGTH * DSSS_GROUP_CACHE_SIZE; + static_assert(group_phase_tolerance > 10); /* FIXME debug, remove */ - /* Insert into group cache if space is available or there is a weaker entry to replace */ - ssize_t found = group_cache_insertion_index(&st->group, st->group_cache, DSSS_GROUP_CACHE_SIZE); - if (!found) - return; /* Nothing changed */ - st->group_cache[found] = st->group; + 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); - float mean_phase = 0.0; - for (size_t i=0; i<DSSS_GROUP_CACHE_SIZE; i++) - mean_phase += (st->group_cache[i].max_ts) % DSSS_CORRELATION_LENGTH; - mean_phase /= DSSS_GROUP_CACHE_SIZE; + float min_score = INFINITY; + ssize_t min_idx = -1; + ssize_t empty_idx = -1; + for (size_t i=0; i<DSSS_MATCHER_CACHE_SIZE; i++) { + /* Search for entries with matching phase */ + + /* This is the score of this group given the cached decoding at [i] */ + int phase_delta = st->matcher_cache[i].last_phase - group_phase; + if (abs(phase_delta) <= group_phase_tolerance) { + + float group_score = score_group(&st->group, phase_delta); + if (st->matcher_cache[i].candidate_score < group_score) { + st->matcher_cache[i].candidate_score = group_score; + st->matcher_cache[i].candidate_phase = group_phase; + st->matcher_cache[i].candidate_data = decode_peak(st->group.max_ch, st->group.max); + } + } - + /* Search for empty entries */ + if (st->matcher_cache[i].last_phase == -1) + empty_idx = i; + + /* Search for weakest entry */ + float score = st->matcher_cache[i].last_score; + if (score < min_score) { + min_idx = i; + min_score = score; + } + } + + /* If we found empty entries, replace one by a new decoding starting at this group */ + if (empty_idx >= 0) { + st->matcher_cache[empty_idx].last_phase = group_phase; + st->matcher_cache[empty_idx].candidate_score = base_score; + st->matcher_cache[empty_idx].last_score = base_score; + st->matcher_cache[empty_idx].candidate_phase = group_phase; + st->matcher_cache[empty_idx].candidate_data = decode_peak(st->group.max_ch, st->group.max); + st->matcher_cache[empty_idx].data_pos = 0; + } + + /* If the weakest decoding in cache is weaker than a new decoding starting here, replace it */ + if (min_score < base_score) { + assert(min_idx >= 0); + st->matcher_cache[min_idx].last_phase = group_phase; + st->matcher_cache[min_idx].candidate_score = base_score; + st->matcher_cache[min_idx].last_score = base_score; + st->matcher_cache[min_idx].candidate_phase = group_phase; + st->matcher_cache[min_idx].candidate_data = decode_peak(st->group.max_ch, st->group.max); + st->matcher_cache[min_idx].data_pos = 0; + } } float run_iir(const float x, const int order, const struct iir_biquad q[order], struct iir_biquad_state st[order]) { diff --git a/controller/fw/src/dsss_demod.h b/controller/fw/src/dsss_demod.h index 30ec3f7..b9e75d2 100644 --- a/controller/fw/src/dsss_demod.h +++ b/controller/fw/src/dsss_demod.h @@ -5,6 +5,9 @@ #define DSSS_GOLD_CODE_COUNT ((1<<DSSS_GOLD_CODE_NBITS) + 1) #define DSSS_CORRELATION_LENGTH (DSSS_GOLD_CODE_LENGTH * DSSS_DECIMATION) +/* 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) struct iir_biquad { @@ -20,22 +23,26 @@ struct cwt_iir_filter_state { struct iir_biquad_state st[3]; }; -struct { +struct group { int len; /* length of group in samples */ float max; /* signed value of largest peak in group on any channel */ uint64_t max_ts; /* absolute position of above peak */ int max_ch; /* channel (gold sequence index) of above peak */ -} group; +}; -struct decoder_state { - int last_phase; +struct matcher_state { + int last_phase; /* 0 .. DSSS_CORRELATION_LENGTH */ int candidate_phase; float last_score; float candidate_score; - uint8_t data[PAYLOAD_DATA_BYTE]; +#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]; int data_pos; + uint8_t candidate_data; }; struct dsss_demod_state { @@ -49,9 +56,13 @@ struct dsss_demod_state { struct group group; - struct group group_cache[DSSS_GROUP_CACHE_SIZE]; + struct matcher_state matcher_cache[DSSS_MATCHER_CACHE_SIZE]; }; + +extern void handle_dsss_received(uint8_t data[TRANSMISSION_SYMBOLS]); + +void dsss_demod_init(struct dsss_demod_state *st); #ifdef SIMULATION void dsss_demod_step(struct dsss_demod_state *st, float new_value, uint64_t ts, int record_channel); #else /* SIMULATION */ diff --git a/controller/fw/tools/dsss_demod_test.c b/controller/fw/tools/dsss_demod_test.c index 4370f80..d76a8bd 100644 --- a/controller/fw/tools/dsss_demod_test.c +++ b/controller/fw/tools/dsss_demod_test.c @@ -12,6 +12,16 @@ #include "dsss_demod.h" +void handle_dsss_received(uint8_t data[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)); + if (i+1 < TRANSMISSION_SYMBOLS) + printf(", "); + } + printf(" ]\n"); +} + void print_usage() { fprintf(stderr, "Usage: dsss_demod_test [test_data.bin] [optional recording channel number]\n"); } @@ -87,7 +97,7 @@ int main(int argc, char **argv) { fprintf(stderr, "Starting simulation.\n"); struct dsss_demod_state demod; - memset(&demod, 0, sizeof(demod)); + dsss_demod_init(&demod); for (size_t i=0; i<n_samples; i++) { //fprintf(stderr, "Iteration %zd/%zd\n", i, n_samples); dsss_demod_step(&demod, buf_f[i], i, record_channel); diff --git a/lab-windows/dsss_experiments-ber.ipynb b/lab-windows/dsss_experiments-ber.ipynb index 2d06233..c483cfc 100644 --- a/lab-windows/dsss_experiments-ber.ipynb +++ b/lab-windows/dsss_experiments-ber.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 121, "metadata": {}, "outputs": [], "source": [ @@ -14,11 +14,11 @@ "from collections import defaultdict\n", "import json\n", "\n", - "\n", "from matplotlib import pyplot as plt\n", "import matplotlib\n", "import numpy as np\n", "from scipy import signal as sig\n", + "from scipy import fftpack as fftpack\n", "import ipywidgets\n", "\n", "from tqdm.notebook import tqdm\n", @@ -29,7 +29,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 108, "metadata": {}, "outputs": [], "source": [ @@ -38,7 +38,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 109, "metadata": {}, "outputs": [], "source": [ @@ -47,7 +47,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 110, "metadata": {}, "outputs": [], "source": [ @@ -73,7 +73,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 111, "metadata": {}, "outputs": [], "source": [ @@ -93,7 +93,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 112, "metadata": {}, "outputs": [], "source": [ @@ -107,7 +107,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 113, "metadata": {}, "outputs": [ { @@ -127,7 +127,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 114, "metadata": {}, "outputs": [], "source": [ @@ -142,7 +142,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 115, "metadata": {}, "outputs": [], "source": [ @@ -162,7 +162,67 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 125, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "4330b0f8ceea4d5d922d2063a81554ca", + "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" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "ax.psd(colorednoise.powerlaw_psd_gaussian(1, 1000))\n", + "None" + ] + }, + { + "cell_type": "code", + "execution_count": 130, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "start 14880 end 24800 rec 29760\n" + ] + } + ], + "source": [ + "test_duration = 32\n", + "test_nbits = 5\n", + "test_signal_amplitude=2.0e-3\n", + "test_decimation=10\n", + "test_signal_amplitude = 200e-3\n", + "noise_level = 10e-3\n", + "\n", + "#test_data = np.random.RandomState(seed=0).randint(0, 2 * (2**test_nbits), test_duration)\n", + "#test_data = np.array([0, 1, 2, 3] * 50)\n", + "test_data = np.array(range(test_duration))\n", + "signal = np.repeat(modulate(test_data, test_nbits, pad=False) * 2.0 - 1, test_decimation) * test_signal_amplitude\n", + "noise = colorednoise.powerlaw_psd_gaussian(1, len(signal)*3) * noise_level\n", + "noise[int(1.5*len(signal)):][:len(signal)] += signal\n", + "print('start', int(1.5*len(signal)), 'end', int(1.5*len(signal))+len(signal), 'rec', len(noise))\n", + "\n", + "with open(f'dsss_test_signals/dsss_test_noiseless_padded.bin', 'wb') as f:\n", + " for e in noise:\n", + " f.write(struct.pack('<f', e))" + ] + }, + { + "cell_type": "code", + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -176,13 +236,13 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "373401133cfe408aa15738e48c58dfaa", + "model_id": "", "version_major": 2, "version_minor": 0 }, @@ -200,6 +260,46 @@ }, { "cell_type": "code", + "execution_count": 104, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "<ipython-input-104-abeb28a85dfa>:5: 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()\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "", + "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" + } + ], + "source": [ + "#nonlinear_distance_impl = lambda x: np.exp(-np.abs(x)/10) * x**4\n", + "nonlinear_distance_impl = lambda x: np.exp(-((x/10 - 0.5)%1 - 0.5)**2 / (2*1.2/10**2))\n", + "\n", + "def plot_distance_func_impl():\n", + " fig, ax = plt.subplots()\n", + " x = np.linspace(-30, 30, 10000)\n", + " ax.plot(x, nonlinear_distance_impl(x))\n", + "\n", + "plot_distance_func_impl()" + ] + }, + { + "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], |