sovamex.c

802.16eCTC的logmap算法仿真

/* MATLAB's C-MEX turbo code soft output viterbi algorithm decoder
 *
 * Copyright Bagawan S. Nugroho, 2006 
 *********************************************/

#include "mex.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 sova( 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, st0, st1, tmpSt, bit;
   int prevBit[1000][64], mlSt[1000], est[1000];
   int Delta = 30;
   double M0, M1, LLR, tmpMax;
   double metric[1000][64], MDiff[1000][64], sym0[2], sym1[2];
   
   /* Length of the decoded information */
   lenData = (int) len/2;
   
   /* Initialize the path metric */
   for( i = 0; i < lenData; i++) {
      for( j = 0; j < nrSt; j++) {
         *(*(metric + i) + j) = -INF;
      } /* for i */
   } /* for j */
   *(*metric) = 0; 
   
   for( i = 0; i < lenData; i++) {
      for( j = 0; j < nrSt; j++) {
         *sym0 = *(prev + j);
         *(sym0 + 1) = *(prev + nrSt + j);
         *sym1 = *(prev + 2*nrSt + j);
         *(sym1 + 1) = *(prev + 3*nrSt + j);
         
         st0 = (int) *(prevSt + j) - 1;
         st1 = (int) *(prevSt + nrSt + j) - 1;
                  
         M0 = (*(x + 2*i) * *sym0) + (*(x + 2*i + 1) * *(sym0 + 1)) 
               - *(LIn + i)/2 + *(*(metric + i) + st0); 
         M1 = (*(x + 2*i) * *sym1) + (*(x + 2*i + 1) * *(sym1 + 1)) 
               + *(LIn + i)/2 + *(*(metric + i) + st1);
         
         if( M0 > M1) {
            *(*(metric + i + 1) + j) = M0;
            *(*(MDiff + i + 1) + j) = M0 - M1;
            *(*(prevBit + i + 1) + j) = 0;          
         } /* if */
         else {
            *(*(metric + i + 1) + j) = M1;
            *(*(MDiff + i + 1) + j) = M1 - M0;
            *(*(prevBit + i + 1) + j) = 1;        
         } /* else */      
      } /* for j */    
   } /* for i */
   
   /* Decoder 1, trace back from all zero states.
    * Decoder 2, trace back from the most likely states */
   if( decIndx == 1) {
      *(mlSt + lenData) = 1;    
   } /* if */
   else if( decIndx == 2) {
      tmpMax = -INF - 1;
      for( j = 0; j < nrSt; j++) {         
         if( *(*(metric + lenData) + j) > tmpMax) {
            *(mlSt + lenData) = j + 1;
            tmpMax = *(*(metric + lenData) + j);
         } /* if */         
      } /* for j */         
   } /* else if */

   /* Trace back to get the estimated bits and the most likely path */
   for( i = lenData - 1; i >= 0; i--) {
      *(est + i) = *(*(prevBit + i + 1) + *(mlSt + i + 1) - 1);
      *(mlSt + i) = (int) *(prevSt + *(mlSt + i + 1) + (*(est + i) * nrSt) - 1);    
   } /* for i */
   
   /* Find the minimum of delta */
   for( i = 0; i < lenData; i++) {
      LLR = INF;
      for( k = 0; k <= Delta; k++) {
         if( i + k < lenData) {
            bit = 1 - *(est + i + k);
            
            tmpSt = (int) *(prevSt + *(mlSt + i + k + 1) + bit * nrSt - 1);
            
            for( j = k - 1; j >= 0; j--) { 
               bit = *(*(prevBit + i + j + 1) + tmpSt - 1);
               tmpSt = (int) *(prevSt + tmpSt + bit * nrSt - 1);
            } /* for j */
            
            if( bit != *(est + i)) { 
               if( LLR < *(*(MDiff + i + k + 1) + *(mlSt + i + k + 1) - 1)) {
                  LLR = LLR;
               } /* if */
               else {
                  LLR = *(*(MDiff + i + k + 1) + *(mlSt + i + k + 1) - 1);
               } /* else */
            } /* if */            
         } /* if */
      } /* for k */
      
      *(LOut + i) = (2**(est + i) - 1) * LLR;
      
   } /* for i */  
   
} /* sova() */

/* ----------- 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 # 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 */
   sova( LOut, x, LIn, decIndx, next, nextSt, prev, prevSt, len, mrows);
}