LDPC decoder fully working
This commit is contained in:
jaseg
parent
89b32316aa
commit
5effadcbaf
3 changed files with 108 additions and 48 deletions
|
|
@ -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++) {
|
||||
|
||||
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))])
|
||||
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;
|
||||
}
|
||||
|
||||
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)
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -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))
|
||||
|
|
|
|||
|
|
@ -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):
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue