From 700e1b3b95116ee156d3af9dc7c09f8347eafbe4 Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Sun, 11 Sep 2016 21:51:39 +0000
Subject: [PATCH] initial checkin of C side of LDPC FSK decoder

git-svn-id: https://svn.code.sf.net/p/freetel/code@2852 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2-dev/src/ldpc_dec.c      | 253 +++++++++++++++
 codec2-dev/src/mpdecode_core.c | 561 +++++++++++++++++++++++++++++++++
 codec2-dev/src/mpdecode_core.h |  74 +++++
 3 files changed, 888 insertions(+)
 create mode 100644 codec2-dev/src/ldpc_dec.c
 create mode 100644 codec2-dev/src/mpdecode_core.c
 create mode 100644 codec2-dev/src/mpdecode_core.h

diff --git a/codec2-dev/src/ldpc_dec.c b/codec2-dev/src/ldpc_dec.c
new file mode 100644
index 00000000..805b4794
--- /dev/null
+++ b/codec2-dev/src/ldpc_dec.c
@@ -0,0 +1,253 @@
+/* 
+  FILE...: ldpc_dec.c
+  AUTHOR.: Matthew C. Valenti, Rohit Iyer Seshadri, David Rowe
+  CREATED: Sep 2016
+
+  Command line C LDPC decoder derived from MpDecode.c in the CML
+  library.  Allows us to run the same decoder in Octave and C.  The
+  code is defined by the parameters and array stored in the include
+  file below, which can be machine generated from the Octave function
+  ldpc_fsk_lib.m:ldpc_decode()
+
+  The include file also contains test input/output vectors for the LDPC
+  decoder for testing this program.
+
+  Build:
+
+    $ gcc -o ldpc_dec ldpc_dec.c mpdecode_core.c -Wall -lm -g
+
+  TODO:
+  [ ] C cmd line encoder
+      [ ] SD output option
+      [ ] Es/No option for testing
+  [ ] decoder
+      [X] test mode or file I/O (incl stdin/stdout)
+      [X] Octave code to generate include file
+          + MAX_ITER as well
+      [ ] check into SVN
+      [ ] enc/dec running on cmd line
+      [ ] fsk_demod modified for soft decisions
+      [ ] drs232 modified for SD
+          + use UW syn cin this program to check BER with coding
+      [ ] revisit CML support, maybe blog post
+*/
+
+#include <assert.h>
+#include <errno.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include "mpdecode_core.h"
+
+/* Machine generated consts, H_rows, H_cols, test input/output data to
+   change LDPC code regenerate this file. */
+
+#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;
+    
+    /* derive some parameters */
+
+    max_iter   = MAX_ITER;
+    CodeLength = CODELENGTH;                    /* length of entire codeword */
+    NumberParityBits = NUMBERPARITYBITS;
+	
+    if (argc < 2) {
+        fprintf(stderr, "usage: %s --test\n", argv[0]);
+        fprintf(stderr, "usage: %s InOneSDSymbolPerDouble OutOneBitPerInt\n", argv[0]);
+        exit(0);
+    }
+
+    int *DecodedBits = calloc( max_iter*CodeLength, sizeof( int ) );
+    int *ParityCheckCount = calloc( max_iter, sizeof(int) );
+
+    if (!strcmp(argv[1],"--test")) {
+
+        /* test mode --------------------------------------------------------*/
+
+        fprintf(stderr, "Starting test using pre-compiled test data .....\n");
+        fprintf(stderr, "Codeword length: %d\n",  CodeLength);
+        fprintf(stderr, "Parity Bits....: %d\n",  NumberParityBits);
+
+        num_runs = 100; num_ok = 0;
+
+        for(r=0; r<num_runs; r++) {
+
+            run_ldpc_decoder(DecodedBits, ParityCheckCount, input);
+
+            /* Check output. Bill has modified decoded to write number of parity checks that are correct
+               to the BitErrors array.  So if decoding is correct the final BitErrors entry will have 
+               NumberParityBits */
+
+            /* find output bit row where decoding has finished */
+
+            int ok = 0;
+            for (i=0;i<max_iter;i++) {
+                if (ParityCheckCount[i] == NumberParityBits) {
+                    for (j=0;j<CodeLength;j++) {
+                        if (output[i + j*max_iter] == DecodedBits[i+j*max_iter])                    
+                            ok++;
+                    }
+                }
+            }
+            if (ok == CodeLength)
+                num_ok++;            
+        }
+
+        fprintf(stderr, "test runs......: %d\n",  num_runs);
+        fprintf(stderr, "test runs OK...: %d\n",  num_ok);
+        if (num_runs == num_ok)
+            fprintf(stderr, "test runs OK...: PASS\n");
+        else
+            fprintf(stderr, "test runs OK...: FAIL\n");
+    }
+    else {
+        FILE *fin, *fout;
+
+        /* File I/O mode ------------------------------------------------*/
+
+        if (strcmp(argv[1], "-")  == 0) fin = stdin;
+        else if ( (fin = fopen(argv[1],"rb")) == NULL ) {
+            fprintf(stderr, "Error opening input SD 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 bit file: %s: %s.\n",
+                    argv[3], strerror(errno));
+            exit(1);
+        }
+
+        double *input_double  =  calloc(CodeLength, sizeof(double));
+
+        while(fread(input_double, sizeof(double), CodeLength, fin) == 1) {
+            run_ldpc_decoder(DecodedBits, ParityCheckCount, input_double);
+            fwrite(DecodedBits, sizeof(int), CodeLength, fout);
+        }
+
+        free(input_double);
+    }
+
+    /* Clean up memory */
+
+    free(ParityCheckCount);
+    free(DecodedBits);
+
+    return 0;
+}
+
+
+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;
+	
+#ifdef TT
+    /* check H_rows[] */
+
+    int length;
+    int gd = 0, bd = 0;
+    length = NumberRowsHcols*max_col_weight;
+    for (i=0;i<length;i++) {
+        if (H_cols[i] == H_cols2[i])
+            gd++;
+        else {
+            bd++;
+            printf("%d H_cols: %f H_cols2: %f\n", i, H_cols[i], H_cols2[i]);
+        }
+        if (bd == 10)
+            exit(0);
+    }
+    printf("gd: %d bd: %d\n", gd, bd);
+#endif
+
+    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) );
+	
+    /* 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
new file mode 100644
index 00000000..b308fd91
--- /dev/null
+++ b/codec2-dev/src/mpdecode_core.c
@@ -0,0 +1,561 @@
+/*
+  FILE...: mpdecode_core.c
+  AUTHOR.: Matthew C. Valenti, Rohit Iyer Seshadri, David Rowe
+  CREATED: Sep 2016
+
+  C-callable core functions moved from MpDecode.c, so they can be used for
+  Octave and C programs.
+*/
+
+#include <math.h>
+#include <stdlib.h>
+#include "mpdecode_core.h"
+
+/* Phi function */
+static float phi0(
+		  float x )
+{
+  float z;
+
+  if (x>10)
+    return( 0 );
+  else if (x< 9.08e-5 )
+    return( 10 );
+  else if (x > 9)
+    return( 1.6881e-4 );
+  /* return( 1.4970e-004 ); */
+  else if (x > 8)
+    return( 4.5887e-4 );
+  /* return( 4.0694e-004 ); */
+  else if (x > 7)
+    return( 1.2473e-3 );
+  /* return( 1.1062e-003 ); */
+  else if (x > 6)
+    return( 3.3906e-3 );
+  /* return( 3.0069e-003 ); */
+  else if (x > 5)
+    return( 9.2168e-3 );
+  /* return( 8.1736e-003 ); */
+  else {
+    z = (float) exp(x);
+    return( (float) log( (z+1)/(z-1) ) ); 
+  }
+}
+
+static float correction(
+			float xinput )
+{
+  if (xinput > 2.625 )
+    return( 0 );
+  else if (xinput < 1 )
+    return( -0.375*xinput + 0.6825 );
+  else 
+    return( -0.1875*xinput + 0.5 );
+
+}
+
+static float LambdaAPPstar(	float mag1,
+				float mag2 )
+{
+  if (mag1 > mag2)
+    return( fabs( mag2 + correction( mag1 + mag2 ) - correction( mag1 - mag2 ) ) );
+  else
+    return( fabs( mag1 + correction( mag1 + mag2 ) - correction( mag2 - mag1 ) ) );
+}
+
+void init_c_v_nodes(struct c_node *c_nodes, 
+                    int     shift, 
+                    int     NumberParityBits, 
+                    int     max_row_weight,
+                    double *H_rows,
+                    int     H1,
+                    int     CodeLength,
+                    struct v_node *v_nodes, 
+                    int     NumberRowsHcols, 
+                    double *H_cols,
+                    int     max_col_weight,
+                    int     dec_type,
+                    double *input)
+{
+    int i, j, k, count, cnt, c_index, v_index;
+
+    /* first determine the degree of each c-node */
+	
+    if (shift ==0){	
+        for (i=0;i<NumberParityBits;i++) {
+            count = 0;
+            for (j=0;j<max_row_weight;j++) {
+                if ( H_rows[i+j*NumberParityBits] > 0 ) {
+                    count++;
+                }
+            }
+            c_nodes[i].degree = count;
+            if (H1){
+                if (i==0){
+                    c_nodes[i].degree=count+1;
+                } 
+                else{
+                    c_nodes[i].degree=count+2;
+                }
+            }
+        }
+    }	
+    else{	
+        cnt=0; 
+        for (i=0;i<(NumberParityBits/shift);i++) {		
+            for (k=0;k<shift;k++){
+                count = 0;
+                for (j=0;j<max_row_weight;j++) {
+                    if ( H_rows[cnt+j*NumberParityBits] > 0 ) {
+                        count++;
+                    }
+                }
+                c_nodes[cnt].degree = count;
+                if ((i==0)||(i==(NumberParityBits/shift)-1)){
+                    c_nodes[cnt].degree=count+1;
+                } 
+                else{
+                    c_nodes[cnt].degree=count+2;
+                }
+                cnt++;
+            }	   
+        }	
+    }
+			
+    if (H1){
+
+        if (shift ==0){
+            for (i=0;i<NumberParityBits;i++) {
+                /* now that we know the size, we can dynamically allocate memory */
+                c_nodes[i].index =  calloc( c_nodes[i].degree, sizeof( int ) );
+                c_nodes[i].message =calloc( c_nodes[i].degree, sizeof( float ) );
+                c_nodes[i].socket = calloc( c_nodes[i].degree, sizeof( int ) );
+			
+                for (j=0;j<c_nodes[i].degree-2;j++) {
+                    c_nodes[i].index[j] = (int) (H_rows[i+j*NumberParityBits] - 1);
+                }			    
+                j=c_nodes[i].degree-2;
+		
+                if (i==0){
+                    c_nodes[i].index[j] = (int) (H_rows[i+j*NumberParityBits] - 1);			        
+                }
+                else {
+                    c_nodes[i].index[j] = (CodeLength-NumberParityBits)+i-1;
+                }
+						
+                j=c_nodes[i].degree-1;			
+                c_nodes[i].index[j] = (CodeLength-NumberParityBits)+i;
+		    
+            }		      
+        }		
+        if (shift >0){
+            cnt=0;
+            for (i=0;i<(NumberParityBits/shift);i++){		  
+		  
+                for (k =0;k<shift;k++){
+                    c_nodes[cnt].index =  calloc( c_nodes[cnt].degree, sizeof( int ) );
+                    c_nodes[cnt].message =calloc( c_nodes[cnt].degree, sizeof( float ) );
+                    c_nodes[cnt].socket = calloc( c_nodes[cnt].degree, sizeof( int ) );
+			 		   
+                    for (j=0;j<c_nodes[cnt].degree-2;j++) {
+                        c_nodes[cnt].index[j] = (int) (H_rows[cnt+j*NumberParityBits] - 1);
+                    }			    
+                    j=c_nodes[cnt].degree-2;
+                    if ((i ==0)||(i==(NumberParityBits/shift-1))){
+                        c_nodes[cnt].index[j] = (int) (H_rows[cnt+j*NumberParityBits] - 1);	
+                    }
+                    else{
+                        c_nodes[cnt].index[j] = (CodeLength-NumberParityBits)+k+shift*(i);
+                    }			
+                    j=c_nodes[cnt].degree-1;    		     
+                    c_nodes[cnt].index[j] = (CodeLength-NumberParityBits)+k+shift*(i+1);
+                    if (i== (NumberParityBits/shift-1))
+                        {
+                            c_nodes[cnt].index[j] = (CodeLength-NumberParityBits)+k+shift*(i);
+                        }
+                    cnt++;			    
+                } 
+            }
+        }
+				
+    } else {
+        for (i=0;i<NumberParityBits;i++) {
+            /* now that we know the size, we can dynamically allocate memory */
+            c_nodes[i].index =  calloc( c_nodes[i].degree, sizeof( int ) );
+            c_nodes[i].message =calloc( c_nodes[i].degree, sizeof( float ) );
+            c_nodes[i].socket = calloc( c_nodes[i].degree, sizeof( int ) );
+            for (j=0;j<c_nodes[i].degree;j++){
+                c_nodes[i].index[j] = (int) (H_rows[i+j*NumberParityBits] - 1);
+            }			
+        }
+    }	
+
+
+    /* determine degree of each v-node */
+
+    for(i=0;i<(CodeLength-NumberParityBits+shift);i++){
+        count=0;		
+        for (j=0;j<max_col_weight;j++) {
+            if ( H_cols[i+j*NumberRowsHcols] > 0 ) {
+                count++;
+            }
+        }
+        v_nodes[i].degree = count;
+    }
+	
+    for(i=CodeLength-NumberParityBits+shift;i<CodeLength;i++){
+        count=0;
+        if (H1){
+            if(i!=CodeLength-1){
+                v_nodes[i].degree=2;
+            }  else{
+                v_nodes[i].degree=1;
+            }	 
+			
+        } else{
+            for (j=0;j<max_col_weight;j++) {
+                if ( H_cols[i+j*NumberRowsHcols] > 0 ) {
+                    count++;
+                }
+            }      
+            v_nodes[i].degree = count;	 	  
+        }	 
+    }  
+	 
+    if (shift>0){
+        v_nodes[CodeLength-1].degree =v_nodes[CodeLength-1].degree+1;	                   
+    }
+    	
+    
+    /* set up v_nodes */
+
+    for (i=0;i<CodeLength;i++) {
+        /* allocate memory according to the degree of the v-node */
+        v_nodes[i].index = calloc( v_nodes[i].degree, sizeof( int ) );
+        v_nodes[i].message = calloc( v_nodes[i].degree, sizeof( float ) );
+        v_nodes[i].sign = calloc( v_nodes[i].degree, sizeof( int ) );
+        v_nodes[i].socket = calloc( v_nodes[i].degree, sizeof( int ) );
+		
+        /* index tells which c-nodes this v-node is connected to */
+        v_nodes[i].initial_value = input[i];
+        count=0;
+
+        for (j=0;j<v_nodes[i].degree;j++) {			
+            if ((H1)&& (i>=CodeLength-NumberParityBits+shift)){
+                v_nodes[i].index[j]=i-(CodeLength-NumberParityBits+shift)+count;
+                if (shift ==0){
+                    count=count+1;
+                }
+                else{
+                    count=count+shift;
+                }
+            } else  {
+                v_nodes[i].index[j] = (int) (H_cols[i+j*NumberRowsHcols] - 1);
+            }			
+						
+            /* search the connected c-node for the proper message value */
+            for (c_index=0;c_index<c_nodes[ v_nodes[i].index[j] ].degree;c_index++)
+                if ( c_nodes[ v_nodes[i].index[j] ].index[c_index] == i ) {
+                    v_nodes[i].socket[j] = c_index;
+                    break;
+                }				
+            /* initialize v-node with received LLR */			
+            if ( dec_type == 1)
+                v_nodes[i].message[j] = fabs(input[i]);
+            else
+                v_nodes[i].message[j] = phi0( fabs(input[i]) );
+				
+            if (input[i] < 0)
+                v_nodes[i].sign[j] = 1;			
+        }
+	
+    }
+	
+	
+	
+    /* now finish setting up the c_nodes */
+    for (i=0;i<NumberParityBits;i++) {		
+        /* index tells which v-nodes this c-node is connected to */
+        for (j=0;j<c_nodes[i].degree;j++) {			
+            /* search the connected v-node for the proper message value */
+            for (v_index=0;v_index<v_nodes[ c_nodes[i].index[j] ].degree;v_index++)
+                if (v_nodes[ c_nodes[i].index[j] ].index[v_index] == i ) {
+                    c_nodes[i].socket[j] = v_index;
+                    break;
+                }
+        }
+    }
+
+}
+
+
+/* function for doing the MP decoding */
+void ApproximateMinStar(	 int	  BitErrors[],
+				 int      DecodedBits[],
+				 struct c_node c_nodes[],
+				 struct v_node v_nodes[],
+				 int	  CodeLength,
+				 int	  NumberParityBits,
+				 int	  max_iter )
+{
+  int i,j, iter;
+  int sign;
+  float temp_sum;
+  float Qi;
+
+  float delta, minval, deltaAPP;
+  int mink;
+
+  for (iter=0;iter<max_iter;iter++) {
+    /* update r */
+    for (j=0;j<NumberParityBits;j++) {	
+      /* start new code for approximate-min-star */
+      mink = 0;
+      sign = v_nodes[ c_nodes[j].index[0] ].sign[ c_nodes[j].socket[0] ];
+      minval = v_nodes[ c_nodes[j].index[0] ].message[ c_nodes[j].socket[0] ];
+		
+      for (i=1;i<c_nodes[j].degree;i++) {
+	/* first find the minimum magnitude input message */
+	if ( v_nodes[ c_nodes[j].index[i] ].message[ c_nodes[j].socket[i] ] < minval ) {
+	  mink = i;
+	  minval = v_nodes[ c_nodes[j].index[i] ].message[ c_nodes[j].socket[i] ];							
+	}
+	/* update the aggregate sign */
+	sign ^= v_nodes[ c_nodes[j].index[i] ].sign[ c_nodes[j].socket[i] ];
+      }
+
+      /* find the magnitude to send out the minimum input magnitude branch */
+      if ( mink == 0 ) {
+	delta = v_nodes[ c_nodes[j].index[1] ].message[ c_nodes[j].socket[1] ];
+	for (i=2;i<c_nodes[j].degree;i++) {
+	  delta = LambdaAPPstar( delta, v_nodes[ c_nodes[j].index[i] ].message[ c_nodes[j].socket[i] ] );
+	}
+      } else {
+	delta = v_nodes[ c_nodes[j].index[0] ].message[ c_nodes[j].socket[0] ];
+	for (i=1;i<c_nodes[j].degree;i++) {
+	  if ( i != mink )
+	    delta = LambdaAPPstar( delta, v_nodes[ c_nodes[j].index[i] ].message[ c_nodes[j].socket[i] ] );
+	}
+      }
+
+      deltaAPP = LambdaAPPstar( delta, v_nodes[ c_nodes[j].index[mink] ].message[ c_nodes[j].socket[mink] ] );
+
+      /* compute outgoing messages */
+      for (i=0;i<c_nodes[j].degree;i++) {
+	if ( i == mink ) {
+	  if ( sign^v_nodes[ c_nodes[j].index[i] ].sign[ c_nodes[j].socket[i] ] )
+	    c_nodes[j].message[i] = - delta;
+	  else
+	    c_nodes[j].message[i] = delta;
+	} else {
+	  if ( sign^v_nodes[ c_nodes[j].index[i] ].sign[ c_nodes[j].socket[i] ] )
+	    c_nodes[j].message[i] = - deltaAPP;
+	  else
+	    c_nodes[j].message[i] = deltaAPP;
+	}
+      }		
+    }
+
+    /* update q */
+    for (i=0;i<CodeLength;i++) {
+
+      /* first compute the LLR */
+      Qi = v_nodes[i].initial_value;
+      for (j=0;j<v_nodes[i].degree;j++) {				
+	Qi += c_nodes[ v_nodes[i].index[j] ].message[ v_nodes[i].socket[j] ];
+      }
+
+      /* make hard decision */			
+      if (Qi < 0) {
+	DecodedBits[iter+max_iter*i] = 1;
+	BitErrors[iter]++;
+      }
+
+      /* now subtract to get the extrinsic information */
+      for (j=0;j<v_nodes[i].degree;j++) {
+	temp_sum = Qi - c_nodes[ v_nodes[i].index[j] ].message[ v_nodes[i].socket[j] ];
+				
+	v_nodes[i].message[j] = fabs( temp_sum );
+	if (temp_sum > 0)
+	  v_nodes[i].sign[j] = 0;
+	else
+	  v_nodes[i].sign[j] = 1;
+      }
+    }
+
+    /* detect errors */
+    if (BitErrors[iter] == 0)
+      break; 
+  }
+}
+
+
+/* function for doing the MP decoding */
+void MinSum(		 int	  BitErrors[],
+				 int      DecodedBits[],
+				 struct c_node c_nodes[],
+				 struct v_node v_nodes[],
+				 int	  CodeLength,
+				 int	  NumberParityBits,
+				 int	  max_iter, 
+				 float    r_scale_factor,
+				 float    q_scale_factor, 
+				 int      data[] )
+{
+  int i,j, iter, i_prime;
+  float min_beta;
+  int sign;
+  float temp_sum;
+  float Qi;
+
+  for (iter=0;iter<max_iter;iter++) {
+
+    /* update r */
+    for (j=0;j<NumberParityBits;j++) {
+      sign = 0;
+      for (i=0;i<c_nodes[j].degree;i++) 
+	sign ^= v_nodes[ c_nodes[j].index[i] ].sign[ c_nodes[j].socket[i] ];
+
+      for (i=0;i<c_nodes[j].degree;i++) {
+	min_beta = 1000;		
+								
+	for (i_prime=0;i_prime<c_nodes[j].degree;i_prime++) 
+	  if ( ( v_nodes[ c_nodes[j].index[i_prime] ].message[c_nodes[j].socket[i_prime]] < min_beta )&&(i_prime != i) )
+	    min_beta = v_nodes[ c_nodes[j].index[i_prime] ].message[c_nodes[j].socket[i_prime]];
+
+	if ( sign^v_nodes[ c_nodes[j].index[i] ].sign[ c_nodes[j].socket[i] ] )
+	  c_nodes[j].message[i] = -min_beta*r_scale_factor;
+	else
+	  c_nodes[j].message[i] = min_beta*r_scale_factor;
+      }
+    }
+
+    /* update q */
+    for (i=0;i<CodeLength;i++) {
+			
+      /* first compute the LLR */
+      Qi = v_nodes[i].initial_value;
+      for (j=0;j<v_nodes[i].degree;j++) {				
+	Qi += c_nodes[ v_nodes[i].index[j] ].message[ v_nodes[i].socket[j] ];
+      }
+
+      /* make hard decision */			
+      if (Qi < 0) {
+	DecodedBits[iter+max_iter*i] = 1;
+      }
+
+      /* now subtract to get the extrinsic information */
+      for (j=0;j<v_nodes[i].degree;j++) {
+	temp_sum = Qi - c_nodes[ v_nodes[i].index[j] ].message[ v_nodes[i].socket[j] ];
+				
+	v_nodes[i].message[j] = fabs( temp_sum )*q_scale_factor;
+	if (temp_sum > 0)
+	  v_nodes[i].sign[j] = 0;
+	else
+	  v_nodes[i].sign[j] = 1;
+      }
+    }
+
+    /* count data bit errors, assuming that it is systematic */
+    for (i=0;i<CodeLength-NumberParityBits;i++)
+      if ( DecodedBits[iter+max_iter*i] != data[i] )
+	BitErrors[iter]++;
+
+    /* detect errors */
+    if (BitErrors[iter] == 0)
+      break; 
+  }
+}
+
+
+/* function for doing the MP decoding */
+void SumProduct(	 int	  BitErrors[],
+			 int      DecodedBits[],
+			 struct c_node c_nodes[],
+			 struct v_node v_nodes[],
+			 int	  CodeLength,
+			 int	  NumberParityBits,
+			 int	  max_iter,
+			 float    r_scale_factor,
+			 float    q_scale_factor, 
+			 int      data[] )
+{
+  int i,j, iter;
+  float phi_sum;
+  int sign;
+  float temp_sum;
+  float Qi;
+  int   ssum; 
+
+  for (iter=0;iter<max_iter;iter++) {
+    /* update r */
+    ssum = 0; 
+    for (j=0;j<NumberParityBits;j++) {		
+      sign = v_nodes[ c_nodes[j].index[0] ].sign[ c_nodes[j].socket[0] ];
+      phi_sum = v_nodes[ c_nodes[j].index[0] ].message[ c_nodes[j].socket[0] ];
+			
+      for (i=1;i<c_nodes[j].degree;i++) {
+	phi_sum += v_nodes[ c_nodes[j].index[i] ].message[ c_nodes[j].socket[i] ];
+	sign ^= v_nodes[ c_nodes[j].index[i] ].sign[ c_nodes[j].socket[i] ];
+      }
+
+      if (sign==0) ssum++; 
+		
+      for (i=0;i<c_nodes[j].degree;i++) {
+	if ( sign^v_nodes[ c_nodes[j].index[i] ].sign[ c_nodes[j].socket[i] ] ) {
+	  c_nodes[j].message[i] = -phi0( phi_sum - v_nodes[ c_nodes[j].index[i] ].message[ c_nodes[j].socket[i] ] )*r_scale_factor;
+	} else
+	  c_nodes[j].message[i] = phi0( phi_sum - v_nodes[ c_nodes[j].index[i] ].message[ c_nodes[j].socket[i] ] )*r_scale_factor;
+      }
+    }
+
+    /* update q */
+    for (i=0;i<CodeLength;i++) {
+
+      /* first compute the LLR */
+      Qi = v_nodes[i].initial_value;
+      for (j=0;j<v_nodes[i].degree;j++) {				
+	Qi += c_nodes[ v_nodes[i].index[j] ].message[ v_nodes[i].socket[j] ];
+      }
+
+      /* make hard decision */			
+      if (Qi < 0) {
+	DecodedBits[iter+max_iter*i] = 1;
+      }
+
+      /* now subtract to get the extrinsic information */
+      for (j=0;j<v_nodes[i].degree;j++) {
+	temp_sum = Qi - c_nodes[ v_nodes[i].index[j] ].message[ v_nodes[i].socket[j] ];
+				
+	v_nodes[i].message[j] = phi0( fabs( temp_sum ) )*q_scale_factor;
+	if (temp_sum > 0)
+	  v_nodes[i].sign[j] = 0;
+	else
+	  v_nodes[i].sign[j] = 1;
+      }
+    }
+
+    /* count data bit errors, assuming that it is systematic */
+    for (i=0;i<CodeLength-NumberParityBits;i++)
+      if ( DecodedBits[iter+max_iter*i] != data[i] )
+	BitErrors[iter]++;
+
+    /* Halt if zero errors */
+    if (BitErrors[iter] == 0)
+      break; 
+
+    // added by Bill -- reuse the BitErrors array to count PCs
+    // count the number of PC satisfied and exit if all OK
+    BitErrors[iter] = ssum;
+    if (ssum==NumberParityBits) break;
+
+
+  }
+   
+  // printf(" ssum is %d \n",   ssum); 
+   
+ 
+    
+}
+
+
diff --git a/codec2-dev/src/mpdecode_core.h b/codec2-dev/src/mpdecode_core.h
new file mode 100644
index 00000000..be430230
--- /dev/null
+++ b/codec2-dev/src/mpdecode_core.h
@@ -0,0 +1,74 @@
+/*
+  FILE...: mpdecode_core.h
+  AUTHOR.: David Rowe
+  CREATED: Sep 2016
+
+  C-callable core functions for MpDecode, so they can be used for
+  Octave and C programs.
+*/
+
+#ifndef __MPDECODE_CORE__
+#define __MPDECODE_CORE__
+
+struct v_node {
+  int degree;
+  float initial_value;
+  int *index;  /* the index of a c_node it is connected to */
+  int *socket; /* socket number at the c_node */
+  float *message;     
+  int *sign;
+};
+
+struct c_node {
+  int degree;
+  int *index;                     
+  float *message;     
+  int *socket; /* socket number at the v_node */
+};
+
+void init_c_v_nodes(struct c_node *c_nodes, 
+                    int     shift, 
+                    int     NumberParityBits, 
+                    int     max_row_weight,
+                    double *H_rows,
+                    int     H1,
+                    int     CodeLength,
+                    struct v_node *v_nodes, 
+                    int     NumberRowsHcols, 
+                    double *H_cols,
+                    int     max_col_weight,
+                    int     dec_type,
+                    double *input);
+
+void ApproximateMinStar(	 int	  BitErrors[],
+				 int      DecodedBits[],
+				 struct c_node c_nodes[],
+				 struct v_node v_nodes[],
+				 int	  CodeLength,
+				 int	  NumberParityBits,
+				 int	  max_iter );
+
+void MinSum(		 int	  BitErrors[],
+				 int      DecodedBits[],
+				 struct c_node c_nodes[],
+				 struct v_node v_nodes[],
+				 int	  CodeLength,
+				 int	  NumberParityBits,
+				 int	  max_iter, 
+				 float    r_scale_factor,
+				 float    q_scale_factor, 
+				 int      data[] );
+
+
+void SumProduct(	 int	  BitErrors[],
+			 int      DecodedBits[],
+			 struct c_node c_nodes[],
+			 struct v_node v_nodes[],
+			 int	  CodeLength,
+			 int	  NumberParityBits,
+			 int	  max_iter,
+			 float    r_scale_factor,
+			 float    q_scale_factor, 
+			 int      data[]);
+
+#endif
-- 
2.25.1