From 89b32316aa46f22a4f967ca7953089a3df0c16e7 Mon Sep 17 00:00:00 2001 From: jaseg Date: Fri, 28 Feb 2020 18:29:41 +0100 Subject: embedded ldpc decoder 50% working --- controller/fw/Makefile | 5 +- controller/fw/ldpc_decoder.c | 197 +++++++++++++++++++++++++------------ controller/fw/ldpc_decoder_test.py | 53 ++++++++++ controller/fw/test_decoder.py | 167 +++++++++++++++++++++++++++++++ controller/fw/test_pyldpc_utils.py | 182 ++++++++++++++++++++++++++++++++++ 5 files changed, 542 insertions(+), 62 deletions(-) create mode 100644 controller/fw/ldpc_decoder_test.py create mode 100644 controller/fw/test_decoder.py create mode 100644 controller/fw/test_pyldpc_utils.py (limited to 'controller') diff --git a/controller/fw/Makefile b/controller/fw/Makefile index e2f0f4c..291f840 100644 --- a/controller/fw/Makefile +++ b/controller/fw/Makefile @@ -50,7 +50,7 @@ CFLAGS += -nostdlib -ffreestanding CFLAGS += -mthumb -mcpu=cortex-m4 -mfloat-abi=hard -mfpu=fpv4-sp-d16 CFLAGS += -fno-common -ffunction-sections -fdata-sections -INT_CFLAGS += -Wall -Wextra -Wshadow -Wimplicit-function-declaration -Wundef +INT_CFLAGS += -Wall -Wextra -Wpedantic -Wshadow -Wimplicit-function-declaration -Wundef INT_CFLAGS += -Wredundant-decls -Wmissing-prototypes -Wstrict-prototypes CXXFLAGS += -Os -g @@ -97,6 +97,9 @@ $(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 + gcc -fPIC -shared -Wall -Wextra -Wpedantic -std=gnu11 -O0 -g -o $@ $^ + clean: -rm -r $(BUILDDIR)/src -rm $(BUILDDIR)/$(BINARY) diff --git a/controller/fw/ldpc_decoder.c b/controller/fw/ldpc_decoder.c index 5200e8f..10b5a08 100644 --- a/controller/fw/ldpc_decoder.c +++ b/controller/fw/ldpc_decoder.c @@ -1,86 +1,160 @@ -"""Decoding module.""" - -/* H: decoding matrix as shape (m, n) array - * m: number of equations represented by H - * n: number of code words - * y: received message in codeword space, length n - * - * out: output array of length n where the solutions in codeword space will be written - * - * maxiter: Maximum iteration number - */ -int decode(uint8_t H[], const int m, n, uint8_t y[], uint8_t out[], int maxiter) { - - /* fixme: preprocess the following. -def bitsandnodes(H): - """Return bits and nodes of a parity-check matrix H.""" - 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 - bits_hist, bits_values, nodes_hist, nodes_values = utils.bitsandnodes(H) - */ +#include +#include +#include +#include +#include + +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]; - for (int n_iter=0; n_iter 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 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 new file mode 100644 index 0000000..6b14532 --- /dev/null +++ b/controller/fw/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) -- cgit