mpdecode.c

1、HSDPA； 2、LTE； 3、turbo code； 4、Mobile WiMAX； 5、LDPC

/* file: MpDecode.c

   Description: Decode a block code using the message passing algorithm. 

   The calling syntax is:

      [output, errors] = MpDecode(input, H_rows, H_cols, [max_iter], [dec_type], [r_scale_factor], [q_scale_factor], [data] )

      output = matrix of dimension maxiter by N that has the decoded code bits for each iteration 
	  errors = (optional) column vector showing the number of (data bit) errors after each iteration.

      input  = the decoder input in LLR form
	  H_cols = a N row matrix specifying the locations of the nonzero entries in each column of the H matrix.
	           The number or columns in the matrix is the max column weight.
	           OR
	           a (K + shift) row matrix specifying locations of the nonzero entries in each coulmn of an extended IRA type 
	           sparse H1 matrix
	           
	  H_rows = a N-K row matrix specifying the locations of the nonzero entries in each row of the H matrix.
	           The number or columns in the matrix is the max row weight, unless this is for an H1 matrix,
			   in which case the last n-k-shift columns of the H matrix are equal to a known H2 matrix.
	  max_iter = (optional) the maximum number of decoder iterations (default = 30).
	  dec_type = (optional) the decoder type:
			   = 0 Sum-product (default)
			   = 1 Min-sum 
			   = 2 Approximate-min-star
	  r_scale_factor = (optional) amount to scale extrinsic output of c-nodes in min-sum decoding (default = 1) 
      q_scale_factor = (optional) amount to scale extrinsic output of v-nodes in min-sum decoding (default = 1)
	  data = (optional) a vector containing the data bits (used for counting errors and for early halting) (default all zeros)
	
   Copyright (C) 2006-2007, Matthew C. Valenti and Rohit Iyer Seshadri

   Last updated on Aug. 8, 2007

   Function MpDecode is part of the Iterative Solutions 
   Coded Modulation Library. The Iterative Solutions Coded Modulation 
   Library is free software; you can redistribute it and/or modify it 
   under the terms of the GNU Lesser General Public License as published 
   by the Free Software Foundation; either version 2.1 of the License, 
   or (at your option) any later version.

   This library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.
  
   You should have received a copy of the GNU Lesser General Public
   License along with this library; if not, write to the Free Software
   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
*/

#include <math.h>
#include <mex.h>
#include <matrix.h>
#include <stdlib.h>

/* Input Arguments */
#define INPUT       prhs[0]
#define HROWS		prhs[1]
#define HCOLS		prhs[2]
#define MAXITER     prhs[3]
#define DECTYPE     prhs[4]
#define RSCALEFACTOR prhs[5]
#define QSCALEFACTOR prhs[6]
#define DATA prhs[7]

/* Output Arguments */
#define OUTPUT      plhs[0]
#define	ERRORS      plhs[1]

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 */
};

/* 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 ) ) );
}

/* function for doing the MP decoding */
static 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 */
static 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, j_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 */
static 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;

	for (iter=0;iter<max_iter;iter++) {
		/* update r */
		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] ];
			}
			
			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]++;