From 5a91daf0439d30914056c3f45c9e09b7e58109a2 Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Fri, 16 Sep 2016 03:09:06 +0000
Subject: [PATCH] refactored LDPC decoder tomake it easier to call from other C
 programs

git-svn-id: https://svn.code.sf.net/p/freetel/code@2865 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2-dev/src/drs232.c        |  10 +--
 codec2-dev/src/ldpc_dec.c      | 143 +++++----------------------------
 codec2-dev/src/mpdecode_core.c | 133 +++++++++++++++++++++++++++++-
 codec2-dev/src/mpdecode_core.h |  21 ++++-
 4 files changed, 174 insertions(+), 133 deletions(-)

diff --git a/codec2-dev/src/drs232.c b/codec2-dev/src/drs232.c
index f0a08ef1..0bf89701 100644
--- a/codec2-dev/src/drs232.c
+++ b/codec2-dev/src/drs232.c
@@ -60,7 +60,7 @@
 #define BITS_PER_BYTE          10
 #define UNPACKED_PACKET_BYTES  ((UW_BYTES+BYTES_PER_PACKET+CRC_BYTES)*BITS_PER_BYTE)
 
-/* UW pattern we look for, inclduing start/stop bits */
+/* UW pattern we look for, including start/stop bits */
 
 uint8_t uw[] = {
     /* 0xb                0xa */
@@ -95,20 +95,20 @@ int main(int argc, char *argv[]) {
     uint16_t    tx_checksum, rx_checksum;
 
     if (argc < 3) {
-	printf("usage: drs232 InputOneBitPerChar OutputPackets\n");
-	exit(1);
+	fprintf(stderr, "usage: drs232 InputOneBitPerChar OutputPackets [--sd]\n");
+ 	exit(1);
     }
 
     if (strcmp(argv[1], "-")  == 0) fin = stdin;
     else if ( (fin = fopen(argv[1],"rb")) == NULL ) {
-	fprintf(stderr, "Error opening input speech file: %s: %s.\n",
+	fprintf(stderr, "Error opening input file: %s: %s.\n",
          argv[2], strerror(errno));
 	exit(1);
     }
 
     if (strcmp(argv[2], "-") == 0) fout = stdout;
     else if ( (fout = fopen(argv[2],"wb")) == NULL ) {
-	fprintf(stderr, "Error opening output speech file: %s: %s.\n",
+	fprintf(stderr, "Error opening output file: %s: %s.\n",
          argv[3], strerror(errno));
 	exit(1);
     }
diff --git a/codec2-dev/src/ldpc_dec.c b/codec2-dev/src/ldpc_dec.c
index ec1c5286..a60323aa 100644
--- a/codec2-dev/src/ldpc_dec.c
+++ b/codec2-dev/src/ldpc_dec.c
@@ -50,13 +50,12 @@
 
 #include "ldpc_code.h"  
 
-void run_ldpc_decoder(int DecodedBits[], int ParityCheckCount[], double input[]);
-
 int main(int argc, char *argv[])
 {    
     int         CodeLength, NumberParityBits, max_iter;
     int         i, j, r, num_ok, num_runs;
     char        out_char[CODELENGTH];
+    struct LDPC ldpc;
 
     /* derive some parameters */
 
@@ -75,6 +74,20 @@ int main(int argc, char *argv[])
         exit(0);
     }
 
+    /* set up LDPC code from include file constants */
+
+    ldpc.max_iter = MAX_ITER;
+    ldpc.dec_type = 0;
+    ldpc.q_scale_factor = 1;
+    ldpc.r_scale_factor = 1;
+    ldpc.CodeLength = CODELENGTH;
+    ldpc.NumberParityBits = NUMBERPARITYBITS;
+    ldpc.NumberRowsHcols = NUMBERROWSHCOLS;
+    ldpc.max_row_weight = MAX_ROW_WEIGHT;
+    ldpc.max_col_weight = MAX_COL_WEIGHT;
+    ldpc.H_rows = H_rows;
+    ldpc.H_cols = H_cols;
+
     int *DecodedBits = calloc( max_iter*CodeLength, sizeof( int ) );
     int *ParityCheckCount = calloc( max_iter, sizeof(int) );
 
@@ -90,7 +103,7 @@ int main(int argc, char *argv[])
 
         for(r=0; r<num_runs; r++) {
 
-            run_ldpc_decoder(DecodedBits, ParityCheckCount, input);
+            run_ldpc_decoder(&ldpc, DecodedBits, ParityCheckCount, input);
 
             /* Check output by comparing every output iteration */
 
@@ -137,41 +150,13 @@ int main(int argc, char *argv[])
                 sdinput = 1;
 
         double *input_double = calloc(CodeLength, sizeof(double));
-        double sum, mean, sign, sumsq, estvar, estEsN0, x;
 
         while(fread(input_double, sizeof(double), CodeLength, fin) == CodeLength) {
             if (sdinput) {
-                /* convert SD samples to LLRs -------------------------------*/
-
-                sum = 0.0;
-                for(i=0; i<CodeLength; i++)
-                    sum += fabs(input_double[i]);
-                mean = sum/CodeLength;
-                
-                /* scale by mean to map onto +/- 1 symbol position */
-
-                for(i=0; i<CodeLength; i++) {
-                    input_double[i] /= mean;
-                }
-
-                /* find variance from +/-1 symbol position */
-
-                sum = sumsq = 0.0; 
-                for(i=0; i<CodeLength; i++) {
-                    sign = (input_double[i] > 0.0) - (input_double[i] < 0.0);
-                    x = (input_double[i] - sign);
-                    sum += x;
-                    sumsq += x*x;
-                }
-                mean = sum/CodeLength;
-                estvar = sumsq/CodeLength - mean*mean;
-
-                estEsN0 = 1.0/(2.0 * estvar + 1E-3); 
-                for(i=0; i<CodeLength; i++)
-                    input_double[i] = 4.0 * estEsN0 * input_double[i];              
+                sd_to_llr(input_double, input_double, CodeLength);
             }
 
-            run_ldpc_decoder(DecodedBits, ParityCheckCount, input_double);
+            run_ldpc_decoder(&ldpc, DecodedBits, ParityCheckCount, input_double);
             
             /* extract output bits from ouput iteration that solved all parity equations, or failing that
                the last iteration. */
@@ -209,97 +194,5 @@ int main(int argc, char *argv[])
 }
 
 
-void run_ldpc_decoder(int DecodedBits[], int ParityCheckCount[], double input[]) {
-    int		max_iter, dec_type;
-    float       q_scale_factor, r_scale_factor;
-    int		max_row_weight, max_col_weight;
-    int         CodeLength, NumberParityBits, NumberRowsHcols, shift, H1;
-    int         i;
-    struct c_node *c_nodes;
-    struct v_node *v_nodes;
-    
-    /* default values */
-
-    max_iter  = MAX_ITER;
-    dec_type  = 0;
-    q_scale_factor = 1;
-    r_scale_factor = 1;
-
-    /* derive some parameters */
-
-    CodeLength = CODELENGTH;                    /* length of entire codeword */
-    NumberParityBits = NUMBERPARITYBITS;
-    NumberRowsHcols = NUMBERROWSHCOLS;
-
-    shift = (NumberParityBits + NumberRowsHcols) - CodeLength;
-    if (NumberRowsHcols == CodeLength) {
-        H1=0;
-        shift=0;
-    } else {
-        H1=1;
-    }
-	
-    max_row_weight = MAX_ROW_WEIGHT;
-    max_col_weight = MAX_COL_WEIGHT;
-	
-    c_nodes = calloc( NumberParityBits, sizeof( struct c_node ) );
-    v_nodes = calloc( CodeLength, sizeof( struct v_node));
-	
-    /* initialize c-node and v-node structures */
-
-    c_nodes = calloc( NumberParityBits, sizeof( struct c_node ) );
-    v_nodes = calloc( CodeLength, sizeof( struct v_node));
-	
-    init_c_v_nodes(c_nodes, shift, NumberParityBits, max_row_weight, H_rows, H1, CodeLength, 
-                   v_nodes, NumberRowsHcols, H_cols, max_col_weight, dec_type, input);
-
-    int DataLength = CodeLength - NumberParityBits;
-    int *data_int = calloc( DataLength, sizeof(int) );
-	
-    /* need to clear these on each call */
-
-    for(i=0; i<max_iter; i++)
-        ParityCheckCount[i] = 0;
-     for(i=0; i<max_iter*CodeLength; i++)
-         DecodedBits[i] = 0;
-
-    /* Call function to do the actual decoding */
-
-    if ( dec_type == 1) {
-        MinSum( ParityCheckCount, DecodedBits, c_nodes, v_nodes, CodeLength, 
-                NumberParityBits, max_iter, r_scale_factor, q_scale_factor, data_int );
-    } else if ( dec_type == 2) {
-        fprintf(stderr, "dec_type = 2 not currently supported");
-        /* ApproximateMinStar( BitErrors, DecodedBits, c_nodes, v_nodes, 
-           CodeLength, NumberParityBits, max_iter, r_scale_factor, q_scale_factor );*/
-    } else {
-        SumProduct( ParityCheckCount, DecodedBits, c_nodes, v_nodes, CodeLength, 
-                    NumberParityBits, max_iter, r_scale_factor, q_scale_factor, data_int ); 
-    }
 
-    /* Clean up memory */
 
-    free( data_int );
-
-    /*  Cleaning c-node elements */
-
-    for (i=0;i<NumberParityBits;i++) {
-        free( c_nodes[i].index );
-        free( c_nodes[i].message );
-        free( c_nodes[i].socket );
-    }
-	
-    /* printf( "Cleaning c-nodes \n" ); */
-    free( c_nodes );
-	
-    /* printf( "Cleaning v-node elements\n" ); */
-    for (i=0;i<CodeLength;i++) {
-        free( v_nodes[i].index);
-        free( v_nodes[i].sign );
-        free( v_nodes[i].message );
-        free( v_nodes[i].socket );
-    }
-	
-    /* printf( "Cleaning v-nodes \n" ); */
-    free( v_nodes );
-}
diff --git a/codec2-dev/src/mpdecode_core.c b/codec2-dev/src/mpdecode_core.c
index b308fd91..746bc5ff 100644
--- a/codec2-dev/src/mpdecode_core.c
+++ b/codec2-dev/src/mpdecode_core.c
@@ -9,6 +9,7 @@
 
 #include <math.h>
 #include <stdlib.h>
+#include <stdio.h>
 #include "mpdecode_core.h"
 
 /* Phi function */
@@ -553,9 +554,137 @@ void SumProduct(	 int	  BitErrors[],
   }
    
   // printf(" ssum is %d \n",   ssum); 
-   
- 
+}
+
+
+/* Convenience function to call LDPC decoder from C programs */
+
+void run_ldpc_decoder(struct LDPC *ldpc, int DecodedBits[], int ParityCheckCount[], double input[]) {
+    int		max_iter, dec_type;
+    float       q_scale_factor, r_scale_factor;
+    int		max_row_weight, max_col_weight;
+    int         CodeLength, NumberParityBits, NumberRowsHcols, shift, H1;
+    int         i;
+    struct c_node *c_nodes;
+    struct v_node *v_nodes;
     
+    /* default values */
+
+    max_iter  = ldpc->max_iter;
+    dec_type  = ldpc->dec_type;
+    q_scale_factor = ldpc->q_scale_factor;
+    r_scale_factor = ldpc->r_scale_factor;
+
+    CodeLength = ldpc->CodeLength;                    /* length of entire codeword */
+    NumberParityBits = ldpc->NumberParityBits;
+    NumberRowsHcols = ldpc->NumberRowsHcols;
+
+    /* derive some parameters */
+
+    shift = (NumberParityBits + NumberRowsHcols) - CodeLength;
+    if (NumberRowsHcols == CodeLength) {
+        H1=0;
+        shift=0;
+    } else {
+        H1=1;
+    }
+	
+    max_row_weight = ldpc->max_row_weight;
+    max_col_weight = ldpc->max_col_weight;
+	
+    c_nodes = calloc( NumberParityBits, sizeof( struct c_node ) );
+    v_nodes = calloc( CodeLength, sizeof( struct v_node));
+	
+    /* initialize c-node and v-node structures */
+
+    c_nodes = calloc( NumberParityBits, sizeof( struct c_node ) );
+    v_nodes = calloc( CodeLength, sizeof( struct v_node));
+	
+    init_c_v_nodes(c_nodes, shift, NumberParityBits, max_row_weight, ldpc->H_rows, H1, CodeLength, 
+                   v_nodes, NumberRowsHcols, ldpc->H_cols, max_col_weight, dec_type, input);
+
+    int DataLength = CodeLength - NumberParityBits;
+    int *data_int = calloc( DataLength, sizeof(int) );
+	
+    /* need to clear these on each call */
+
+    for(i=0; i<max_iter; i++)
+        ParityCheckCount[i] = 0;
+     for(i=0; i<max_iter*CodeLength; i++)
+         DecodedBits[i] = 0;
+
+    /* Call function to do the actual decoding */
+
+    if ( dec_type == 1) {
+        MinSum( ParityCheckCount, DecodedBits, c_nodes, v_nodes, CodeLength, 
+                NumberParityBits, max_iter, r_scale_factor, q_scale_factor, data_int );
+    } else if ( dec_type == 2) {
+        fprintf(stderr, "dec_type = 2 not currently supported");
+        /* ApproximateMinStar( BitErrors, DecodedBits, c_nodes, v_nodes, 
+           CodeLength, NumberParityBits, max_iter, r_scale_factor, q_scale_factor );*/
+    } else {
+        SumProduct( ParityCheckCount, DecodedBits, c_nodes, v_nodes, CodeLength, 
+                    NumberParityBits, max_iter, r_scale_factor, q_scale_factor, data_int ); 
+    }
+
+    /* Clean up memory */
+
+    free( data_int );
+
+    /*  Cleaning c-node elements */
+
+    for (i=0;i<NumberParityBits;i++) {
+        free( c_nodes[i].index );
+        free( c_nodes[i].message );
+        free( c_nodes[i].socket );
+    }
+	
+    /* printf( "Cleaning c-nodes \n" ); */
+    free( c_nodes );
+	
+    /* printf( "Cleaning v-node elements\n" ); */
+    for (i=0;i<CodeLength;i++) {
+        free( v_nodes[i].index);
+        free( v_nodes[i].sign );
+        free( v_nodes[i].message );
+        free( v_nodes[i].socket );
+    }
+	
+    /* printf( "Cleaning v-nodes \n" ); */
+    free( v_nodes );
 }
 
 
+void sd_to_llr(double llr[], double sd[], int n) {
+    double sum, mean, sign, sumsq, estvar, estEsN0, x;
+    int i;
+
+    /* convert SD samples to LLRs -------------------------------*/
+
+    sum = 0.0;
+    for(i=0; i<n; i++)
+        sum += fabs(sd[i]);
+    mean = sum/n;
+                
+    /* scale by mean to map onto +/- 1 symbol position */
+
+    for(i=0; i<n; i++) {
+        sd[i] /= mean;
+    }
+
+    /* find variance from +/-1 symbol position */
+
+    sum = sumsq = 0.0; 
+    for(i=0; i<n; i++) {
+        sign = (sd[i] > 0.0) - (sd[i] < 0.0);
+        x = (sd[i] - sign);
+        sum += x;
+        sumsq += x*x;
+    }
+    mean = sum/n;
+    estvar = sumsq/n - mean*mean;
+
+    estEsN0 = 1.0/(2.0 * estvar + 1E-3); 
+    for(i=0; i<n; i++)
+        llr[i] = 4.0 * estEsN0 * sd[i];              
+}
diff --git a/codec2-dev/src/mpdecode_core.h b/codec2-dev/src/mpdecode_core.h
index be430230..44a6b886 100644
--- a/codec2-dev/src/mpdecode_core.h
+++ b/codec2-dev/src/mpdecode_core.h
@@ -4,12 +4,31 @@
   CREATED: Sep 2016
 
   C-callable core functions for MpDecode, so they can be used for
-  Octave and C programs.
+  Octave and C programs.  Also some convenience functions to help use
+  the C-callable LDPC decoder in C programs.
 */
 
 #ifndef __MPDECODE_CORE__
 #define __MPDECODE_CORE__
 
+struct LDPC {
+    int max_iter;
+    int dec_type;
+    int q_scale_factor;
+    int r_scale_factor;
+    int CodeLength;
+    int NumberParityBits;
+    int NumberRowsHcols;
+    int max_row_weight;
+    int max_col_weight;
+    double *H_rows;
+    double *H_cols;
+};
+
+void run_ldpc_decoder(struct LDPC *ldpc, int DecodedBits[], int ParityCheckCount[], double input[]);
+
+void sd_to_llr(double llr[], double sd[], int n);
+
 struct v_node {
   int degree;
   float initial_value;
-- 
2.25.1