path: root/controller/fw/src/dsss_demod.c

#include <unistd.h>
#include <stdbool.h>
#include <math.h>
#include <stdlib.h>
#include <assert.h>

#include <arm_math.h>

#include "freq_meas.h"
#include "sr_global.h"
#include "dsss_demod.h"
#include "simulation.h"

#include "generated/dsss_gold_code.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};

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) {
    if (!debug)
        return;

    if (index) {
        DEBUG_PRINTN("    %16s  [", "");
        for (size_t i=0; i<len; i++)
            DEBUG_PRINTN("%8d  ", i);
        DEBUG_PRINTN("]\n");
    }

    DEBUG_PRINTN("    %16s: [", name);
    for (size_t i=0; i<len; i++)
        DEBUG_PRINTN("%8.5f, ", data[i*stride]);
    DEBUG_PRINTN("]\n");
}
#else 
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
void dsss_demod_step(struct dsss_demod_state *st, float new_value) {
#endif

    //const float hole_patching_threshold = 0.01 * DSSS_CORRELATION_LENGTH;

    st->signal[st->signal_wpos] = new_value;
    st->signal_wpos = (st->signal_wpos + 1) % ARRAY_LENGTH(st->signal);

    /* use new, incremented wpos for gold_correlate_step as first element of old data in ring buffer */
    for (size_t i=0; i<DSSS_GOLD_CODE_COUNT; i++)
        st->correlation[i][st->correlation_wpos] = gold_correlate_step(i, st->signal, st->signal_wpos, false);

    debug_print_vector("correlation",
            DSSS_GOLD_CODE_COUNT, &st->correlation[0][st->correlation_wpos], DSSS_WAVELET_LUT_SIZE, true, debug);

    st->correlation_wpos = (st->correlation_wpos + 1) % ARRAY_LENGTH(st->correlation[0]);

    float cwt[DSSS_GOLD_CODE_COUNT];
    for (size_t i=0; i<DSSS_GOLD_CODE_COUNT; i++)
        cwt[i] = cwt_convolve_step(st->correlation[i], st->correlation_wpos);

    debug_print_vector("cwt", DSSS_GOLD_CODE_COUNT, cwt, 1, false, debug);

    float avg[DSSS_GOLD_CODE_COUNT];
    for (size_t i=0; i<DSSS_GOLD_CODE_COUNT; i++)
        avg[i] = run_iir(fabs(cwt[i]), ARRAY_LENGTH(cwt_filter_bq), cwt_filter_bq, st->cwt_filter[i].st);

    debug_print_vector("avg", DSSS_GOLD_CODE_COUNT, avg, 1, false, debug);

    if (record_channel != -1)
        DEBUG_PRINTN("%f, %f, %f\n",
                st->correlation[record_channel][st->correlation_wpos], cwt[record_channel], avg[record_channel]);

    float max_val = st->group.max;
    int max_ch = st->group.max_ch;
    int max_ts = st->group.max_ts;
    bool found = false;
    for (size_t i=0; i<DSSS_GOLD_CODE_COUNT; i++) {
        float val = cwt[i] / avg[i];

        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++;
        st->group.max = max_val;
        st->group.max_ch = max_ch;
        st->group.max_ts = max_ts;
        return;
    }
    
    if (st->group.len == 0)
        /* We're between groups */
        return;

    /* A group ended. Process result. */
    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;
    st->group.max_ts = 0;
    st->group.max_ch = 0;
    st->group.max = 0.0f;
}

/* 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);
}

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 */
}

/* 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));
}


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) {
    static_assert(group_phase_tolerance > 10); /* FIXME debug, remove */

    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 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]) {
    float intermediate = x;
    for (int i=0; i<(order+1)/2; i++)
        intermediate = run_biquad(intermediate, &q[i], &st[i]);
    return intermediate;
}

float run_biquad(float x, const struct iir_biquad *const q, struct iir_biquad_state *const restrict st) {
    /* direct form 2, see https://en.wikipedia.org/wiki/Digital_biquad_filter */
    float intermediate = x + st->reg[0] * -q->a[0] + st->reg[1] * -q->a[1];
    float out = intermediate * q->b[0] + st->reg[0] * q->b[1] + st->reg[1] * q->b[2];
    st->reg[1] = st->reg[0];
    st->reg[0] = intermediate;
    return out;
}

float cwt_convolve_step(const float v[DSSS_WAVELET_LUT_SIZE], size_t offx) {
    float sum = 0.0f;
    for (ssize_t j=0; j<DSSS_WAVELET_LUT_SIZE; j++) {
        /* Our wavelet is symmetric so convolution and correlation are identical. Use correlation here for ease of
         * implementation */
        sum += v[(offx + j) % DSSS_WAVELET_LUT_SIZE] * dsss_cwt_wavelet_table[j];
        //DEBUG_PRINT("        j=%d v=%f w=%f", j, v[(offx + j) % DSSS_WAVELET_LUT_SIZE], dsss_cwt_wavelet_table[j]);
    }
    return sum;
}

/* Compute last element of correlation for input [a] and hard-coded gold sequences.
 *
 * This is intened to be used once for each new incoming sample in [a]. It expects [a] to be of length
 * [dsss_correlation_length] and produces the one sample where both the reference sequence and the input fully overlap.
 * This is equivalent to "valid" mode in numpy's terminology[0].
 *
 * [0] https://docs.scipy.org/doc/numpy/reference/generated/numpy.correlate.html
 */
float gold_correlate_step(const size_t ncode, const float a[DSSS_CORRELATION_LENGTH], size_t offx, bool debug) {

    float acc_outer = 0.0f;
    uint8_t table_byte = 0;
    if (debug) DEBUG_PRINTN("Correlate n=%d: ", ncode);
    for (size_t i=0; i<DSSS_GOLD_CODE_LENGTH; i++) {

        if ((i&7) == 0) {
            table_byte = dsss_gold_code_table[ncode][i>>3]; /* Fetch sequence table item */
            if (debug) DEBUG_PRINTN("|");
        }
        int bv = table_byte & (0x80>>(i&7)); /* Extract bit */
        bv = !!bv*2 - 1; /* Map 0, 1 -> -1, 1 */
        if (debug) DEBUG_PRINTN("%s%d\033[0m", bv == 1 ? "\033[92m" : "\033[91m", (bv+1)/2);

        float acc_inner = 0.0f;
        for (size_t j=0; j<DSSS_DECIMATION; j++)
            acc_inner += a[(offx + i*DSSS_DECIMATION + j) % DSSS_CORRELATION_LENGTH]; /* Multiply item */
        //if (debug) DEBUG_PRINTN("%.2f ", acc_inner);
        acc_outer += acc_inner * bv;
    }
    if (debug) DEBUG_PRINTN("\n");
    return acc_outer / DSSS_CORRELATION_LENGTH;
}