LDPC decoder fully working

1 files changed, 67 insertions, 27 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)
-*/