diff options
Diffstat (limited to 'controller/fw')
| -rw-r--r-- | controller/fw/ldpc_decoder.c | 94 | ||||
| -rw-r--r-- | controller/fw/ldpc_decoder_test.py | 35 | ||||
| -rw-r--r-- | controller/fw/test_decoder.py | 27 |
3 files changed, 108 insertions, 48 deletions
diff --git a/controller/fw/ldpc_decoder.c b/controller/fw/ldpc_decoder.c index 10b5a08..fe59d77 100644 --- a/controller/fw/ldpc_decoder.c +++ b/controller/fw/ldpc_decoder.c @@ -5,6 +5,9 @@ #include <math.h> #include <stdio.h> + +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[], @@ -60,6 +63,7 @@ int decode(size_t n, size_t nodes_count, size_t bits_count, uint32_t bits[], int 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}; @@ -86,6 +90,7 @@ int decode(size_t n, size_t nodes_count, size_t bits_count, uint32_t bits[], int } printf("\n]\n"); } + */ for (size_t i=0; i<n; i++) out[i] = L_posteriori[i] <= 0.0f; @@ -130,18 +135,18 @@ void inner_logbp( /* step 1 : Horizontal */ unsigned int bits_counter = 0; for (size_t i=0; i<m; i++) { - printf("=== i=%zu\n", i); + //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); + //printf("\033[38;5;240mj=%04zd ", j); float x = 1; if (n_iter == 0) { 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); + //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]]); } } @@ -153,7 +158,7 @@ void inner_logbp( } } - printf("\n==== i=%03zd p=%01zd x=%08f\n", i, p-bits_counter, x); + //printf("\n==== i=%03zd p=%01zd x=%08f\n", i, p-bits_counter, x); float num = 1 + x; float denom = 1 - x; @@ -197,31 +202,66 @@ void inner_logbp( } } -/* TODO -def get_message(tG, x): - """Compute the original `n_bits` message from a `n_code` codeword `x`. +/* Compute the original (k) bit message from a (n) bit codeword x. + * + * tG: (n, k)-matrix + * x: (n)-vector + * out: (k)-vector + */ +void get_message(size_t n, size_t k, int8_t *tG, int8_t *x, int8_t *out) { - Parameters - ---------- - tG: array (n_code, n_bits) coding matrix tG. - x: array (n_code,) decoded codeword of length `n_code`. + gausselimination(n, k, tG, x); - Returns - ------- - message: array (n_bits,). Original binary message. + out[k - 1] = x[k - 1]; + for (ssize_t i=k-2; i>=0; i--) { + out[i] = x[i]; - """ - n, k = tG.shape + uint8_t sum = 0; + for (size_t j=i+1; j<k; j++) + sum ^= tG[i*k + j] * out[j]; - rtG, rx = utils.gausselimination(tG, x) + out[i] = !!(out[i] - sum); + } +} - message = np.zeros(k).astype(int) +/* Solve linear system in Z/2Z via Gauss Gauss elimination. + * + * A: (n, k)-matrix + * b: (n)-vector + */ +void gausselimination(size_t n, size_t k, int8_t *A, int8_t *b) { + ssize_t d = k<n ? k : n; + for (ssize_t j=0; j<d; j++) { + + ssize_t pivot = -1; + for (size_t i=j; i<n; i++) { + if (A[i*k + j]) { + pivot = i; + break; + } + } + if (pivot == -1) + continue; + + if (pivot != j) { + for (size_t i=0; i<k; i++) { + int8_t tmp = A[j*k + i]; + A[j*k + i] = A[pivot*k + i]; + A[pivot*k + i] = tmp; + } - 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))]) + int8_t tmp = b[j]; + b[j] = b[pivot]; + b[pivot] = tmp; + } + + for (size_t i=j+1; i<n; i++) { + if (A[i*k + j]) { + for (size_t p=0; p<k; p++) + A[i*k + p] = !!(A[i*k + p] - A[j*k + p]); + b[i] = !!(b[i] - b[j]); + } + } + } +} - return abs(message) -*/ diff --git a/controller/fw/ldpc_decoder_test.py b/controller/fw/ldpc_decoder_test.py index d5cbb24..3b91bba 100644 --- a/controller/fw/ldpc_decoder_test.py +++ b/controller/fw/ldpc_decoder_test.py @@ -16,8 +16,8 @@ args = parser.parse_args() lib = C.CDLL('./ldpc_decoder_test.so') -n = 384 -nodes, bits = 6, 8 +n = 5*19 +nodes, bits = 17, 19 H, G = pyldpc.make_ldpc(n, nodes, bits, systematic=False, seed=0) _1, bits_pos, _2 = scipy.sparse.find(H) _, k = G.shape @@ -26,18 +26,23 @@ st = np.random.RandomState(seed=0) test_data = st.randint(0, 2, k) d = np.dot(G, test_data) % 2 x = (-1) ** d +x[29:] = 0 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] +lib.get_message.argtypes = [C.c_size_t, C.c_size_t, C.POINTER(C.c_int8), C.POINTER(C.c_int8), C.POINTER(C.c_int8)] 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) + print('msg reconstruction:', test_decoder.get_message(G, ref_out)) + print('reference decoder: ', np.all(np.equal(test_decoder.get_message(G, ref_out), test_data)), flush=True) + np.set_printoptions(linewidth=220) + print(test_data) + print(test_decoder.get_message(G, ref_out)) + print(test_decoder.get_message(G, ref_out) ^ test_data) else: out = np.zeros(n, dtype=np.uint8) @@ -48,6 +53,20 @@ else: 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) + print('msg reconstruction:', test_decoder.get_message(G, out.astype(np.int64))) + print('decoder under test:', np.all(np.equal(test_decoder.get_message(G, out.astype(np.int64)), test_data))) + np.set_printoptions(linewidth=220) + print(test_data) + print(test_decoder.get_message(G, out.astype(np.int64))) + G = G.astype(np.int8) + msg = np.zeros(k, dtype=np.int8) + lib.get_message( + n, k, + G.ctypes.data_as(C.POINTER(C.c_int8)), + out.astype(np.int8).ctypes.data_as(C.POINTER(C.c_int8)), + msg.ctypes.data_as(C.POINTER(C.c_int8))) + print(msg) + print(msg ^ test_data) + +print('codeword length:', len(x)) +print('data length:', len(test_data)) diff --git a/controller/fw/test_decoder.py b/controller/fw/test_decoder.py index 806125b..8be5b02 100644 --- a/controller/fw/test_decoder.py +++ b/controller/fw/test_decoder.py @@ -43,23 +43,24 @@ def decode(H, y, snr, maxiter=100): Lr = np.zeros(shape=(m, n, n_messages)) for n_iter in range(maxiter): - print(f'============================ iteration {n_iter} ============================') + #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("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() + #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: @@ -83,13 +84,13 @@ def _logbp_numba(bits_hist, bits_values, nodes_hist, nodes_values, Lc, Lq, Lr, bits_counter = 0 nodes_counter = 0 for i in range(m): - print(f'=== i={i}') + #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=' ') + #print(f'\033[38;5;240mj={j:04d}', end=' ') X = np.ones(n_messages) if n_iter == 0: @@ -97,13 +98,13 @@ def _logbp_numba(bits_hist, bits_values, nodes_hist, nodes_values, Lc, Lq, Lr, 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=' ') + #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}') + #print(f'\n==== {i:03d} {j_iter:01d} {X[0]:> 8f}') num = 1 + X denom = 1 - X for ll in range(n_messages): |