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