summaryrefslogtreecommitdiff
path: root/controller/fw/src/ldpc_decoder.c
diff options
context:
space:
mode:
authorjaseg <git-bigdata-wsl-arch@jaseg.de>2020-03-02 19:42:36 +0100
committerjaseg <git-bigdata-wsl-arch@jaseg.de>2020-03-02 19:42:36 +0100
commitca01d52a8658ed01301f14a4d294097d4db533cf (patch)
tree28b32499b6a849beda66c6ebc70239fbd9f21530 /controller/fw/src/ldpc_decoder.c
parent5effadcbaf66f476c8fffefa2c349676f41c3f52 (diff)
downloadmaster-thesis-ca01d52a8658ed01301f14a4d294097d4db533cf.tar.gz
master-thesis-ca01d52a8658ed01301f14a4d294097d4db533cf.tar.bz2
master-thesis-ca01d52a8658ed01301f14a4d294097d4db533cf.zip
Finishing up freq meas
Diffstat (limited to 'controller/fw/src/ldpc_decoder.c')
-rw-r--r--controller/fw/src/ldpc_decoder.c267
1 files changed, 267 insertions, 0 deletions
diff --git a/controller/fw/src/ldpc_decoder.c b/controller/fw/src/ldpc_decoder.c
new file mode 100644
index 0000000..fe59d77
--- /dev/null
+++ b/controller/fw/src/ldpc_decoder.c
@@ -0,0 +1,267 @@
+
+#include <stdint.h>
+#include <unistd.h>
+#include <stdbool.h>
+#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[],
+ 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];
+
+ /* 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<bits_count; j++)
+ sum ^= out[bits[i*bits_count + j]];
+ if (sum)
+ continue;
+ }
+
+ fflush(stdout);
+ return n_iter;
+ }
+
+ fflush(stdout);
+ return -1;
+}
+
+/* Perform inner ext LogBP solver */
+void inner_logbp(
+ 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 */
+ unsigned int bits_counter = 0;
+ for (size_t i=0; i<m; 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);
+
+ 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);
+ x *= tanhf(0.5f * Lc[bits_values[q]]);
+ }
+ }
+
+ } else {
+ for (size_t q=bits_counter; q<bits_counter+bits_count; q++) {
+ if (bits_values[q] != j)
+ 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)
+ Lr[i*n + j] = -1.0f;
+ else if (denom == 0)
+ Lr[i*n + j] = 1.0f;
+ 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++) {
+ for (size_t p=bits_counter; p<nodes_counter+nodes_count; p++) {
+ size_t i = nodes_values[p];
+
+ Lq[i*n + j] = Lc[j];
+
+ for (size_t q=bits_counter; q<nodes_counter+nodes_count; q++) {
+ if (nodes_values[q] != i)
+ Lq[i*n + j] += Lr[nodes_values[q]*n + j];
+ }
+ }
+
+ nodes_counter += nodes_count;
+ }
+
+ /* LLR a posteriori */
+ nodes_counter = 0;
+ for (size_t j=0; j<n; j++) {
+ float sum = 0;
+ for (size_t k=bits_counter; k<nodes_counter+nodes_count; k++)
+ sum += Lr[nodes_values[k]*n + j];
+ nodes_counter += nodes_count;
+
+ L_posteriori_out[j] = Lc[j] + sum;
+ }
+}
+
+/* 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) {
+
+ gausselimination(n, k, tG, x);
+
+ out[k - 1] = x[k - 1];
+ for (ssize_t i=k-2; i>=0; i--) {
+ out[i] = x[i];
+
+ uint8_t sum = 0;
+ for (size_t j=i+1; j<k; j++)
+ sum ^= tG[i*k + j] * out[j];
+
+ out[i] = !!(out[i] - sum);
+ }
+}
+
+/* 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;
+ }
+
+ 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]);
+ }
+ }
+ }
+}
+