/* MATLAB's C-MEX space-time trellis decoder for QPSK
 *
 * Copyright Bagawan S. Nugroho, 2006 
 *********************************************/

#include "mex.h"
#include "matrix.h"
#include <stdlib.h>
/* m-Ary, M = 4 for QPSK */
#define M            4
/* Number of incoming branch */
#define branch       4
/* Maximum number of state */
#define ST           256
/* Maximum length of data */
#define LEN          2000

/* ----------- The main function ---------- */
void stViterbi( double *sRe, double *sIm, double nrSt, int inLen, double *st1, double *st2, 
                  double *bin, int row, int col, double *chRe, double *chIm, double *y)
/* sRe, the real part of the received signal s
 * sIm, the imaginary part of the received signal s
 * nrSt, number of state
 * inLen, length of the received signal s
 * st1, current states trellis
 * st2, next states trellis
 * bin, decimal to binary map
 * row, number of row of matrix bin
 * col, number of col of matrix bin
 * chRe, the real part of the fading channel
 * chIm, the imaginary part of the fading channel
 *
 * y, the decoder output (0/1). */
                                    
{
   int i, k, l, m, n, L, stMap1, stMap2, delay, prevSt[ST][LEN], xTmp[LEN], decSt[LEN];
   double dist, distRe, distIm, minDist, metric[ST], nextMetric[ST]; 
   
   div_t Z;
   
   /* Some initializations... */
   int constelRe[4] = {1, 0, -1, 0};
   int constelIm[4] = {0, 1, 0, -1};
   
   Z = div( col, 2);
   delay = Z.quot + Z.rem;
   
   L = (int) nrSt/M;
    
   for( i = 0; i < nrSt; i++) {
      *(metric + i) = 0;
      for( k = 0; k < inLen; k++) {
         *(*(prevSt + i) + k) = 0;
      } /* for j */
   } /* for i */
   
   /* Viterbi decoding; find the surviving paths over time interval inLen */
   for( k = 1; k < inLen; k++) {
      for( l = 0; l < L; l++) {
         for( m = 0; m < M; m++) {
            minDist = 1e10;
            for( n = 0; n < branch; n++) {
               /* stMap1 = current states map; stMap2 = next states map */ 
               stMap1 = (int) *(st1 + L*n + l);
               stMap2 = (int) *(st2 + L*n + l + (int) nrSt*m);
               
               /* Find the eucledean distance of the received vector over all possible branches */
               
               /* The real part */
               distRe = (*chRe**(constelRe + stMap1) - *chIm**(constelIm + stMap1))
                  + (*(chRe + 1)**(constelRe + stMap2) - *(chIm + 1)**(constelIm + stMap2))
                  - *(sRe + k);
               /* The imaginary part */
               distIm = (*chIm**(constelRe + stMap1) + *chRe**(constelIm + stMap1))
                  + (*(chIm + 1)**(constelRe + stMap2) + *(chRe + 1)**(constelIm + stMap2)) 
                  - *(sIm + k);
               /* The magnitude */
               dist = distRe*distRe + distIm*distIm;
               
               /* Determine the minimum eucledean distance and the suviving paths */
               if( dist + *(metric + L*n + l) < minDist) {
                  minDist = dist + *(metric + L*n + l);
                  *(nextMetric + m + 4*l) = minDist;
                  *(*(prevSt + m + 4*l) + k) = L*n + l;                  
               } /* if */
            } /* for n */
         } /* for m */
      } /* for l */
      /* Update the branch metric vector */
      for( i = 0; i < nrSt; i++) {
         *(metric + i) = *(nextMetric + i);
      } /* for i */
   } /* for k */
   
   /* Trace back to find the estimated states */
   *(xTmp + inLen - 1) = 0;
   for( k = inLen - 1; k >= 0; k--) {      
      Z = div(*(xTmp + k + 1), 4);      
      *(xTmp + k) = *(*(prevSt + *(xTmp + k + 1)) + k);
      *(decSt + k) = M**(xTmp + k) + Z.rem;      
   } /* for i */

   /* Find the estimated bits sequence y */
   for( i = 0; i < inLen - delay; i++) {
      /* Even index sequence */
      *(y + 2*i) = *(bin + *(decSt + i) + row*(col - 1));
      /* Odd index sequence */
      *(y + 2*i + 1) = *(bin + *(decSt + i) + row*(col - 2));
   } /* for i */

} /* main */

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

/* input: s    = prhs[0], input vector; complex number   
 *        st1  = prhs[1], current state matrix
 *        st2  = prhs[2], next state matrix
 *        nrSt = prhs[3], number of state
 *        bin  = prhs[4], decimal to binary map
 *        ch   = prhs[5], fading channel; complex number    
 *
 * output: y   = plhs[0], output vector (0/1). */
{    
   double *sRe, *sIm, nrSt, *st1, *st2, *bin, *chRe, *chIm, *y;
   int x, inLen, row, col;
   
   div_t Z;
   
   if( nrhs != 6) {
      mexErrMsgTxt( "Missing or wrong input argument(s)!");
   }
        
   /* Create pointer to the complex input vector s */
   sRe = mxGetPr( prhs[0]);
   sIm = mxGetPi( prhs[0]);

   /* Get the length of the vector input s */
   inLen = mxGetN( prhs[0]);
          
   /* Create pointer to the input matrix st1 */
   st1 = mxGetPr( prhs[1]);
   
   /* Create pointer to the input matrix st2 */
   st2 = mxGetPr( prhs[2]);
   
   /* Get the scalar input nrSt */
   nrSt = mxGetScalar( prhs[3]);
   
   /* Get the pointer to the input matrix bin */
   bin = mxGetPr( prhs[4]);
   
   /* Get dimension of the matrix bin */
   row = mxGetM( prhs[4]);
   col = mxGetN( prhs[4]);
   
   /* Get the pointer to the complex input vector ch */
   chRe = mxGetPr( prhs[5]);
   chIm = mxGetPi( prhs[5]);
   
   /* Create output the matrix y with length of x */
   Z = div( col, 2);
   x = Z.quot + Z.rem;   
   x = 2*(inLen - x);
   
   plhs[0] = mxCreateDoubleMatrix( 1, x, mxREAL);
    
   /* Create a pointer to a copy ot the output matrix */
   y = mxGetPr( plhs[0]);
    
   /* Call the main subroutine */
   stViterbi( sRe, sIm, nrSt, inLen, st1, st2, bin, row, col, chRe, chIm, y);   
}