logmapmex.c

802.16eCTC的logmap算法仿真

/* MATLAB's C-MEX turbo code Log-MAP decoder
 *
 * Copyright Bagawan S. Nugroho, 2006 
 *********************************************/

#include "mex.h"
#include "matrix.h"
#include <math.h>
#include <stdlib.h>
/* Define an infinity number */
#define INF 1e10
/* Define a small number (i.e. the smallest IEEE 754 double-precision) */
#define SMALL 2.2e-308

/* ----------- The main function ----------- */
void logMap( double *LOut, double *x, double *LIn, int decIndx, double *next, 
              double *nextSt, double *prev, double *prevSt, int len, int nrSt)

/* x, receive symbols
 * LIn, priory information
 * decIndx, decoder index (1/2)
 * next, next trellis output (-1/1)
 * nextSt, next state
 * prev, previous trellis output (-1/1) 
 * prevSt, previous state
 * len, the length of input vector x
 * nrSt, number of state 
 * 
 * LOut, log-likelihood ratio output vector. */
              
{
   int lenData, i, j, k;
   double X, Y, sum, Gamma0, Gamma1, tmp, tmp0, tmp1;
   double Alpha[1000][64], Beta[1000][64];   
   double Gamma[64], tmpMax[2000];
   
   /* Length of the decoded information */
   lenData = (int) len/2;
    
   /* Initialize Alpha */
   *(*Alpha) = 0;
   for( j = 1; j < nrSt; j++) {
      *(*Alpha + j) = -INF;
   } /* for j */
   
   /* Initialize Beta */
   if( decIndx == 1) {
      *(*(Beta + lenData - 1)) = 0;
      for( i = 1; i < nrSt; i++) {
         *(*(Beta + lenData - 1) + i) = -INF;
      } /* for i */
   } /* if */
   else if ( decIndx == 2) {
      for( i = 0; i < nrSt; i++) {
         *(*(Beta + lenData - 1) + i) = 0;
      } /* for i */
   } /* else if */
   
   /* Compute Alpha */
   for( i = 0; i < lenData; i++) {
      
      *(tmpMax + i) = -INF - 1;
      
      for( j = 0; j < nrSt; j++) {
         /* Initialize Gamma */
         for( k = 0; k < nrSt; k++) {
            *(Gamma + k) = -INF;
         } /* for k */ 
         
         /* Avoid error of log of zero or negative */
         if( (1 + exp( *(LIn + i))) < SMALL) {
            X = -INF;
         }
         else {
            X = log( 1 + exp( *(LIn + i)));
         }
         
         /* Compute Gamma */
         *(Gamma + (int) *(prevSt + j) - 1) = (-*(x + i*2) 
                                              + *(x + i*2 + 1)**(prev + nrSt + j)) - X;                                          
         *(Gamma + (int) *(prevSt + j + nrSt) - 1) = (*(x + i*2) 
                                                   + *(x + i*2 + 1)**(prev + 3*nrSt + j)) 
                                                   + *(LIn + i) - X;
         
         sum = 0;
         for( k = 0; k < nrSt; k++) {
            sum += exp( *(Gamma + k) + *(*(Alpha + i) + k));
         } /* for k */
         
         if( sum < SMALL) {
            *(*(Alpha + i + 1) + j) = -INF;
         } /* if */
         else {
            *(*(Alpha + i + 1) + j) = log( sum);
         } /* else */
                
         /* Find maximum */
         if( *(*(Alpha + i + 1) + j) > *(tmpMax + i)) {
            *(tmpMax + i) = *(*(Alpha + i + 1) + j);
         } /* if */
      } /* for j */
      
      for( j = 0; j < nrSt; j++) {
         *(*(Alpha + i + 1) + j) = *(*(Alpha + i + 1) + j) - *(tmpMax + i);       
      } /* for j */

   } /* for i */
   
   /* Compute Beta */
   for( i = lenData - 2; i >= 0; i--) {
      
      for(j = 0; j < nrSt; j++) {
         /* Initialize Gamma */
         for( k = 0; k < nrSt; k++) {
            *(Gamma + k) = -INF;
         } /* for k */
         
         /* Avoid error of log of zero or negative */
         if( (1 + exp( *(LIn + i + 1))) < SMALL) {
            X = -INF;
         }
         else {
            X = log( 1 + exp( *(LIn + i + 1)));
         }
      
         /* Compute Gamma */
         *(Gamma + (int) *(nextSt + j) - 1) = (-*(x + i*2 + 2) 
                                              + *(x + i*2 + 3)**(next + nrSt + j)) - X;
         *(Gamma + (int) *(nextSt + j + nrSt) - 1) = (*(x + i*2 + 2) 
                                                   + *(x + i*2 + 3)**(next + 3*nrSt + j)) 
                                                   + *(LIn + i + 1) - X;
         
         sum = 0;
         for( k = 0; k < nrSt; k++) {
            sum += exp( *(Gamma + k) + *(*(Beta + i + 1) + k));
         } /* for k */
         
         if( sum < SMALL) {
            *(*(Beta + i) + j) = -INF;
         } /* if */
         else {
            *(*(Beta + i) + j) = log( sum);
         } /* else */  
      
      } /* for j */
      
      for( j = 0; j < nrSt; j++) {
         *(*(Beta + i) + j) = *(*(Beta + i) + j) - *(tmpMax + i);         
      } /* for j */
      
   } /* for i */
   
   /* Compute the log-likelihood ratio of the symbols */
   for( i = 0; i < lenData; i++) {
    
      tmp0 = 0; tmp1 = 0;
      for( j = 0; j < nrSt; j++) {
         
         /* Avoid error of log of zero or negative */
         if( 1 + exp( *(LIn + i)) < SMALL) {
            X = -INF;
         }
         else {
            X = log( 1 + exp( *(LIn + i)));
         }

         /* Compute Gamma */
         Gamma0 = ( -*(x + 2*i) + *(x + 2*i + 1) * *(prev + nrSt + j)) - X;
                  
         Gamma1 = ( *(x + 2*i) + *(x + 2*i + 1) * *(prev + 3*nrSt + j)) 
                  + *(LIn + i) - X;
         
         tmp0 += exp( Gamma0 + *(*(Alpha + i) + (int) *(prevSt + j) - 1) 
                + *(*(Beta + i) + j));
         tmp1 += exp( Gamma1 + *(*(Alpha + i) + (int) *(prevSt + j + nrSt) - 1) 
                + *(*(Beta + i) + j));              
      } /* for j */
      
      /* Avoid error of log of zero or negative */      
      if( tmp0 < SMALL) {
         X = -INF;
      }
      else {
         X = log( tmp0);
      }
      if( tmp1 < SMALL) {
         Y = -INF;        
      }
      else {
         Y = log( tmp1);
      }

      /* The log-likelihood ratio output */
      *(LOut + i) = Y - X;
   } /* for i */

} /* LogMap() */


/* ----------- Gateway function ---------- */
void mexFunction( int nlhs, mxArray *plhs[],
                    int nrhs, const mxArray *prhs[] )

/* input : x = prhs[0], received vector
 *         LIn = prhs[1], priori information for the current decoder
 *         decIndx = prhs[2], decoder index
 *         next = prhs[3], next trellis output
 *         nextSt = prhs[4], next state
 *         prev = prhs[5], previous trellis output
 *         prevSt = prhs[6], previous state    
 *
 * output: LOut = plhs[0], log-likelihood ratio of the symbols. */  
{
   
   double *LOut, *x, *LIn, *next, *nextSt, *prev, *prevSt; 
   int decIndx, mrows, len;
   
   if( nrhs != 7) {
      mexErrMsgTxt( "Missing input argument(s)!");
   }
   
   /* Create pointer to the inputs */
   x = mxGetPr( prhs[0]);
   LIn = mxGetPr( prhs[1]);   
   decIndx = mxGetScalar( prhs[2]);   
   next = mxGetPr( prhs[3]);   
   nextSt = mxGetPr( prhs[4]);   
   prev = mxGetPr( prhs[5]);   
   prevSt = mxGetPr( prhs[6]);
   
   /* Get the number of row of the vector input next */
   mrows = mxGetM( prhs[3]);

   /* Get the length of the vector input x */
   len = mxGetN( prhs[0]);
        
   /* Set the output pointer to the output vector */
   plhs[0] = mxCreateDoubleMatrix( 1, (int) len/2, mxREAL);
   
   /* Create a pointer to a copy ot the output vector LOut */
   LOut = mxGetPr( plhs[0]);
    
   /* Call the main subroutine */
   logMap( LOut, x, LIn, decIndx, next, nextSt, prev, prevSt, len, mrows);  
}