trellismex.c

802.16eCTC的logmap算法仿真

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

/* Prototypes */
int nextStToInt(int msb, int m, int st, int I);
int *binState( int st, int *state, int m);
void stTrellis( int *result, double *G, int m, int st, int bit);


int nextStToInt( int msb, int m, int st, int I)
/* Next state into integer.
 * Convert m bits vector of RSC's state st into integer of the next state.
 * msb is the feedback bit */ 
{
   int i, state[8] = {0,0,0,0,0,0,0,0};
   
   for( i = 0; i < m - 1; i++) {
      I = I + *(binState( st, state, m) + i) * pow( 2, m - i - 2);
   } /* for */
   I = I + msb * pow( 2, m - 1);
   return I;
}
 

int *binState( int st, int *state, int m)
/* Convert integer into m bits binary vector
 * Maximum up to 8 bit = 256 */ 
{   
   int i;
   div_t N;
   
   for( i = m - 1; i >= 0; i--) {
      N = div( st, 2);
      *(state + i) = N.rem;
      st = N.quot;            
   }
   return state;
}


void stTrellis( int *result, double *G, int m, int st, int bit)
/* Find next state, previous state, next output and previous output */ 
{   
   int i, N1, N2, state[8] = {0,0,0,0,0,0,0,0};
   div_t N;
      
   N1 = 0;
   /* The upper matrix generator */
   for( i = 0; i < m; i++) {
      N = div( (N1 + *(binState( st, state, m) + i) * *(G + 2*i + 2)), 2);
      N1 = N.rem;
   } /* for */
   N1 = N.rem;
      
   /* The feedback */
   N = div( (bit + N1), 2);
   N1 = N.rem;   
   result[0] = N1;
 
   N2 = 0;
   /* The lower matrix generator */
   for( i = 0; i < m; i++) {
      N = div( (N2 + *(binState( st, state, m)+i) * *(G + 2*i + 3)), 2);
      N2 = N.rem;          
   } /* for */
   N2 = N.rem;
   N = div( (N2 + N1 * *(G + 1)), 2);
   N2 = N.rem;   
   result[1] = N2;   
}


void trellis( double *next, double *nextState, double *prev, 
               double *prevState, double *G, int m, int maxSt)
/* ----------- The main subroutine ---------- */
{
   int i, result[2];

   for( i = 0; i < maxSt; i++) {
      
   /* ---------- Input = "0" ---------- */      
      stTrellis( result, G, m, i, 0);
      
      /* Next state integer */
      *(nextState + i) = nextStToInt(*result, m, i, 0) + 1;
      
      /* Previous state integer */
      *(prevState + (int) *(nextState + i) - 1) = i + 1;
      
      /* Next and previous output, systematic
       * "0" is represented by -1 */
      *(next + i) = -1;
      *(prev + i) = -1;
      
      /* Next output non systematic - parity */         
      *(next + 1 * maxSt + i) = 2 * *(result + 1) - 1;
      
      /* Previous output non systematic - parity */
      *(prev + maxSt + (int) *(nextState + i)-1) = *(next + maxSt + i); 
      
      
   /* ---------- Input = "1" ---------- */      
      stTrellis( result, G, m, i, 1);
      
      /* Next state integer */
      *(nextState + maxSt + i) = nextStToInt(*result, m, i, 0) + 1;
      
      /* Previous state integer */
      *(prevState + maxSt + (int) *(nextState + maxSt + i) - 1) = i + 1;
      
      /* Next and previous output, systematic */
      *(next + 2 * maxSt + i) = 1;
      *(prev + 2 * maxSt + i) = 1;

      /* Next output non systematic - parity */
      *(next + 3*maxSt + i) = 2 * *(result + 1) - 1;
      
      /* Previous output non systematic - parity */
      *(prev + 3*maxSt + (int) *(nextState + i)-1) = *(next + 3*maxSt + i);
            
   } /* for */         
}


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

/* input : G = prhs[0], generator matrix
 *
 * output: next = plhs[0], next output
 *         nextState = plhs[1], next state
 *         prev = plhs[2], previous output
 *         prevState = plhs[3], previous state. */
{                     
   double *G, *next, *nextState, *prev, *prevState;
   int ncols, m, maxState;
   
                     
   /* Create pointer to the input matrix G */
   G = mxGetPr( prhs[0]);
   
   /* Get the # of column of the matrix generator G */
   ncols = mxGetN( prhs[0]);
   
   /* # of memory */
   m = ncols - 1;
      
   /* # of state */
   maxState = (int) pow( 2,m);
   
   /* Set the output pointer to the output vector */
   plhs[0] = mxCreateDoubleMatrix( maxState, 4, mxREAL);
   plhs[1] = mxCreateDoubleMatrix( maxState, 2, mxREAL);
   plhs[2] = mxCreateDoubleMatrix( maxState, 4, mxREAL);
   plhs[3] = mxCreateDoubleMatrix( maxState, 2, mxREAL);
   
   /* Create a pointer to a copy ot the output matrix */
   next = mxGetPr( plhs[0]);
   nextState = mxGetPr( plhs[1]);
   prev = mxGetPr( plhs[2]);
   prevState = mxGetPr( plhs[3]);
   
   /* Call the main subroutine */
   trellis( next, nextState, prev, prevState, G, m, maxState);   
}