mpdecode.c

iterative 解码库,matlab 和 c写的

		/* Halt if zero errors */
		if (BitErrors[iter] == 0)
			break; 

	}
}

/* main function that interfaces with MATLAB */
void mexFunction(
				 int            nlhs,
				 mxArray       *plhs[],
				 int            nrhs,
				 const mxArray *prhs[] )
{	int		max_iter, dec_type;
	int		max_row_weight, max_col_weight;
	int		NumberParityBits, CodeLength;
	double  *H_rows, *H_cols;		/* Parity check matrix info */
	double  *input;		
	int     i, j, count, v_index, c_index;
	int		*DecodedBits;	/* Output of the decoder.  Is an array of size iter by CodeLength */
	int		*BitErrors;		/* Number of errors at each iteration */
	double  *errors_p, *output_p;
	struct c_node *c_nodes;
	struct v_node *v_nodes;
	float   q_scale_factor, r_scale_factor;
	double  *data;
	int     *data_int;
	int     DataLength;
	int     NumberRowsHcols;
	int		H1;
	int     shift;
	int     cnt;
	int     k;

	/* default values */
	max_iter  = 30;
	dec_type  = 0;
	q_scale_factor = 1;
	r_scale_factor = 1;

	/* Check for proper number of arguments */
	if ( (nrhs < 3 )|(nlhs  > 2) ) {
		mexErrMsgTxt("Usage: [output, errors] = MpDecode(input, H_rows, H_cols, max_iter, dec_type, r_scale_factor, q_scale_factor, data )");
	} else {
		/* first input is the received data in LLR form */
		input = mxGetPr(INPUT);	

		/* second input is H_rows matrix */
		H_rows = mxGetPr( HROWS );

		/* third input is H_cols matrix */
		H_cols = mxGetPr( HCOLS );

		/* derive some parameters */
		CodeLength = mxGetN(INPUT); /* number of coded bits */
		NumberParityBits = mxGetM( HROWS );
		NumberRowsHcols=mxGetM( HCOLS );
		shift=(NumberParityBits+ NumberRowsHcols)-CodeLength;
		
		
		
		if (NumberRowsHcols ==CodeLength){
		      H1=0;
              shift=0;
		} else {
		       H1=1;
		}
	
		if (( CodeLength != NumberRowsHcols ) && (CodeLength-NumberParityBits + shift!= NumberRowsHcols))
		 	 mexErrMsgTxt("Error: Number of rows in H_cols must equal number of received bits or number of data bits");
		
		max_row_weight = mxGetN( HROWS );
		max_col_weight = mxGetN( HCOLS );
	} 
	
	/* initialize c-node structures */
    c_nodes = calloc( NumberParityBits, sizeof( struct c_node ) );
	
	/* 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);
			}			
		}
	}	

	/* initialize v-node structures */
	v_nodes = calloc( CodeLength, sizeof( struct v_node));
	
	/* 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;	                   
    	}
    	
    

	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;
				}
		}
	}
	
	
	if (nrhs > 3 ) {
		/* fourth input (optional) is maximum number of iterations */
		max_iter   = (int) *mxGetPr(MAXITER);
	} if (nrhs > 4 ) {
		/* fifth input (optional) is the decoder type */
		dec_type  = (int) *mxGetPr(DECTYPE);
	} if (nrhs > 5 ) {
		/* next input is the factor for extrinsic info scaling */
		r_scale_factor = (float) *mxGetPr(RSCALEFACTOR);
	}  if (nrhs > 6 ) {
		/* next input is the factor for extrinsic info scaling */
		q_scale_factor = (float) *mxGetPr(QSCALEFACTOR);
	}  
	
	DataLength = CodeLength - NumberParityBits;
	data_int = calloc( DataLength, sizeof(int) );
	
	if (nrhs > 7 ) {
		/* next input is the data */
		data = mxGetPr(DATA);	
		if ( DataLength != mxGetN(DATA) ) /* number of data bits */
			mexErrMsgTxt("Error: Incorrect number of data bits");
		
		/* cast the input into a vector of integers */
		for (i=0;i<DataLength;i++) {
			data_int[i] = (int) data[i];
		}
		
	}

	/* create matrices for the decoded data */		
	OUTPUT = mxCreateDoubleMatrix(max_iter, CodeLength, mxREAL );
	output_p = mxGetPr(OUTPUT);	

	/* Decode */
	DecodedBits = calloc( max_iter*CodeLength, sizeof( int ) );
	BitErrors = calloc( max_iter, sizeof(int) );

    /* Call function to do the actual decoding */
	if ( dec_type == 1) {
		MinSum( BitErrors, DecodedBits, c_nodes, v_nodes, CodeLength, 
		   NumberParityBits, max_iter, r_scale_factor, q_scale_factor, data_int );
	} else if ( dec_type == 2) {
		mexErrMsgTxt("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( BitErrors, DecodedBits, c_nodes, v_nodes, CodeLength, 
		   NumberParityBits, max_iter, r_scale_factor, q_scale_factor, data_int ); 
	}

	/* cast to output */
    for (i=0;i<max_iter;i++) {
		for (j=0;j<CodeLength;j++) {
			output_p[i + j*max_iter] = DecodedBits[i+j*max_iter];
		}
	}
	
	if (nlhs > 1 ) {
		/* second output is a count of the number of errors */
		ERRORS = mxCreateDoubleMatrix(max_iter, 1, mxREAL);
		errors_p = mxGetPr(ERRORS);
		
		/* cast to output */    
		for (i=0;i<max_iter;i++) {
			errors_p[i] = BitErrors[i];
		}
	}
	
	/* Clean up memory */
	free( BitErrors );
	free( DecodedBits );
	free( data_int );
	
	/* printf( "Cleaning c-node elements\n" ); */
	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 );
	
	return;
}