Working on DSSS demodulator sim

This commit is contained in:
jaseg 2020-03-05 19:15:28 +01:00
parent d9b26d16c0
commit 4b419bd1ad
14 changed files with 94312 additions and 62 deletions

View file

@ -6,7 +6,7 @@ LIBSODIUM_DIR ?= libsodium
TINYAES_DIR ?= tinyaes
MUSL_DIR ?= musl
C_SOURCES := src/main.c src/mspdebug_wrapper.c src/spi_flash.c src/freq_meas.c
C_SOURCES := src/main.c src/mspdebug_wrapper.c src/spi_flash.c src/freq_meas.c src/dsss_demod.c
C_SOURCES += $(MSPDEBUG_DIR)/drivers/jtaglib.c
CMSIS_SOURCES += $(CMSIS_DIR)/CMSIS/DSP/Source/TransformFunctions/arm_rfft_fast_init_f32.c
@ -41,13 +41,17 @@ OPENCM3_LIB := opencm3_stm32f4
PREFIX ?= arm-none-eabi-
GOLD_CODE_NBITS ?= 5
FMEAS_ADC_SAMPLING_RATE ?= 1000
FMEAS_ADC_MAX ?= 4096
FMEAS_FFT_LEN ?= 256
FMEAS_FFT_WINDOW ?= gaussian
FMEAS_FFT_WINDOW_SIGMA ?= 16.0
DSSS_GOLD_CODE_NBITS ?= 5
DSSS_DECIMATION ?= 10
DSSS_THESHOLD_FACTOR ?= 4.0f
DSSS_WAVELET_WIDTH ?= 7.3
DSSS_WAVELET_LUT_SIZE ?= 69
CC := $(PREFIX)gcc
CXX := $(PREFIX)g++
@ -78,6 +82,7 @@ MUSL_DIR_ABS := $(abspath $(MUSL_DIR))
COMMON_CFLAGS += -I$(OPENCM3_DIR_ABS)/include -Imspdebug/util -Imspdebug/drivers -Ilevmarq
COMMON_CFLAGS += -I$(CMSIS_DIR_ABS)/CMSIS/DSP/Include -I$(CMSIS_DIR_ABS)/CMSIS/Core/Include
CFLAGS += -I$(abspath musl_include_shims)
COMMON_CFLAGS += -I$(BUILDDIR) -Isrc
COMMON_CFLAGS += -Os -std=gnu11 -g -DSTM32F4
CFLAGS += -mthumb -mcpu=cortex-m4 -mfloat-abi=hard -mfpu=fpv4-sp-d16
@ -85,13 +90,18 @@ CFLAGS += -mthumb -mcpu=cortex-m4 -mfloat-abi=hard -mfpu=fpv4-sp-d16
# Note: libopencm3 requires some standard libc definitions from stdint.h and stdbool.h, so we don't pass -nostdinc here.
CFLAGS += -nostdlib -ffreestanding
CFLAGS += -fno-common -ffunction-sections -fdata-sections
COMMON_CFLAGS += -DGOLD_CODE_NBITS=$(GOLD_CODE_NBITS) -DFMEAS_FFT_LEN=$(FMEAS_FFT_LEN) -DFMEAS_ADC_MAX=$(FMEAS_ADC_MAX)
COMMON_CFLAGS += -DDSSS_GOLD_CODE_NBITS=$(DSSS_GOLD_CODE_NBITS) -DFMEAS_FFT_LEN=$(FMEAS_FFT_LEN) -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_WAVELET_WIDTH=$(DSSS_WAVELET_WIDTH)
COMMON_CFLAGS += -DDSSS_WAVELET_LUT_SIZE=$(DSSS_WAVELET_LUT_SIZE)
# for musl
CFLAGS += -Dhidden=
SIM_CFLAGS += -Isrc -lm -DSIMULATION
SIM_CFLAGS += -lm -DSIMULATION
INT_CFLAGS += -Wall -Wextra -Wpedantic -Wshadow -Wimplicit-function-declaration -Wundef
INT_CFLAGS += -Wredundant-decls -Wmissing-prototypes -Wstrict-prototypes
@ -112,7 +122,24 @@ LDFLAGS += -L$(OPENCM3_DIR_ABS)/lib -l$(OPENCM3_LIB) $(shell $(CC) -print-libg
all: $(BUILDDIR)/$(BINARY)
tests: $(BUILDDIR)/tools/freq_meas_test
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 $< $@
.PRECIOUS: $(BUILDDIR)/generated/gold_code_%.c
$(BUILDDIR)/generated/gold_code_%.c $(BUILDDIR)/generated/gold_code_%.h: | $(BUILDDIR)/generated
$(PYTHON3) tools/gold_code_header_gen.py -v dsss_gold_code_table -$(patsubst .%,%,$(suffix $@)) $* > $@
.PRECIOUS: $(BUILDDIR)/generated/fmeas_fft_window.c
$(BUILDDIR)/generated/fmeas_fft_window.c: | $(BUILDDIR)/generated
$(PYTHON3) tools/fft_window_header_gen.py -v fmeas_fft_window_table $(FMEAS_FFT_WINDOW) $(FMEAS_FFT_LEN) $(FMEAS_FFT_WINDOW_SIGMA) > $@
.PRECIOUS: $(BUILDDIR)/generated/dsss_cwt_wavelet.c
$(BUILDDIR)/generated/dsss_cwt_wavelet.c: | $(BUILDDIR)/generated
$(PYTHON3) tools/cwt_wavelet_header_gen.py -v dsss_cwt_wavelet_table $(DSSS_WAVELET_LUT_SIZE) $(DSSS_WAVELET_WIDTH) > $@
$(BUILDDIR)/generated: ; mkdir -p $@
OBJS := $(addprefix $(BUILDDIR)/,$(C_SOURCES:.c=.o) $(CXX_SOURCES:.cpp=.o))
@ -121,16 +148,23 @@ ALL_OBJS += $(OPENCM3_DIR)/lib/lib$(OPENCM3_LIB).a
ALL_OBJS += $(BUILDDIR)/libsodium/src/libsodium/.libs/libsodium.a
ALL_OBJS += $(BUILDDIR)/tinyaes/aes.o
ALL_OBJS += $(BUILDDIR)/levmarq/levmarq.o
ALL_OBJS += $(BUILDDIR)/generated/gold_code_$(GOLD_CODE_NBITS).o
ALL_OBJS += $(BUILDDIR)/generated/gold_code_$(DSSS_GOLD_CODE_NBITS).o
ALL_OBJS += $(BUILDDIR)/generated/fmeas_fft_window.o
ALL_OBJS += $(BUILDDIR)/generated/dsss_cwt_wavelet.o
$(BUILDDIR)/$(BINARY): $(ALL_OBJS)
$(LD) -T$(LDSCRIPT) $(COMMON_LDFLAGS) $(LDFLAGS) -o $@ -Wl,-Map=$(BUILDDIR)/src/$*.map $^
tools: $(BUILDDIR)/tools/freq_meas_test
$(BUILDDIR)/tools/freq_meas_test: tools/freq_meas_test.c src/freq_meas.c levmarq/levmarq.c $(BUILDDIR)/generated/fmeas_fft_window.c $(CMSIS_SOURCES)
mkdir -p $(@D)
$(HOST_CC) $(COMMON_CFLAGS) $(SIM_CFLAGS) -o $@ $^
tools: $(BUILDDIR)/tools/dsss_demod_test
$(BUILDDIR)/tools/dsss_demod_test: tools/dsss_demod_test.c src/dsss_demod.c $(BUILDDIR)/generated/dsss_cwt_wavelet.c $(BUILDDIR)/generated/gold_code_$(DSSS_GOLD_CODE_NBITS).c
mkdir -p $(@D)
$(HOST_CC) $(COMMON_CFLAGS) $(SIM_CFLAGS) -o $@ $^
$(BUILDDIR)/src/%.o: src/%.c
mkdir -p $(@D)
$(CC) $(COMMON_CFLAGS) $(CFLAGS) $(INT_CFLAGS) -o $@ -c $<
@ -158,14 +192,6 @@ $(BUILDDIR)/tinyaes/aes.o:
mkdir -p $(@D)
make -C $(@D) -f $(TINYAES_DIR_ABS)/Makefile VPATH=$(TINYAES_DIR_ABS) CFLAGS="$(COMMON_CFLAGS) $(CFLAGS) -c" CC=$(CC) LD=$(LD) AR=$(AR) aes.o
$(BUILDDIR)/generated/gold_code_%.c:
mkdir -p $(@D)
$(PYTHON3) tools/gold_code_header_gen.py $* > $@
$(BUILDDIR)/generated/fmeas_fft_window.c:
mkdir -p $(@D)
$(PYTHON3) tools/fft_window_header_gen.py -v fmeas_fft_window_table $(FMEAS_FFT_WINDOW) $(FMEAS_FFT_LEN) $(FMEAS_FFT_WINDOW_SIGMA) > $@
$(BUILDDIR)/musl/lib/libc.a:
mkdir -p $(BUILDDIR)/musl
cd $(BUILDDIR)/musl && CFLAGS="$(SIM_CFLAGS) $(COMMON_CFLAGS)" CC=$(CC) LD=$(LD) AR=$(AR) $(MUSL_DIR_ABS)/configure && $(MAKE) TARGET=arm-linux-musleabihf GCC_CONFIG="-mcpu=cortex-m4 -mfloat-abi=soft" -j $(shell nproc)
@ -174,13 +200,13 @@ build/ldpc_decoder_test.so: src/ldpc_decoder.c
gcc -fPIC -shared -Wall -Wextra -Wpedantic -std=gnu11 -O0 -g -o $@ $^
clean:
-rm -r $(BUILDDIR)/src
-rm -r $(BUILDDIR)/generated
-rm $(BUILDDIR)/$(BINARY)
-rm $(BUILDDIR)/tools/freq_meas_test
rm -rf $(BUILDDIR)/src
rm -rf $(BUILDDIR)/generated
rm -f $(BUILDDIR)/$(BINARY)
rm -f $(BUILDDIR)/tools/freq_meas_test
mrproper: clean
-rm -r build
rm -rf build
-$(MAKE) -C $(OPENCM3_DIR) clean
.PHONY: clean mrproper

View file

@ -0,0 +1,180 @@
#include <unistd.h>
#include <stdbool.h>
#include <math.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"
/* Generated CWT wavelet LUT */
extern const float * const dsss_cwt_wavelet_table;
struct iir_biquad cwt_filter_bq[3] = {
{.a = {-1.993939440f, 0.993949280f}, .b = {0.2459934300683e-5, 0.4919868601367e-5, 0.2459934300683e-5}},
{.a = {-1.995557124f, 0.995566972f}, .b = {0.2461930046414e-5, 0.4923860092828e-5, 0.2461930046414e-5}},
{.a = {-1.998365254f, 0.998375115f}, .b = {0.2465394452097e-5, 0.4930788904195e-5, 0.2465394452097e-5}},
};
float gold_correlate_step(const size_t ncode, const float a[DSSS_CORRELATION_LENGTH], size_t offx);
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);
#ifdef SIMULATION
void dsss_demod_step(struct dsss_demod_state *st, float new_value, size_t sim_pos) {
#else /* SIMULATION */
void dsss_demod_step(struct dsss_demod_state *st, float new_value) {
#endif /* SIMULATION */
//#define DEBUG_PRINT(...) ((void)0)
//#define DEBUG_PRINTN(...) ((void)0)
bool debug = sim_pos % DSSS_CORRELATION_LENGTH == 0;
if (debug) DEBUG_PRINT("Iteration %zd", sim_pos);
//const float peak_group_threshold = 0.05 * DSSS_CORRELATION_LENGTH;
//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);
if (debug) DEBUG_PRINT(" signal: %f", new_value);
/* 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->correlation_wpos);
/* debug */
/*
DEBUG_PRINTN(" correlation: [");
for (size_t i=0; i<DSSS_GOLD_CODE_COUNT; i++)
DEBUG_PRINTN("%f, ", st->correlation[i][st->correlation_wpos]);
DEBUG_PRINTN("]\n");
*/
/* end */
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 */
if (debug) DEBUG_PRINTN(" cwt: [");
for (size_t i=0; i<DSSS_GOLD_CODE_COUNT; i++)
if (debug) DEBUG_PRINTN("%f, ", cwt[i]);
if (debug) DEBUG_PRINTN("]\n");
/* end */
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 */
/*
DEBUG_PRINTN(" avg: [");
for (size_t i=0; i<DSSS_GOLD_CODE_COUNT; i++)
DEBUG_PRINTN("%f, ", avg[i]);
DEBUG_PRINTN("]\n");
*/
/* end */
float max_val = st->group.max;
int max_ch = st->group.max_ch;
int max_idx = st->group.max_idx + 1;
bool found = false;
if (debug) DEBUG_PRINTN(" rel: [");
for (size_t i=0; i<DSSS_GOLD_CODE_COUNT; i++) {
float val = cwt[i] / avg[i];
if (debug) DEBUG_PRINTN("%f, ", val);
if (fabs(val) > DSSS_THESHOLD_FACTOR)
found = true;
if (fabs(val) > fabs(max_val)) {
max_val = val;
max_ch = i;
max_idx = st->group.len;
}
}
if (debug) DEBUG_PRINTN("]\n");
/* debug */
if (debug) DEBUG_PRINT(" found=%d len=%d idx=%d ch=%d max=%f",
found, st->group.len, st->group.max_idx, st->group.max_ch, st->group.max);
/* end */
if (found) {
/* Continue ongoing group */
st->group.len++;
st->group.max = max_val;
st->group.max_ch = max_ch;
st->group.max_idx = max_idx;
return;
}
if (st->group.len == 0)
/* We're between groups */
return;
/* A group ended. Process result. */
DEBUG_PRINT("GROUP FOUND: %8d len=%3d max=%f ch=%d offx=%d",
sim_pos, st->group.len, st->group.max, st->group.max_ch, st->group.max_idx);
/* reset grouping state */
st->group.len = 0;
st->group.max_idx = 0;
st->group.max_ch = 0;
st->group.max = 0.0f;
}
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 / DSSS_WAVELET_LUT_SIZE;
}
/* 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) {
float acc_outer = 0.0f;
uint8_t table_byte = 0;
for (size_t i=0, pos=0; i<DSSS_GOLD_CODE_LENGTH; i++, pos += DSSS_DECIMATION) {
float acc_inner = 0.0f;
for (size_t j=0; j<DSSS_DECIMATION; j++) {
if ((pos&7) == 0)
table_byte = dsss_gold_code_table[ncode][pos>>3]; /* Fetch sequence table item */
int bv = table_byte & (1<<(pos&7)); /* Extract bit */
bv = !!bv*2 - 1; /* Map 0, 1 -> -1, 1 */
acc_inner += a[(offx + i + j) % DSSS_CORRELATION_LENGTH] * bv; /* Multiply item */
}
acc_outer += acc_inner;
}
return acc_outer / DSSS_CORRELATION_LENGTH;
}

View file

@ -0,0 +1,44 @@
#ifndef __DSSS_DEMOD_H__
#define __DSSS_DEMOD_H__
#define DSSS_GOLD_CODE_LENGTH ((1<<DSSS_GOLD_CODE_NBITS) - 1)
#define DSSS_GOLD_CODE_COUNT ((1<<DSSS_GOLD_CODE_NBITS) + 1)
#define DSSS_CORRELATION_LENGTH (DSSS_GOLD_CODE_LENGTH * DSSS_DECIMATION)
struct iir_biquad {
float a[2];
float b[3];
};
struct iir_biquad_state {
float reg[2];
};
struct cwt_iir_filter_state {
struct iir_biquad_state st[3];
};
struct dsss_demod_state {
float signal[DSSS_CORRELATION_LENGTH];
size_t signal_wpos;
float correlation[DSSS_GOLD_CODE_COUNT][DSSS_WAVELET_LUT_SIZE];
size_t correlation_wpos;
struct cwt_iir_filter_state cwt_filter[DSSS_GOLD_CODE_COUNT];
struct {
int len; /* length of group in samples */
float max; /* signed value of largest peak in group on any channel */
int max_idx; /* position of above peak counted from start of group */
int max_ch; /* channel (gold sequence index) of above peak */
} group;
};
#ifdef SIMULATION
void dsss_demod_step(struct dsss_demod_state *st, float new_value, size_t sim_pos);
#else /* SIMULATION */
void dsss_demod_step(struct dsss_demod_state *st, float new_value);
#endif /* SIMULATION */
#endif /* __DSSS_DEMOD_H__ */

View file

@ -2,6 +2,7 @@
#define __SR_GLOBAL_H__
#define UNUSED(x) ((void) x)
#define ARRAY_LENGTH(x) (sizeof(x) / sizeof(x[0]))
static inline uint16_t htole(uint16_t val) { return val; }

View file

@ -0,0 +1,29 @@
#!/usr/bin/env python3
import textwrap
import scipy.signal as sig
import numpy as np
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('n', type=int, help='Window size')
parser.add_argument('w', type=float, help='Wavelet width')
parser.add_argument('-v', '--variable', default='cwt_ricker_table', help='Name for alias variable pointing to generated wavelet LUT')
args = parser.parse_args()
print(f'/* CWT Ricker wavelet LUT for {args.n} sample window of width {args.w}. */')
varname = f'cwt_ricker_{args.n}_window_{str(args.w).replace(".", "F")}'
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)
print(textwrap.fill(par,
initial_indent=' '*4, subsequent_indent=' '*4,
width=120,
replace_whitespace=False, drop_whitespace=False))
print('};')
print()
print(f'const float * const {args.variable} __attribute__((weak)) = {varname};')

View file

@ -0,0 +1,84 @@
#include <stdint.h>
#include <math.h>
#include <unistd.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <sys/fcntl.h>
#include "dsss_demod.h"
void print_usage() {
fprintf(stderr, "Usage: dsss_demod_test [test_data.bin]\n");
}
int main(int argc, char **argv) {
if (argc != 2) {
fprintf(stderr, "Error: Invalid arguments.\n");
print_usage();
return 1;
}
int fd = open(argv[1], O_RDONLY);
struct stat st;
if (fstat(fd, &st)) {
fprintf(stderr, "Error querying test data file size: %s\n", strerror(errno));
return 2;
}
if (st.st_size < 0 || st.st_size > 1000000) {
fprintf(stderr, "Error reading test data: too much test data (size=%zd)\n", st.st_size);
return 2;
}
if (st.st_size % sizeof(float) != 0) {
fprintf(stderr, "Error reading test data: file size is not divisible by %zd (size=%zd)\n", sizeof(float), st.st_size);
return 2;
}
char *buf = malloc(st.st_size);
if (!buf) {
fprintf(stderr, "Error allocating memory");
return 2;
}
fprintf(stderr, "Reading %zd samples test data...", st.st_size/sizeof(float));
size_t nread = 0;
while (nread < st.st_size) {
ssize_t rc = read(fd, buf, st.st_size - nread);
if (rc == -EINTR || rc == -EAGAIN)
continue;
if (rc < 0) {
fprintf(stderr, "\nError reading test data: %s\n", strerror(errno));
return 2;
}
if (rc == 0) {
fprintf(stderr, "\nError reading test data: Unexpected end of file\n");
return 2;
}
nread += rc;
}
fprintf(stderr, " done.\n");
const size_t n_samples = st.st_size / sizeof(float);
float *buf_f = (float *)buf;
fprintf(stderr, "Starting simulation.\n");
struct dsss_demod_state demod;
memset(&demod, 0, sizeof(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);
}
return 0;
}

View file

@ -0,0 +1,39 @@
#!/usr/bin/env python3
import os
from os import path
import subprocess
import json
import numpy as np
np.set_printoptions(linewidth=240)
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')
args = parser.parse_args()
bin_files = [ path.join(args.dir, d) for d in os.listdir(args.dir) if d.lower().endswith('.bin') ]
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)
# 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)))
print(f'{path.basename(p):<60}: mean={mean:<8.4f}Hz rms={rms*1000:.3f}mHz')
if args.dump:
with open(args.dump, 'w') as f:
json.dump(savedata, f)

View file

@ -70,7 +70,7 @@ int main(int argc, char **argv) {
}
fprintf(stderr, " done.\n");
size_t n_samples = st.st_size / sizeof(float);
const size_t n_samples = st.st_size / sizeof(float);
float *buf_f = (float *)buf;
int16_t *sim_adc_buf = calloc(sizeof(int16_t), n_samples);

View file

@ -1,6 +1,9 @@
#!/usr/bin/env python3
import sys
import math
import textwrap
import contextlib
import numpy as np
import scipy.signal as sig
@ -24,22 +27,50 @@ def gold(n):
(seq0, _st0), (seq1, _st1) = sig.max_len_seq(n, taps=t0), sig.max_len_seq(n, taps=t1)
return gen_gold(seq0, seq1)
@contextlib.contextmanager
def print_include_guards(macro_name):
print(f'#ifndef {macro_name}')
print(f'#define {macro_name}')
yield
print(f'#endif /* {macro_name} */')
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser()
parser = argparse.ArgumentParser(add_help=False)
parser.add_argument('n', type=int, choices=preferred_pairs, help='bit width of shift register. Generate 2**n + 1 sequences of length 2**n - 1.')
parser.add_argument('-v', '--variable', default='gold_code_table', help='Name for weak alias of generated table')
parser.add_argument('-h', '--header', action='store_true', help='Generate header file')
parser.add_argument('-c', '--source', action='store_true', help='Generate table source file')
args = parser.parse_args()
print('#include <unistd.h>')
print()
print(f'/* {args.n} bit gold sequences: {2**args.n+1} sequences of length {2**args.n-1} bit.')
print(f' *')
print(f' * Each code is packed left-aligned into {2**args.n // 8} bytes in big-endian byte order.')
print(f' */')
print(f'const uint8_t gold_code_{args.n}bit[] = {{')
for i, code in enumerate(gold(args.n)):
par = ' '.join(f'0x{d:02x},' for d in np.packbits(code)) + f' /* {i: 3d} "{"".join(str(x) for x in code)}" */'
print(textwrap.fill(par, initial_indent=' '*4, subsequent_indent=' '*4, width=120))
print('};')
print()
print(f'const uint8_t * const gold_code_table __attribute__((weak)) = gold_code_{args.n}bit;')
if not args.header != args.source:
print('Exactly one of --header and --source must be given.', file=sys.stderr)
sys.exit(1)
nbytes = math.ceil((2**args.n-1)/8)
if args.source:
print('/* THIS IS A GENERATED FILE. DO NOT EDIT! */')
print('#include <unistd.h>')
print('#include <stdint.h>')
print()
print(f'/* {args.n} bit gold sequences: {2**args.n+1} sequences of length {2**args.n-1} bit.')
print(f' *')
print(f' * Each code is packed left-aligned into {nbytes} bytes in big-endian byte order.')
print(f' */')
print(f'const uint8_t gold_code_{args.n}bit[{2**args.n+1}][{nbytes}] = {{')
for i, code in enumerate(gold(args.n)):
par = '{' + ' '.join(f'0x{d:02x},' for d in np.packbits(code)) + f'}}, /* {i: 3d} "{"".join(str(x) for x in code)}" */'
print(textwrap.fill(par, initial_indent=' '*4, subsequent_indent=' '*4, width=120))
print('};')
print()
print(f'const uint8_t * const {args.variable} __attribute__((weak)) = (uint8_t *const)gold_code_{args.n}bit;')
print(f'const size_t {args.variable}_nbits __attribute__((weak)) = {args.n};')
else:
print('/* THIS IS A GENERATED FILE. DO NOT EDIT! */')
with print_include_guards(f'__GOLD_CODE_GENERATED_HEADER_{args.n}__'):
print(f'extern const uint8_t gold_code_{args.n}bit[{2**args.n+1}][{nbytes}];')
with print_include_guards(f'__GOLD_CODE_GENERATED_HEADER_GLOBAL_SYM_{args.variable.upper()}__'):
print(f'extern const uint8_t {args.variable}[{2**args.n+1}][{nbytes}];')
print(f'extern const size_t {args.variable}_nbits;')

View file

@ -137,6 +137,25 @@
" return test_data, signal + noise"
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [],
"source": [
"test_duration = 300\n",
"test_nbits = 5\n",
"test_signal_amplitude=2.0e-3\n",
"test_decimation=10\n",
"\n",
"for test_signal_amplitude in [2.0e-3, 20e-3, 200e-3, 2]:\n",
" test_data = np.random.RandomState(seed=0).randint(0, 2 * (2**test_nbits), test_duration)\n",
" signal = np.repeat(modulate(test_data, test_nbits) * 2.0 - 1, test_decimation) * test_signal_amplitude\n",
" with open(f'dsss_test_signals/dsss_test_noiseless_{test_signal_amplitude*1000:.0f}mHz.bin', 'wb') as f:\n",
" for e in signal:\n",
" f.write(struct.pack('<f', e))"
]
},
{
"cell_type": "code",
"execution_count": 9,
@ -159,7 +178,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "c5fba7c7cec24f09a4f23f7c7d87eb90",
"model_id": "373401133cfe408aa15738e48c58dfaa",
"version_major": 2,
"version_minor": 0
},
@ -184,6 +203,129 @@
"noprint = lambda *args, **kwargs: None"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-44-38927d86fcc7>:1: 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": "2015d1f2951340d49397c6c289096b01",
"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"
},
{
"data": {
"text/plain": [
"Text(0.5, 1.0, 'Ricker wavelet, w=69 a=7.3')"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"w = 69\n",
"a = 7.3\n",
"ax.plot(range(-w//2+1, w//2+1), sig.ricker(w, a))\n",
"ax.grid()\n",
"ax.axvline(0, color='orange')\n",
"ax.set_title(f'Ricker wavelet, w={w} a={a}')"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-56-8e404350b85f>:1: 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": "3cd0c25b914b41df835005ef6fb51b34",
"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",
"r = list(range(60, 120))\n",
"ax.plot(r, [sum(sig.ricker(w, a)) for w in r])\n",
"ax.set_yscale('log')\n",
"ax.grid()"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-67-2097444ed7d2>:1: 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": "c5d1587fa0c944378aed31941ab66ff6",
"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",
"sw = 256\n",
"w = sig.ricker(sw, a)\n",
"r = list(range(1, sw//2 - 10))\n",
"d = [-sum(w[:i]) - sum(w[-i:]) for i in r]\n",
"ax.plot([sw-2*x for x in r], d)\n",
"ax.set_yscale('log')\n",
"ax.grid()"
]
},
{
"cell_type": "code",
"execution_count": 12,
@ -279,7 +421,7 @@
" chain_candidates = []\n",
" for chain_score, chain in candidates:\n",
" pos, ampl, _idx = chain[-1]\n",
" score_fun = lambda pos, npos, npol: pol_score_factor*abs(npol)/avg_peak + nonlinear_distance((npos-pos)/bit_period)\n",
" score_fun = lambda pos, npos, npol: pol_score_eebfactor*abs(npol)/avg_peak + nonlinear_distance((npos-pos)/bit_period)\n",
" 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)\n",
"\n",
" print(f' candidates for {pos}, {ampl}:')\n",
@ -366,13 +508,13 @@
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "5106e62521e246a0b0766cce3f9e556f",
"model_id": "672fae23dc99473a967c46b4b90d09d7",
"version_major": 2,
"version_minor": 0
},
@ -386,10 +528,10 @@
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fc200d359d0>]"
"[<matplotlib.lines.Line2D at 0x7fc38cb743d0>]"
]
},
"execution_count": 14,
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
@ -444,7 +586,7 @@
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
@ -454,13 +596,13 @@
},
{
"cell_type": "code",
"execution_count": 16,
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "15e67925b1bf47b5bbde289b353ed201",
"model_id": "985b1935796343d19c27e0b0ae31363c",
"version_major": 2,
"version_minor": 0
},
@ -474,7 +616,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "6432dbccc7db43d7a227384bca772716",
"model_id": "418a9dd281534ffa8e0aa5f250424127",
"version_major": 2,
"version_minor": 0
},
@ -491,12 +633,12 @@
"text": [
"nbits=5\n",
"nbits=6\n",
"signal_amplitude=0.00046: ser=0.97396 ±0.00737, br=14.64844\n",
"signal_amplitude=0.00015: ser=0.98958 ±0.00737, br=5.85938\n",
"signal_amplitude=0.00022: ser=0.96875 ±0.01276, br=17.57812\n",
"signal_amplitude=0.00010: ser=0.97396 ±0.00737, br=14.64844\n",
"signal_amplitude=0.00032: ser=0.97396 ±0.00737, br=14.64844\n",
"signal_amplitude=0.00046: ser=0.97396 ±0.00737, br=14.64844\n",
"signal_amplitude=0.00010: ser=0.97396 ±0.00737, br=14.64844\n",
"signal_amplitude=0.00068: ser=0.94792 ±0.01949, br=29.29688\n",
"signal_amplitude=0.00022: ser=0.96875 ±0.01276, br=17.57812\n",
"signal_amplitude=0.00100: ser=0.78125 ±0.02552, br=123.04688\n",
"signal_amplitude=0.00147: ser=0.28646 ±0.03683, br=401.36719\n",
"signal_amplitude=0.00215: ser=0.06250 ±0.03375, br=527.34375\n",
@ -505,8 +647,8 @@
"signal_amplitude=0.00010: ser=0.98958 ±0.00737, br=3.51562\n",
"signal_amplitude=0.00022: ser=0.98438 ±0.00000, br=5.27344\n",
"signal_amplitude=0.00032: ser=0.98438 ±0.00000, br=5.27344\n",
"signal_amplitude=0.00068: ser=0.68229 ±0.09051, br=107.22656\n",
"signal_amplitude=0.00046: ser=0.97917 ±0.01949, br=7.03125\n",
"signal_amplitude=0.00068: ser=0.68229 ±0.09051, br=107.22656\n",
"signal_amplitude=0.00100: ser=0.15104 ±0.05156, br=286.52344\n",
"signal_amplitude=0.00147: ser=0.01562 ±0.00000, br=332.22656\n",
"signal_amplitude=0.00215: ser=0.01562 ±0.00000, br=332.22656\n",
@ -519,10 +661,10 @@
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7fc23a223550>"
"<matplotlib.legend.Legend at 0x7fc382e9b7f0>"
]
},
"execution_count": 16,
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
@ -595,13 +737,13 @@
},
{
"cell_type": "code",
"execution_count": 17,
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "bda9f8070b0142928936a7752f318e97",
"model_id": "8b9593f48de94ec399fd8e61dc11d1de",
"version_major": 2,
"version_minor": 0
},
@ -615,10 +757,10 @@
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7fc20046b730>"
"<matplotlib.legend.Legend at 0x7fc3c4c54b20>"
]
},
"execution_count": 17,
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
@ -651,13 +793,13 @@
},
{
"cell_type": "code",
"execution_count": 18,
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "33b1f777dda149eda3a2f031087372c9",
"model_id": "15dc2e1a7cdd460e8cbde0352ea18270",
"version_major": 2,
"version_minor": 0
},
@ -861,7 +1003,7 @@
},
{
"cell_type": "code",
"execution_count": 19,
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
@ -870,13 +1012,13 @@
},
{
"cell_type": "code",
"execution_count": 25,
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "406505452bb84597a6c39bc12267dcc6",
"model_id": "e381ce02f19d4b85bbbaa9ea99e65623",
"version_major": 2,
"version_minor": 0
},
@ -891,7 +1033,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-25-50856f807a2a>:20: RuntimeWarning: divide by zero encountered in log10\n",
"<ipython-input-19-50856f807a2a>:20: RuntimeWarning: divide by zero encountered in log10\n",
" cm_func = lambda x: cmap(np.log10(x - min(decimations)) / (np.log10(max(decimations)) - np.log10(min(decimations))))\n"
]
}
@ -1011,9 +1153,20 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 20,
"metadata": {},
"outputs": [],
"outputs": [
{
"data": {
"text/plain": [
"41.6"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"13 * 2**5 / 10"
]

File diff suppressed because one or more lines are too long

File diff suppressed because it is too large Load diff

File diff suppressed because one or more lines are too long