--- /dev/null
+/* \r
+ FILE...: ldpc_dec.c\r
+ AUTHOR.: Matthew C. Valenti, Rohit Iyer Seshadri, David Rowe\r
+ CREATED: Sep 2016\r
+\r
+ Command line C LDPC decoder derived from MpDecode.c in the CML\r
+ library. Allows us to run the same decoder in Octave and C. The\r
+ code is defined by the parameters and array stored in the include\r
+ file below, which can be machine generated from the Octave function\r
+ ldpc_fsk_lib.m:ldpc_decode()\r
+\r
+ The include file also contains test input/output vectors for the LDPC\r
+ decoder for testing this program.\r
+\r
+ Build:\r
+\r
+ $ gcc -o ldpc_dec ldpc_dec.c mpdecode_core.c -Wall -lm -g\r
+\r
+ TODO:\r
+ [ ] C cmd line encoder\r
+ [ ] SD output option\r
+ [ ] Es/No option for testing\r
+ [ ] decoder\r
+ [X] test mode or file I/O (incl stdin/stdout)\r
+ [X] Octave code to generate include file\r
+ + MAX_ITER as well\r
+ [ ] check into SVN\r
+ [ ] enc/dec running on cmd line\r
+ [ ] fsk_demod modified for soft decisions\r
+ [ ] drs232 modified for SD\r
+ + use UW syn cin this program to check BER with coding\r
+ [ ] revisit CML support, maybe blog post\r
+*/\r
+\r
+#include <assert.h>\r
+#include <errno.h>\r
+#include <math.h>\r
+#include <stdlib.h>\r
+#include <string.h>\r
+#include <stdio.h>\r
+#include "mpdecode_core.h"\r
+\r
+/* Machine generated consts, H_rows, H_cols, test input/output data to\r
+ change LDPC code regenerate this file. */\r
+\r
+#include "ldpc_code.h" \r
+\r
+void run_ldpc_decoder(int DecodedBits[], int ParityCheckCount[], double input[]);\r
+\r
+int main(int argc, char *argv[])\r
+{ \r
+ int CodeLength, NumberParityBits, max_iter;\r
+ int i, j, r, num_ok, num_runs;\r
+ \r
+ /* derive some parameters */\r
+\r
+ max_iter = MAX_ITER;\r
+ CodeLength = CODELENGTH; /* length of entire codeword */\r
+ NumberParityBits = NUMBERPARITYBITS;\r
+ \r
+ if (argc < 2) {\r
+ fprintf(stderr, "usage: %s --test\n", argv[0]);\r
+ fprintf(stderr, "usage: %s InOneSDSymbolPerDouble OutOneBitPerInt\n", argv[0]);\r
+ exit(0);\r
+ }\r
+\r
+ int *DecodedBits = calloc( max_iter*CodeLength, sizeof( int ) );\r
+ int *ParityCheckCount = calloc( max_iter, sizeof(int) );\r
+\r
+ if (!strcmp(argv[1],"--test")) {\r
+\r
+ /* test mode --------------------------------------------------------*/\r
+\r
+ fprintf(stderr, "Starting test using pre-compiled test data .....\n");\r
+ fprintf(stderr, "Codeword length: %d\n", CodeLength);\r
+ fprintf(stderr, "Parity Bits....: %d\n", NumberParityBits);\r
+\r
+ num_runs = 100; num_ok = 0;\r
+\r
+ for(r=0; r<num_runs; r++) {\r
+\r
+ run_ldpc_decoder(DecodedBits, ParityCheckCount, input);\r
+\r
+ /* Check output. Bill has modified decoded to write number of parity checks that are correct\r
+ to the BitErrors array. So if decoding is correct the final BitErrors entry will have \r
+ NumberParityBits */\r
+\r
+ /* find output bit row where decoding has finished */\r
+\r
+ int ok = 0;\r
+ for (i=0;i<max_iter;i++) {\r
+ if (ParityCheckCount[i] == NumberParityBits) {\r
+ for (j=0;j<CodeLength;j++) {\r
+ if (output[i + j*max_iter] == DecodedBits[i+j*max_iter]) \r
+ ok++;\r
+ }\r
+ }\r
+ }\r
+ if (ok == CodeLength)\r
+ num_ok++; \r
+ }\r
+\r
+ fprintf(stderr, "test runs......: %d\n", num_runs);\r
+ fprintf(stderr, "test runs OK...: %d\n", num_ok);\r
+ if (num_runs == num_ok)\r
+ fprintf(stderr, "test runs OK...: PASS\n");\r
+ else\r
+ fprintf(stderr, "test runs OK...: FAIL\n");\r
+ }\r
+ else {\r
+ FILE *fin, *fout;\r
+\r
+ /* File I/O mode ------------------------------------------------*/\r
+\r
+ if (strcmp(argv[1], "-") == 0) fin = stdin;\r
+ else if ( (fin = fopen(argv[1],"rb")) == NULL ) {\r
+ fprintf(stderr, "Error opening input SD file: %s: %s.\n",\r
+ argv[2], strerror(errno));\r
+ exit(1);\r
+ }\r
+ \r
+ if (strcmp(argv[2], "-") == 0) fout = stdout;\r
+ else if ( (fout = fopen(argv[2],"wb")) == NULL ) {\r
+ fprintf(stderr, "Error opening output bit file: %s: %s.\n",\r
+ argv[3], strerror(errno));\r
+ exit(1);\r
+ }\r
+\r
+ double *input_double = calloc(CodeLength, sizeof(double));\r
+\r
+ while(fread(input_double, sizeof(double), CodeLength, fin) == 1) {\r
+ run_ldpc_decoder(DecodedBits, ParityCheckCount, input_double);\r
+ fwrite(DecodedBits, sizeof(int), CodeLength, fout);\r
+ }\r
+\r
+ free(input_double);\r
+ }\r
+\r
+ /* Clean up memory */\r
+\r
+ free(ParityCheckCount);\r
+ free(DecodedBits);\r
+\r
+ return 0;\r
+}\r
+\r
+\r
+void run_ldpc_decoder(int DecodedBits[], int ParityCheckCount[], double input[]) {\r
+ int max_iter, dec_type;\r
+ float q_scale_factor, r_scale_factor;\r
+ int max_row_weight, max_col_weight;\r
+ int CodeLength, NumberParityBits, NumberRowsHcols, shift, H1;\r
+ int i;\r
+ struct c_node *c_nodes;\r
+ struct v_node *v_nodes;\r
+ \r
+ /* default values */\r
+\r
+ max_iter = MAX_ITER;\r
+ dec_type = 0;\r
+ q_scale_factor = 1;\r
+ r_scale_factor = 1;\r
+\r
+ /* derive some parameters */\r
+\r
+ CodeLength = CODELENGTH; /* length of entire codeword */\r
+ NumberParityBits = NUMBERPARITYBITS;\r
+ NumberRowsHcols = NUMBERROWSHCOLS;\r
+\r
+ shift = (NumberParityBits + NumberRowsHcols) - CodeLength;\r
+ if (NumberRowsHcols == CodeLength) {\r
+ H1=0;\r
+ shift=0;\r
+ } else {\r
+ H1=1;\r
+ }\r
+ \r
+ max_row_weight = MAX_ROW_WEIGHT;\r
+ max_col_weight = MAX_COL_WEIGHT;\r
+ \r
+#ifdef TT\r
+ /* check H_rows[] */\r
+\r
+ int length;\r
+ int gd = 0, bd = 0;\r
+ length = NumberRowsHcols*max_col_weight;\r
+ for (i=0;i<length;i++) {\r
+ if (H_cols[i] == H_cols2[i])\r
+ gd++;\r
+ else {\r
+ bd++;\r
+ printf("%d H_cols: %f H_cols2: %f\n", i, H_cols[i], H_cols2[i]);\r
+ }\r
+ if (bd == 10)\r
+ exit(0);\r
+ }\r
+ printf("gd: %d bd: %d\n", gd, bd);\r
+#endif\r
+\r
+ c_nodes = calloc( NumberParityBits, sizeof( struct c_node ) );\r
+ v_nodes = calloc( CodeLength, sizeof( struct v_node));\r
+ \r
+ /* initialize c-node and v-node structures */\r
+\r
+ c_nodes = calloc( NumberParityBits, sizeof( struct c_node ) );\r
+ v_nodes = calloc( CodeLength, sizeof( struct v_node));\r
+ \r
+ init_c_v_nodes(c_nodes, shift, NumberParityBits, max_row_weight, H_rows, H1, CodeLength, \r
+ v_nodes, NumberRowsHcols, H_cols, max_col_weight, dec_type, input);\r
+\r
+ int DataLength = CodeLength - NumberParityBits;\r
+ int *data_int = calloc( DataLength, sizeof(int) );\r
+ \r
+ /* Call function to do the actual decoding */\r
+\r
+ if ( dec_type == 1) {\r
+ MinSum( ParityCheckCount, DecodedBits, c_nodes, v_nodes, CodeLength, \r
+ NumberParityBits, max_iter, r_scale_factor, q_scale_factor, data_int );\r
+ } else if ( dec_type == 2) {\r
+ fprintf(stderr, "dec_type = 2 not currently supported");\r
+ /* ApproximateMinStar( BitErrors, DecodedBits, c_nodes, v_nodes, \r
+ CodeLength, NumberParityBits, max_iter, r_scale_factor, q_scale_factor );*/\r
+ } else {\r
+ SumProduct( ParityCheckCount, DecodedBits, c_nodes, v_nodes, CodeLength, \r
+ NumberParityBits, max_iter, r_scale_factor, q_scale_factor, data_int ); \r
+ }\r
+\r
+ /* Clean up memory */\r
+\r
+ free( data_int );\r
+\r
+ /* Cleaning c-node elements */\r
+\r
+ for (i=0;i<NumberParityBits;i++) {\r
+ free( c_nodes[i].index );\r
+ free( c_nodes[i].message );\r
+ free( c_nodes[i].socket );\r
+ }\r
+ \r
+ /* printf( "Cleaning c-nodes \n" ); */\r
+ free( c_nodes );\r
+ \r
+ /* printf( "Cleaning v-node elements\n" ); */\r
+ for (i=0;i<CodeLength;i++) {\r
+ free( v_nodes[i].index);\r
+ free( v_nodes[i].sign );\r
+ free( v_nodes[i].message );\r
+ free( v_nodes[i].socket );\r
+ }\r
+ \r
+ /* printf( "Cleaning v-nodes \n" ); */\r
+ free( v_nodes );\r
+}\r
--- /dev/null
+/*
+ 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);
+
+
+
+}
+
+
projects / freetel-svn-tracking.git / commitdiff