From ca01d52a8658ed01301f14a4d294097d4db533cf Mon Sep 17 00:00:00 2001 From: jaseg Date: Mon, 2 Mar 2020 19:42:36 +0100 Subject: Finishing up freq meas --- controller/fw/Makefile | 46 ++++- controller/fw/ldpc_decoder.c | 267 --------------------------- controller/fw/ldpc_decoder_test.py | 72 -------- controller/fw/levmarq | 1 + controller/fw/src/freq_meas.c | 94 ++++++++++ controller/fw/src/freq_meas.h | 7 + controller/fw/src/gold_code.h | 4 + controller/fw/src/ldpc_decoder.c | 267 +++++++++++++++++++++++++++ controller/fw/src/ldpc_decoder_test.py | 72 ++++++++ controller/fw/src/test_decoder.py | 168 +++++++++++++++++ controller/fw/src/test_pyldpc_utils.py | 182 ++++++++++++++++++ controller/fw/test_decoder.py | 168 ----------------- controller/fw/test_pyldpc_utils.py | 182 ------------------ controller/fw/tools/fft_window_header_gen.py | 59 ++++++ controller/fw/tools/gold_code_header_gen.py | 45 +++++ 15 files changed, 940 insertions(+), 694 deletions(-) delete mode 100644 controller/fw/ldpc_decoder.c delete mode 100644 controller/fw/ldpc_decoder_test.py create mode 160000 controller/fw/levmarq create mode 100644 controller/fw/src/freq_meas.c create mode 100644 controller/fw/src/freq_meas.h create mode 100644 controller/fw/src/gold_code.h create mode 100644 controller/fw/src/ldpc_decoder.c create mode 100644 controller/fw/src/ldpc_decoder_test.py create mode 100644 controller/fw/src/test_decoder.py create mode 100644 controller/fw/src/test_pyldpc_utils.py delete mode 100644 controller/fw/test_decoder.py delete mode 100644 controller/fw/test_pyldpc_utils.py create mode 100644 controller/fw/tools/fft_window_header_gen.py create mode 100644 controller/fw/tools/gold_code_header_gen.py (limited to 'controller') diff --git a/controller/fw/Makefile b/controller/fw/Makefile index 291f840..c01b01c 100644 --- a/controller/fw/Makefile +++ b/controller/fw/Makefile @@ -6,13 +6,17 @@ LIBSODIUM_DIR ?= libsodium TINYAES_DIR ?= tinyaes MUSL_DIR ?= musl -C_SOURCES := src/main.c src/mspdebug_wrapper.c src/spi_flash.c +C_SOURCES := src/main.c src/mspdebug_wrapper.c src/spi_flash.c src/freq_meas.c C_SOURCES += $(MSPDEBUG_DIR)/drivers/jtaglib.c C_SOURCES += $(CMSIS_DIR)/CMSIS/DSP/Source/TransformFunctions/arm_rfft_f32.c C_SOURCES += $(CMSIS_DIR)/CMSIS/DSP/Source/TransformFunctions/arm_bitreversal.c C_SOURCES += $(CMSIS_DIR)/CMSIS/DSP/Source/TransformFunctions/arm_cfft_radix4_f32.c C_SOURCES += $(MUSL_DIR)/src/math/tanhf.c $(MUSL_DIR)/src/math/atanhf.c C_SOURCES += $(MUSL_DIR)/src/math/expm1f.c $(MUSL_DIR)/src/math/log1pf.c +C_SOURCES += $(MUSL_DIR)/src/math/expf.c $(MUSL_DIR)/src/math/exp2f_data.c +C_SOURCES += $(MUSL_DIR)/src/math/__math_oflowf.c +C_SOURCES += $(MUSL_DIR)/src/math/__math_uflowf.c +C_SOURCES += $(MUSL_DIR)/src/math/__math_xflowf.c CXX_SOURCES += src/ldpc_wrapper.cpp @@ -24,6 +28,14 @@ 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 ?= 8.0 + + CC := $(PREFIX)gcc CXX := $(PREFIX)g++ LD := $(PREFIX)gcc @@ -32,6 +44,7 @@ AS := $(PREFIX)as OBJCOPY := $(PREFIX)objcopy OBJDUMP := $(PREFIX)objdump GDB := $(PREFIX)gdb +PYTHON3 ?= python3 OPENCM3_DIR_ABS := $(abspath $(OPENCM3_DIR)) CMSIS_DIR_ABS := $(abspath $(CMSIS_DIR)) @@ -41,14 +54,17 @@ TINYAES_DIR_ABS := $(abspath $(TINYAES_DIR)) MUSL_DIR_ABS := $(abspath $(MUSL_DIR)) CFLAGS += -I$(OPENCM3_DIR_ABS)/include -Imspdebug/util -Imspdebug/drivers -CFLAGS += -I$(CMSIS_DIR_ABS)/CMSIS/DSP/Include -I$(CMSIS_DIR_ABS)/CMSIS/Core/Include -CFLAGS += -I$(abspath musl_include_shims) +CFLAGS += -I$(CMSIS_DIR_ABS)/CMSIS/DSP/Include -I$(CMSIS_DIR_ABS)/CMSIS/Core/Include -Dhidden= +CFLAGS += -I$(abspath musl_include_shims) -Ilevmarq -D CFLAGS += -Os -std=gnu11 -g -DSTM32F4 # Note: libopencm3 requires some standard libc definitions from stdint.h and stdbool.h, so we don't pass -nostdinc here. CFLAGS += -nostdlib -ffreestanding CFLAGS += -mthumb -mcpu=cortex-m4 -mfloat-abi=hard -mfpu=fpv4-sp-d16 CFLAGS += -fno-common -ffunction-sections -fdata-sections +CFLAGS += -DGOLD_CODE_NBITS=$(GOLD_CODE_NBITS) -DFMEAS_FFT_LEN=$(FMEAS_FFT_LEN) -DFMEAS_ADC_MAX=$(FMEAS_ADC_MAX) +CFLAGS += -DFMEAS_ADC_SAMPLING_RATE=$(FMEAS_ADC_SAMPLING_RATE) -DFMEAS_FFT_WINDOW=$(FMEAS_FFT_WINDOW) +CFLAGS += -DFMEAS_FFT_WINDOW_SIGMA=$(FMEAS_FFT_WINDOW_SIGMA) INT_CFLAGS += -Wall -Wextra -Wpedantic -Wshadow -Wimplicit-function-declaration -Wundef INT_CFLAGS += -Wredundant-decls -Wmissing-prototypes -Wstrict-prototypes @@ -71,7 +87,13 @@ all: $(BUILDDIR)/$(BINARY) OBJS := $(addprefix $(BUILDDIR)/,$(C_SOURCES:.c=.o) $(CXX_SOURCES:.cpp=.o)) -$(BUILDDIR)/$(BINARY): $(OBJS) $(OPENCM3_DIR)/lib/lib$(OPENCM3_LIB).a $(BUILDDIR)/libsodium/src/libsodium/.libs/libsodium.a $(BUILDDIR)/tinyaes/aes.o +$(BUILDDIR)/$(BINARY): $(OBJS) \ + $(OPENCM3_DIR)/lib/lib$(OPENCM3_LIB).a \ + $(BUILDDIR)/libsodium/src/libsodium/.libs/libsodium.a \ + $(BUILDDIR)/tinyaes/aes.o \ + $(BUILDDIR)/levmarq/levmarq.o \ + $(BUILDDIR)/generated/gold_code_$(GOLD_CODE_NBITS).o \ + $(BUILDDIR)/generated/fmeas_fft_window.o $(LD) -T$(LDSCRIPT) $(LDFLAGS) -o $@ -Wl,-Map=$(BUILDDIR)/src/$*.map $^ $(BUILDDIR)/src/%.o: src/%.c @@ -82,6 +104,10 @@ $(BUILDDIR)/src/%.o: src/%.cpp mkdir -p $(@D) $(CXX) $(CXXFLAGS) -o $@ -c $< +$(BUILDDIR)/generated/%.o: $(BUILDDIR)/generated/%.c + mkdir -p $(@D) + $(CC) $(CFLAGS) $(INT_CFLAGS) -o $@ -c $< + $(BUILDDIR)/%.o: %.c mkdir -p $(@D) $(CC) $(CFLAGS) $(EXT_CFLAGS) -o $@ -c $< @@ -97,11 +123,21 @@ $(BUILDDIR)/tinyaes/aes.o: mkdir -p $(@D) make -C $(@D) -f $(TINYAES_DIR_ABS)/Makefile VPATH=$(TINYAES_DIR_ABS) CFLAGS="$(CFLAGS) -c" CC=$(CC) LD=$(LD) AR=$(AR) aes.o -ldpc_decoder_test.so: ldpc_decoder.c +$(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) > $@ + + +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) mrproper: clean diff --git a/controller/fw/ldpc_decoder.c b/controller/fw/ldpc_decoder.c deleted file mode 100644 index fe59d77..0000000 --- a/controller/fw/ldpc_decoder.c +++ /dev/null @@ -1,267 +0,0 @@ - -#include -#include -#include -#include -#include - - -void gausselimination(size_t n, size_t k, int8_t *A, int8_t *b); - -void inner_logbp( - size_t m, size_t n, - size_t bits_count, size_t nodes_count, const uint32_t bits_values[], const uint32_t nodes_values[], - int8_t Lc[], - float Lq[], float Lr[], - unsigned int n_iter, - float L_posteriori_out[]); - -//decode(384, 6, 8, ...) -int decode(size_t n, size_t nodes_count, size_t bits_count, uint32_t bits[], int8_t y[], int8_t out[], unsigned int maxiter) { - const size_t m = n * nodes_count / bits_count; - float Lq[m*n]; - float Lr[m*n]; - float L_posteriori[n]; - - /* Calculate column bit positions from row bit positions */ - int32_t bits_transposed[nodes_count * n]; - for (size_t i=0; i 1000 && i/w > head_tail && i/w < m*n/w-head_tail) { - if (!ellipsis) { - ellipsis = true; - printf("\n ..."); - } - continue; - } - if (i) - printf(", "); - if (i%w == 0) - printf("\n "); - float outf = arrs[j][i]; - char *s = outf < 0 ? "\033[91m" : (outf > 0 ? "\033[92m" : "\033[94m"); - printf("%s% 012.6g\033[38;5;240m", s, outf); - } - printf("\n]\n"); - } - */ - - for (size_t i=0; i 0 ? "\033[92m" : "\033[94m"); - //printf("nij=%04u Lc=%s%3d\033[38;5;240m ", bits_values[q], s, lcv); - x *= tanhf(0.5f * Lc[bits_values[q]]); - } - } - - } else { - for (size_t q=bits_counter; q=0; i--) { - out[i] = x[i]; - - uint8_t sum = 0; - for (size_t j=i+1; j +#include + +#include +#include + +#include "freq_meas.h" +#include "sr_global.h" + + +/* FTT window lookup table defined in generated/fmeas_fft_window.c */ +extern const float * const fmeas_fft_window_table; + +/* jury-rig some definitions for these functions since the ARM headers only export an over-generalized variable bin size + * variant. */ +extern arm_status arm_rfft_32_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_64_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_128_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_256_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_512_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_1024_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_2048_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_4096_fast_init_f32(arm_rfft_fast_instance_f32 * S); + +#define CONCAT(A, B, C) A ## B ## C +#define arm_rfft_init_name(nbits) CONCAT(arm_rfft_, nbits, _fast_init_f32) + +float func_gauss_grad(float *out, float *params, int x, void *userdata); +float func_gauss(float *params, int x, void *userdata); + +int adc_buf_measure_freq(uint16_t adc_buf[FMEAS_FFT_LEN], float *out) { + int rc; + float in_buf[FMEAS_FFT_LEN]; + float out_buf[FMEAS_FFT_LEN]; + for (size_t i=0; i a_max) { + a_max = out_buf[i]; + i_max = i; + } + } + + float par[3] = { + a_max, i_max, 1.0f + }; + + if (levmarq(3, ¶ms, nbins, out_buf, NULL, func_gauss, func_gauss_grad, NULL, &lmstat)) + return -1; + + *out = (params[1] + first_bin) * binsize; + + return 0; +} + +float func_gauss(float *params, int x, void *userdata) { + UNUSED(userdata); + float a = params[0]; + float mu = params[1]; + float sigma = params[2]; + return a*expf(-arm_power_f32((x-mu), 2.0f/(2.0f*(sigma*sigma)))); +} + +float func_gauss_grad(float *out, float *params, int x, void *userdata) { + UNUSED(userdata); + float a = params[0]; + float mu = params[1]; + float sigma = params[2]; + return -(x-mu) / ( sigma*sigma*sigma * 2.5066282746310002f) * a*expf(-arm_power_f32((x-mu), 2.0f/(2.0f*(sigma*sigma)))); +} diff --git a/controller/fw/src/freq_meas.h b/controller/fw/src/freq_meas.h new file mode 100644 index 0000000..1c083f8 --- /dev/null +++ b/controller/fw/src/freq_meas.h @@ -0,0 +1,7 @@ + +#ifndef __FREQ_MEAS_H__ +#define __FREQ_MEAS_H__ + +int adc_buf_measure_freq(uint16_t adc_buf[FMEAS_FFT_LEN], float *out); + +#endif /* __FREQ_MEAS_H__ */ diff --git a/controller/fw/src/gold_code.h b/controller/fw/src/gold_code.h new file mode 100644 index 0000000..739b477 --- /dev/null +++ b/controller/fw/src/gold_code.h @@ -0,0 +1,4 @@ + +/* header file for generated gold code tables */ + +extern const uint8_t * const gold_code_table; diff --git a/controller/fw/src/ldpc_decoder.c b/controller/fw/src/ldpc_decoder.c new file mode 100644 index 0000000..fe59d77 --- /dev/null +++ b/controller/fw/src/ldpc_decoder.c @@ -0,0 +1,267 @@ + +#include +#include +#include +#include +#include + + +void gausselimination(size_t n, size_t k, int8_t *A, int8_t *b); + +void inner_logbp( + size_t m, size_t n, + size_t bits_count, size_t nodes_count, const uint32_t bits_values[], const uint32_t nodes_values[], + int8_t Lc[], + float Lq[], float Lr[], + unsigned int n_iter, + float L_posteriori_out[]); + +//decode(384, 6, 8, ...) +int decode(size_t n, size_t nodes_count, size_t bits_count, uint32_t bits[], int8_t y[], int8_t out[], unsigned int maxiter) { + const size_t m = n * nodes_count / bits_count; + float Lq[m*n]; + float Lr[m*n]; + float L_posteriori[n]; + + /* Calculate column bit positions from row bit positions */ + int32_t bits_transposed[nodes_count * n]; + for (size_t i=0; i 1000 && i/w > head_tail && i/w < m*n/w-head_tail) { + if (!ellipsis) { + ellipsis = true; + printf("\n ..."); + } + continue; + } + if (i) + printf(", "); + if (i%w == 0) + printf("\n "); + float outf = arrs[j][i]; + char *s = outf < 0 ? "\033[91m" : (outf > 0 ? "\033[92m" : "\033[94m"); + printf("%s% 012.6g\033[38;5;240m", s, outf); + } + printf("\n]\n"); + } + */ + + for (size_t i=0; i 0 ? "\033[92m" : "\033[94m"); + //printf("nij=%04u Lc=%s%3d\033[38;5;240m ", bits_values[q], s, lcv); + x *= tanhf(0.5f * Lc[bits_values[q]]); + } + } + + } else { + for (size_t q=bits_counter; q=0; i--) { + out[i] = x[i]; + + uint8_t sum = 0; + for (size_t j=i+1; j 0 else '\033[94m') + #print(f'nij={nij[kk]:04d} Lc={lcc}{lcv:> 8f}\033[38;5;240m', end=' ') + X *= np.tanh(0.5 * Lc[nij[kk]]) + else: + for kk in range(len(nij)): + if nij[kk] != j: + X *= np.tanh(0.5 * Lq[i, nij[kk]]) + #print(f'\n==== {i:03d} {j_iter:01d} {X[0]:> 8f}') + num = 1 + X + denom = 1 - X + for ll in range(n_messages): + if num[ll] == 0: + Lr[i, j, ll] = -1 + elif denom[ll] == 0: + Lr[i, j, ll] = 1 + else: + Lr[i, j, ll] = np.log(num[ll] / denom[ll]) + # step 2 : Vertical + + for j in range(n): + ff = nodes_hist[j] + mj = nodes_values[bits_counter: nodes_counter + ff] + nodes_counter += ff + for i in mj: + mji = mj[:] + Lq[i, j] = Lc[j] + + for kk in range(len(mji)): + if mji[kk] != i: + Lq[i, j] += Lr[mji[kk], j] + + # LLR a posteriori: + L_posteriori = np.zeros((n, n_messages)) + nodes_counter = 0 + for j in range(n): + ff = nodes_hist[j] + mj = nodes_values[bits_counter: nodes_counter + ff] + nodes_counter += ff + L_posteriori[j] = Lc[j] + Lr[mj, j].sum(axis=0) + + return Lq, Lr, L_posteriori + + +def get_message(tG, x): + """Compute the original `n_bits` message from a `n_code` codeword `x`. + + Parameters + ---------- + tG: array (n_code, n_bits) coding matrix tG. + x: array (n_code,) decoded codeword of length `n_code`. + + Returns + ------- + message: array (n_bits,). Original binary message. + + """ + n, k = tG.shape + + rtG, rx = utils.gausselimination(tG, x) + + message = np.zeros(k).astype(int) + + message[k - 1] = rx[k - 1] + for i in reversed(range(k - 1)): + message[i] = rx[i] + message[i] -= utils.binaryproduct(rtG[i, list(range(i+1, k))], + message[list(range(i+1, k))]) + + return abs(message) diff --git a/controller/fw/src/test_pyldpc_utils.py b/controller/fw/src/test_pyldpc_utils.py new file mode 100644 index 0000000..6b14532 --- /dev/null +++ b/controller/fw/src/test_pyldpc_utils.py @@ -0,0 +1,182 @@ +"""Conversion tools.""" +import math +import numbers +import numpy as np +import scipy +from scipy.stats import norm +pi = math.pi + + +def int2bitarray(n, k): + """Change an array's base from int (base 10) to binary (base 2).""" + binary_string = bin(n) + length = len(binary_string) + bitarray = np.zeros(k, 'int') + for i in range(length - 2): + bitarray[k - i - 1] = int(binary_string[length - i - 1]) + + return bitarray + + +def bitarray2int(bitarray): + """Change array's base from binary (base 2) to int (base 10).""" + bitstring = "".join([str(i) for i in bitarray]) + + return int(bitstring, 2) + + +def binaryproduct(X, Y): + """Compute a matrix-matrix / vector product in Z/2Z.""" + A = X.dot(Y) + try: + A = A.toarray() + except AttributeError: + pass + return A % 2 + + +def gaussjordan(X, change=0): + """Compute the binary row reduced echelon form of X. + + Parameters + ---------- + X: array (m, n) + change : boolean (default, False). If True returns the inverse transform + + Returns + ------- + if `change` == 'True': + A: array (m, n). row reduced form of X. + P: tranformations applied to the identity + else: + A: array (m, n). row reduced form of X. + + """ + A = np.copy(X) + m, n = A.shape + + if change: + P = np.identity(m).astype(int) + + pivot_old = -1 + for j in range(n): + filtre_down = A[pivot_old+1:m, j] + pivot = np.argmax(filtre_down)+pivot_old+1 + + if A[pivot, j]: + pivot_old += 1 + if pivot_old != pivot: + aux = np.copy(A[pivot, :]) + A[pivot, :] = A[pivot_old, :] + A[pivot_old, :] = aux + if change: + aux = np.copy(P[pivot, :]) + P[pivot, :] = P[pivot_old, :] + P[pivot_old, :] = aux + + for i in range(m): + if i != pivot_old and A[i, j]: + if change: + P[i, :] = abs(P[i, :]-P[pivot_old, :]) + A[i, :] = abs(A[i, :]-A[pivot_old, :]) + + if pivot_old == m-1: + break + + if change: + return A, P + return A + + +def binaryrank(X): + """Compute rank of a binary Matrix using Gauss-Jordan algorithm.""" + A = np.copy(X) + m, n = A.shape + + A = gaussjordan(A) + + return sum([a.any() for a in A]) + + +def f1(y, sigma): + """Compute normal density N(1,sigma).""" + f = norm.pdf(y, loc=1, scale=sigma) + return f + + +def fm1(y, sigma): + """Compute normal density N(-1,sigma).""" + + f = norm.pdf(y, loc=-1, scale=sigma) + return f + + +def bitsandnodes(H): + """Return bits and nodes of a parity-check matrix H.""" + if type(H) != scipy.sparse.csr_matrix: + bits_indices, bits = np.where(H) + nodes_indices, nodes = np.where(H.T) + else: + bits_indices, bits = scipy.sparse.find(H)[:2] + nodes_indices, nodes = scipy.sparse.find(H.T)[:2] + bits_histogram = np.bincount(bits_indices) + nodes_histogram = np.bincount(nodes_indices) + + return bits_histogram, bits, nodes_histogram, nodes + + +def incode(H, x): + """Compute Binary Product of H and x.""" + return (binaryproduct(H, x) == 0).all() + + +def gausselimination(A, b): + """Solve linear system in Z/2Z via Gauss Gauss elimination.""" + if type(A) == scipy.sparse.csr_matrix: + A = A.toarray().copy() + else: + A = A.copy() + b = b.copy() + n, k = A.shape + + for j in range(min(k, n)): + listedepivots = [i for i in range(j, n) if A[i, j]] + if len(listedepivots): + pivot = np.min(listedepivots) + else: + continue + if pivot != j: + aux = (A[j, :]).copy() + A[j, :] = A[pivot, :] + A[pivot, :] = aux + + aux = b[j].copy() + b[j] = b[pivot] + b[pivot] = aux + + for i in range(j+1, n): + if A[i, j]: + A[i, :] = abs(A[i, :]-A[j, :]) + b[i] = abs(b[i]-b[j]) + + return A, b + + +def check_random_state(seed): + """Turn seed into a np.random.RandomState instance + Parameters + ---------- + seed : None | int | instance of RandomState + If seed is None, return the RandomState singleton used by np.random. + If seed is an int, return a new RandomState instance seeded with seed. + If seed is already a RandomState instance, return it. + Otherwise raise ValueError. + """ + if seed is None or seed is np.random: + return np.random.mtrand._rand + if isinstance(seed, numbers.Integral): + return np.random.RandomState(seed) + if isinstance(seed, np.random.RandomState): + return seed + raise ValueError('%r cannot be used to seed a numpy.random.RandomState' + ' instance' % seed) diff --git a/controller/fw/test_decoder.py b/controller/fw/test_decoder.py deleted file mode 100644 index 8be5b02..0000000 --- a/controller/fw/test_decoder.py +++ /dev/null @@ -1,168 +0,0 @@ -"""Decoding module.""" -import numpy as np -import warnings -import test_pyldpc_utils as utils - -from numba import njit, int64, types, float64 - -np.set_printoptions(linewidth=180, threshold=1000, edgeitems=20) - -def decode(H, y, snr, maxiter=100): - """Decode a Gaussian noise corrupted n bits message using BP algorithm. - - Decoding is performed in parallel if multiple codewords are passed in y. - - Parameters - ---------- - H: array (n_equations, n_code). Decoding matrix H. - y: array (n_code, n_messages) or (n_code,). Received message(s) in the - codeword space. - maxiter: int. Maximum number of iterations of the BP algorithm. - - Returns - ------- - x: array (n_code,) or (n_code, n_messages) the solutions in the - codeword space. - - """ - m, n = H.shape - - bits_hist, bits_values, nodes_hist, nodes_values = utils.bitsandnodes(H) - - var = 10 ** (-snr / 10) - - if y.ndim == 1: - y = y[:, None] - # step 0: initialization - - Lc = 2 * y / var - _, n_messages = y.shape - - Lq = np.zeros(shape=(m, n, n_messages)) - - Lr = np.zeros(shape=(m, n, n_messages)) - - for n_iter in range(maxiter): - #print(f'============================ iteration {n_iter} ============================') - Lq, Lr, L_posteriori = _logbp_numba(bits_hist, bits_values, nodes_hist, - nodes_values, Lc, Lq, Lr, n_iter) - #print("Lq=", Lq.flatten()) - #print("Lr=", Lr.flatten()) - #print("L_posteriori=", L_posteriori.flatten()) - #print('L_posteriori=[') - #for row in L_posteriori.reshape([-1, 16]): - # for val in row: - # cc = '\033[91m' if val < 0 else ('\033[92m' if val > 0 else '\033[94m') - # print(f"{cc}{val: 012.6g}\033[38;5;240m", end=', ') - # print() - x = np.array(L_posteriori <= 0).astype(int) - - product = utils.incode(H, x) - - if product: - print(f'found, n_iter={n_iter}') - break - - if n_iter == maxiter - 1: - warnings.warn("""Decoding stopped before convergence. You may want - to increase maxiter""") - return x.squeeze() - - -output_type_log2 = types.Tuple((float64[:, :, :], float64[:, :, :], - float64[:, :])) - - -#@njit(output_type_log2(int64[:], int64[:], int64[:], int64[:], float64[:, :], -# float64[:, :, :], float64[:, :, :], int64), cache=True) -def _logbp_numba(bits_hist, bits_values, nodes_hist, nodes_values, Lc, Lq, Lr, - n_iter): - """Perform inner ext LogBP solver.""" - m, n, n_messages = Lr.shape - # step 1 : Horizontal - - bits_counter = 0 - nodes_counter = 0 - for i in range(m): - #print(f'=== i={i}') - ff = bits_hist[i] - ni = bits_values[bits_counter: bits_counter + ff] - bits_counter += ff - for j_iter, j in enumerate(ni): - nij = ni[:] - #print(f'\033[38;5;240mj={j:04d}', end=' ') - - X = np.ones(n_messages) - if n_iter == 0: - for kk in range(len(nij)): - if nij[kk] != j: - lcv = Lc[nij[kk],0] - lcc = '\033[91m' if lcv < 0 else ('\033[92m' if lcv > 0 else '\033[94m') - #print(f'nij={nij[kk]:04d} Lc={lcc}{lcv:> 8f}\033[38;5;240m', end=' ') - X *= np.tanh(0.5 * Lc[nij[kk]]) - else: - for kk in range(len(nij)): - if nij[kk] != j: - X *= np.tanh(0.5 * Lq[i, nij[kk]]) - #print(f'\n==== {i:03d} {j_iter:01d} {X[0]:> 8f}') - num = 1 + X - denom = 1 - X - for ll in range(n_messages): - if num[ll] == 0: - Lr[i, j, ll] = -1 - elif denom[ll] == 0: - Lr[i, j, ll] = 1 - else: - Lr[i, j, ll] = np.log(num[ll] / denom[ll]) - # step 2 : Vertical - - for j in range(n): - ff = nodes_hist[j] - mj = nodes_values[bits_counter: nodes_counter + ff] - nodes_counter += ff - for i in mj: - mji = mj[:] - Lq[i, j] = Lc[j] - - for kk in range(len(mji)): - if mji[kk] != i: - Lq[i, j] += Lr[mji[kk], j] - - # LLR a posteriori: - L_posteriori = np.zeros((n, n_messages)) - nodes_counter = 0 - for j in range(n): - ff = nodes_hist[j] - mj = nodes_values[bits_counter: nodes_counter + ff] - nodes_counter += ff - L_posteriori[j] = Lc[j] + Lr[mj, j].sum(axis=0) - - return Lq, Lr, L_posteriori - - -def get_message(tG, x): - """Compute the original `n_bits` message from a `n_code` codeword `x`. - - Parameters - ---------- - tG: array (n_code, n_bits) coding matrix tG. - x: array (n_code,) decoded codeword of length `n_code`. - - Returns - ------- - message: array (n_bits,). Original binary message. - - """ - n, k = tG.shape - - rtG, rx = utils.gausselimination(tG, x) - - message = np.zeros(k).astype(int) - - message[k - 1] = rx[k - 1] - for i in reversed(range(k - 1)): - message[i] = rx[i] - message[i] -= utils.binaryproduct(rtG[i, list(range(i+1, k))], - message[list(range(i+1, k))]) - - return abs(message) diff --git a/controller/fw/test_pyldpc_utils.py b/controller/fw/test_pyldpc_utils.py deleted file mode 100644 index 6b14532..0000000 --- a/controller/fw/test_pyldpc_utils.py +++ /dev/null @@ -1,182 +0,0 @@ -"""Conversion tools.""" -import math -import numbers -import numpy as np -import scipy -from scipy.stats import norm -pi = math.pi - - -def int2bitarray(n, k): - """Change an array's base from int (base 10) to binary (base 2).""" - binary_string = bin(n) - length = len(binary_string) - bitarray = np.zeros(k, 'int') - for i in range(length - 2): - bitarray[k - i - 1] = int(binary_string[length - i - 1]) - - return bitarray - - -def bitarray2int(bitarray): - """Change array's base from binary (base 2) to int (base 10).""" - bitstring = "".join([str(i) for i in bitarray]) - - return int(bitstring, 2) - - -def binaryproduct(X, Y): - """Compute a matrix-matrix / vector product in Z/2Z.""" - A = X.dot(Y) - try: - A = A.toarray() - except AttributeError: - pass - return A % 2 - - -def gaussjordan(X, change=0): - """Compute the binary row reduced echelon form of X. - - Parameters - ---------- - X: array (m, n) - change : boolean (default, False). If True returns the inverse transform - - Returns - ------- - if `change` == 'True': - A: array (m, n). row reduced form of X. - P: tranformations applied to the identity - else: - A: array (m, n). row reduced form of X. - - """ - A = np.copy(X) - m, n = A.shape - - if change: - P = np.identity(m).astype(int) - - pivot_old = -1 - for j in range(n): - filtre_down = A[pivot_old+1:m, j] - pivot = np.argmax(filtre_down)+pivot_old+1 - - if A[pivot, j]: - pivot_old += 1 - if pivot_old != pivot: - aux = np.copy(A[pivot, :]) - A[pivot, :] = A[pivot_old, :] - A[pivot_old, :] = aux - if change: - aux = np.copy(P[pivot, :]) - P[pivot, :] = P[pivot_old, :] - P[pivot_old, :] = aux - - for i in range(m): - if i != pivot_old and A[i, j]: - if change: - P[i, :] = abs(P[i, :]-P[pivot_old, :]) - A[i, :] = abs(A[i, :]-A[pivot_old, :]) - - if pivot_old == m-1: - break - - if change: - return A, P - return A - - -def binaryrank(X): - """Compute rank of a binary Matrix using Gauss-Jordan algorithm.""" - A = np.copy(X) - m, n = A.shape - - A = gaussjordan(A) - - return sum([a.any() for a in A]) - - -def f1(y, sigma): - """Compute normal density N(1,sigma).""" - f = norm.pdf(y, loc=1, scale=sigma) - return f - - -def fm1(y, sigma): - """Compute normal density N(-1,sigma).""" - - f = norm.pdf(y, loc=-1, scale=sigma) - return f - - -def bitsandnodes(H): - """Return bits and nodes of a parity-check matrix H.""" - if type(H) != scipy.sparse.csr_matrix: - bits_indices, bits = np.where(H) - nodes_indices, nodes = np.where(H.T) - else: - bits_indices, bits = scipy.sparse.find(H)[:2] - nodes_indices, nodes = scipy.sparse.find(H.T)[:2] - bits_histogram = np.bincount(bits_indices) - nodes_histogram = np.bincount(nodes_indices) - - return bits_histogram, bits, nodes_histogram, nodes - - -def incode(H, x): - """Compute Binary Product of H and x.""" - return (binaryproduct(H, x) == 0).all() - - -def gausselimination(A, b): - """Solve linear system in Z/2Z via Gauss Gauss elimination.""" - if type(A) == scipy.sparse.csr_matrix: - A = A.toarray().copy() - else: - A = A.copy() - b = b.copy() - n, k = A.shape - - for j in range(min(k, n)): - listedepivots = [i for i in range(j, n) if A[i, j]] - if len(listedepivots): - pivot = np.min(listedepivots) - else: - continue - if pivot != j: - aux = (A[j, :]).copy() - A[j, :] = A[pivot, :] - A[pivot, :] = aux - - aux = b[j].copy() - b[j] = b[pivot] - b[pivot] = aux - - for i in range(j+1, n): - if A[i, j]: - A[i, :] = abs(A[i, :]-A[j, :]) - b[i] = abs(b[i]-b[j]) - - return A, b - - -def check_random_state(seed): - """Turn seed into a np.random.RandomState instance - Parameters - ---------- - seed : None | int | instance of RandomState - If seed is None, return the RandomState singleton used by np.random. - If seed is an int, return a new RandomState instance seeded with seed. - If seed is already a RandomState instance, return it. - Otherwise raise ValueError. - """ - if seed is None or seed is np.random: - return np.random.mtrand._rand - if isinstance(seed, numbers.Integral): - return np.random.RandomState(seed) - if isinstance(seed, np.random.RandomState): - return seed - raise ValueError('%r cannot be used to seed a numpy.random.RandomState' - ' instance' % seed) diff --git a/controller/fw/tools/fft_window_header_gen.py b/controller/fw/tools/fft_window_header_gen.py new file mode 100644 index 0000000..7df2ee3 --- /dev/null +++ b/controller/fw/tools/fft_window_header_gen.py @@ -0,0 +1,59 @@ +#!/usr/bin/env python3 + +import textwrap + +import scipy.signal as sig +import numpy as np + +WINDOW_TYPES = [ + 'boxcar', + 'triang', + 'blackman', + 'hamming', + 'hann', + 'bartlett', + 'flattop', + 'parzen', + 'bohman', + 'blackmanharris', + 'nuttall', + 'barthann', + 'kaiser', + 'gaussian', + 'general_gaussian', + 'slepian', + 'dpss', + 'chebwin', + 'exponential', + 'tukey', + ] + +if __name__ == '__main__': + import argparse + parser = argparse.ArgumentParser() + parser.add_argument('window', choices=WINDOW_TYPES, help='Type of window function to use') + parser.add_argument('n', type=int, help='Width of window in samples') + parser.add_argument('window_args', nargs='*', type=float, + help='''Window argument(s) if required. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html#scipy.signal.get_window for details.''') + parser.add_argument('-v', '--variable', default='fft_window_table', help='Name for alias variable pointing to generated window') + args = parser.parse_args() + + print(f'/* FTT window table for {args.n} sample {args.window} window.') + if args.window_args: + print(f' * Window arguments were: ({" ,".join(str(arg) for arg in args.window_args)})') + print(f' */') + winargs = ''.join(f'_{arg:.4g}'.replace('.', 'F') for arg in args.window_args) + varname = f'fft_{args.n}_window_{args.window}{winargs}' + print(f'const float {varname}[{args.n}] = {{') + + win = sig.get_window(args.window if not args.window_args else (args.window, *args.window_args), + Nx=args.n, fftbins=True) + par = ' '.join(f'{f:>013.8g},' 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};') + diff --git a/controller/fw/tools/gold_code_header_gen.py b/controller/fw/tools/gold_code_header_gen.py new file mode 100644 index 0000000..e2bc210 --- /dev/null +++ b/controller/fw/tools/gold_code_header_gen.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python3 + +import textwrap + +import numpy as np +import scipy.signal as sig + +# 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) + +if __name__ == '__main__': + import argparse + parser = argparse.ArgumentParser() + 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.') + args = parser.parse_args() + + print('#include ') + 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;') -- cgit