jaseg/master-thesis
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

embedded ldpc decoder 50% working

Browse source
This commit is contained in:
jaseg 2020-02-28 18:29:41 +01:00
parent 9cf1fee2e9
commit 89b32316aa
6 changed files with 923 additions and 93 deletions
Download patch file Download diff file Expand all files Collapse all files

5
controller/fw/Makefile
View file

@ -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)

193
controller/fw/ldpc_decoder.c
View file

@ -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) {
#include <stdint.h>
#include <unistd.h>
#include <stdbool.h>
#include <math.h>
#include <stdio.h>
/* 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)
*/
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<maxiter; n_iter++) {
Lq, Lr, L_posteriori = inner_logbp(bits_hist, bits_values, nodes_hist,
nodes_values, y, Lq, Lr, n_iter, L_posteriori)
/* Calculate column bit positions from row bit positions */
int32_t bits_transposed[nodes_count * n];
for (size_t i=0; i<nodes_count * n; i++)
bits_transposed[i] = -1;
for (size_t i=0; i<m; i++) {
for (size_t j=0; j<bits_count; j++) {
int32_t *base = bits_transposed + bits[i*bits_count + j] * nodes_count;
for (; *base != -1; base++)
;
*base = i;
}
}
/*
printf("Row positions: [");
for (size_t i=0; i<m*bits_count; i++) {
if (i)
printf(", ");
if (i%32 == 0)
printf("\n ");
printf("%4d", bits[i]);
}
printf("\n]\n");
printf("Column positions: [");
for (size_t i=0; i<n*nodes_count; i++) {
if (i)
printf(", ");
if (i%32 == 0)
printf("\n ");
printf("%4d", bits_transposed[i]);
}
printf("\n]\n");
*/
/* Run iterative optimization algorithm */
for (unsigned int n_iter=0; n_iter<maxiter; n_iter++) {
inner_logbp(m, n, bits_count, nodes_count, bits, (uint32_t*)bits_transposed, y, Lq, Lr, n_iter, L_posteriori);
float *arrs[3] = {Lq, Lr, L_posteriori};
const char *names[3] = {"Lq", "Lr", "L_posteriori"};
size_t lens[3] = {m*n, m*n, n};
const size_t head_tail = 10;
for (int j=0; j<3; j++) {
printf("%s=[", names[j]);
bool ellipsis = false;
const int w = 16;
for (size_t i=0; i<lens[j]; i++) {
if (lens[j] > 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<n; i++)
out[i] = L_posteriori[i] <= 0.0f;
for (size_t i=0; i<m; i++) {
bool sum = 0;
for (size_t j=0; j<n; j++)
sum ^= H[i*n + j] * out[j];
for (size_t j=0; j<bits_count; j++)
sum ^= out[bits[i*bits_count + j]];
if (sum)
continue;
}
return 1;
fflush(stdout);
return n_iter;
}
fflush(stdout);
return -1;
}
/* Perform inner ext LogBP solver */
void inner_logbp(
size_t m, n,
int bits_hist[], bits_values[], nodes_hist[], nodes_values[],
uint8_t Lc[],
float Lq[], Lr[],
int n_iter,
size_t m, size_t n,
size_t bits_count, size_t nodes_count, uint32_t const bits_values[], const uint32_t nodes_values[],
int8_t Lc[],
float Lq[], float Lr[],
unsigned int n_iter,
float L_posteriori_out[]) {
/*
printf("Input data: [");
for (size_t i=0; i<n; i++) {
if (i)
printf(", ");
if (i%32 == 0)
printf("\n ");
printf("%4d", Lc[i]);
}
printf("\n]\n");
*/
/* step 1 : Horizontal */
int bits_counter = 0;
int nodes_counter = 0;
unsigned int bits_counter = 0;
for (size_t i=0; i<m; i++) {
int ff = bits_hist[i];
bits_counter += ff;
for (size_t p=bits_counter; p<bits_counter+ff; p++) {
int j = bits_values[p];
printf("=== i=%zu\n", i);
for (size_t p=bits_counter; p<bits_counter+bits_count; p++) {
size_t j = bits_values[p];
printf("\033[38;5;240mj=%04zd ", j);
float x = 1;
if (n_iter == 0) {
for (size_t q=bits_counter; q<bits_counter+ff; q++) {
if (bits_values[q] != j)
for (size_t q=bits_counter; q<bits_counter+bits_count; q++) {
if (bits_values[q] != j) {
int lcv = Lc[bits_values[q]];
char *s = lcv < 0 ? "\033[91m" : (lcv > 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<bits_counter+ff; q++) {
for (size_t q=bits_counter; q<bits_counter+bits_count; q++) {
if (bits_values[q] != j)
x *= tanhf(0.5f * Lq[i, bits_values[q]]);
x *= tanhf(0.5f * Lq[i*n + bits_values[q]]);
}
}
printf("\n==== i=%03zd p=%01zd x=%08f\n", i, p-bits_counter, x);
float num = 1 + x;
float denom = 1 - x;
if (num == 0)
@ -90,40 +164,40 @@ void inner_logbp(
else
Lr[i*n + j] = logf(num/denom);
}
bits_counter += bits_count;
}
/* step 2 : Vertical */
unsigned int nodes_counter = 0;
for (size_t j=0; j<n; j++) {
ff = nodes_hist[j];
for (size_t p=bits_counter; p<nodes_counter+nodes_count; p++) {
size_t i = nodes_values[p];
for (size_t p=bits_counter; p<nodes_counter+ff; p++) {
int i = nodes_values[p];
Lq[i*n + j] = Lc[j];
Lq[i, j] = Lc[j]
for (size_t q=bits_counter; q<nodes_counter+ff; q++) {
for (size_t q=bits_counter; q<nodes_counter+nodes_count; q++) {
if (nodes_values[q] != i)
Lq[i, j] += Lr[nodes_values[q], j];
Lq[i*n + j] += Lr[nodes_values[q]*n + j];
}
}
nodes_counter += ff;
nodes_counter += nodes_count;
}
/* LLR a posteriori */
size_t nodes_counter = 0;
nodes_counter = 0;
for (size_t j=0; j<n; j++) {
size_t ff = nodes_hist[j];
float sum = 0;
for (size_t k=bits_counter; k<nodes_counter+ff; k++)
for (size_t k=bits_counter; k<nodes_counter+nodes_count; k++)
sum += Lr[nodes_values[k]*n + j];
nodes_counter += ff;
nodes_counter += nodes_count;
L_posteriori_out[j] = Lc[j] + sum;
}
}
/* TODO
def get_message(tG, x):
"""Compute the original `n_bits` message from a `n_code` codeword `x`.
@ -150,3 +224,4 @@ def get_message(tG, x):
message[list(range(i+1, k))])
return abs(message)
*/

53
controller/fw/ldpc_decoder_test.py Normal file
View file

@ -0,0 +1,53 @@
import pyldpc
import scipy.sparse
import numpy as np
import test_decoder
import os, sys
import ctypes as C
import argparse
if __name__ != '__main__':
raise RuntimeError("Please don't import this module, this is a command-line program.")
parser = argparse.ArgumentParser()
parser.add_argument('-r', '--reference', action='store_true', default=False, help='Run reference decoder instead of C implemention')
args = parser.parse_args()
lib = C.CDLL('./ldpc_decoder_test.so')
n = 384
nodes, bits = 6, 8
H, G = pyldpc.make_ldpc(n, nodes, bits, systematic=False, seed=0)
_1, bits_pos, _2 = scipy.sparse.find(H)
_, k = G.shape
st = np.random.RandomState(seed=0)
test_data = st.randint(0, 2, k)
d = np.dot(G, test_data) % 2
x = (-1) ** d
bits_pos = bits_pos.astype(np.uint32)
x = x.astype(np.int8)
lib.decode.argtypes = [C.c_size_t, C.c_size_t, C.c_size_t, C.POINTER(C.c_size_t), C.POINTER(C.c_int8), C.POINTER(C.c_int8), C.c_uint]
if args.reference:
ref_out = test_decoder.decode(H, x, 3)
print('decoder output:', ref_out, flush=True)
print('msg reconstruction:', pyldpc.get_message(G, ref_out))
print('reference decoder: ', np.all(np.equal(pyldpc.get_message(G, ref_out), test_data)), flush=True)
print(ref_out.dtype)
else:
out = np.zeros(n, dtype=np.uint8)
# print('python data:', x, flush=True)
print('decoder iterations:', lib.decode(n, nodes, bits,
bits_pos.ctypes.data_as(C.POINTER(C.c_ulong)),
x.ctypes.data_as(C.POINTER(C.c_int8)),
out.ctypes.data_as(C.POINTER(C.c_int8)),
25), flush=True)
print('decoder output:', out)
print('msg reconstruction:', pyldpc.get_message(G, out))
print('decoder under test:', np.all(np.equal(pyldpc.get_message(G, out.astype(np.int64)), test_data)))
print(out.dtype)

167
controller/fw/test_decoder.py Normal file
View file

@ -0,0 +1,167 @@
"""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:
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)

182
controller/fw/test_pyldpc_utils.py Normal file
View file

@ -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)