convviterbi.c

卷积码的编码及解码代码

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

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

#define  INF 1e10 

/* ----------- The main function ----------- */
void convViterbi( double *x, double *y, short R, short L, double *c0, double *c1, 
                     short nrSt, short D)
/* x, input matrix
 * R, denominator of the code rate (eq. 1/R)
 * L, output length
 * c0, all possible codewords when input is '0'
 * c1, all possible codewords when input is '1'
 * nrSt, number of state
 * D, number of memory = constraint length - 1 
 * 
 * y, output vector. */
                     
{
   short i, j, k, halfSt, state, path[256][1024];
   double tmp, distEven, distOdd, metric[256], nextMetric[256], tmpOut[1024];
   
   /* Half of the number of state */
   halfSt = (short) nrSt/2;
   
   /* Initialize the metric vector */
   for( i = 0; i < nrSt; i++) {
      *(metric + i) = INF;
   }
   *(metric) = 0; 
   
   /* Iteration over L - K bits */
   for( i = 0; i < L - D; i++) {        
      /* Iteration over states */      
      for( j = 0; j < halfSt; j++) {         
         /* Find the distance between the input and codewords of '0' */
         distEven = 0; distOdd = 0;
         for( k = 0; k < R; k++) {                     
            /* Even states codewords when input is '0' */
            tmp = *(x + k + i*R)**(c0 + k + 2*j*R);
            if( tmp < 0) {
               tmp = -tmp;                       
            }
            else {
               tmp = 0;
            }
            distEven += tmp;
                              
            /* Odd states codewords when input is '0' */
            tmp = *(x + k + i*R)**(c0 + k + (2*j + 1)*R);
            if( tmp < 0) {
               tmp = -tmp;                       
            }
            else {
               tmp = 0;
            }
            distOdd += tmp;
         } /* for k */
         /* Accumulate the distance with the previous branch metric */            
         distEven = *(metric + 2*j) + distEven;
         distOdd = *(metric + 2*j + 1) + distOdd;                             
         /* Find the next metric and the surviving paths */
         if( distEven <= distOdd) {
            /* An even state branch is selected, update the next metric */
            *(nextMetric + j) = distEven;
            *(*(path + j) + i) = 2*j;
         }
         else {
            /* An odd state branch is selected, update the next metric */
            *(nextMetric + j) = distOdd;
            *(*(path + j) + i) = 2*j + 1;
         }
         
         /* Find the distance between the input and codewords of '1' */               
         distEven = 0; distOdd = 0;
         for( k = 0; k < R; k++) {
            /* Even states codewords when input is '1' */
            tmp = *(x + k + i*R)**(c1 + k + 2*j*R);
            if( tmp < 0) {
               tmp = -tmp;                       
            }
            else {
               tmp = 0;
            }
            distEven += tmp;
               
            /* Odd states codewords when input is '1' */
            tmp = *(x + k + i*R)**(c1 + k + (2*j + 1)*R);
            if( tmp < 0) {
               tmp = -tmp;                       
            }
            else {
               tmp = 0;
            }
            distOdd += tmp;
         } /* for k */
         /* Accumulate the distance with the previous branch metric */      
         distEven = *(metric + 2*j) + distEven;
         distOdd = *(metric + 2*j + 1) + distOdd;           
         /* Find the next metric and the surviving paths */
         if( distEven <= distOdd) {
            /* An even state branch is selected, update the next metric */
            *(nextMetric + j + halfSt) = distEven;
            *(*(path + j + halfSt) + i) = 2*j;
         }
         else {
            /* An odd state branch is selected, update the next metric */
            *(nextMetric + j + halfSt) = distOdd;
            *(*(path + j + halfSt) + i) = 2*j + 1;
         }
                          
      } /* for j */
      /* Update the branch metric vector */
      for( j = 0; j < nrSt; j++) {
         *(metric + j) = *(nextMetric + j);
      } /* for j */         
   } /* for i */
   
   /* Iteration over L - D to L bits, which known has codewords '0' */
   for( i = L - D; i < L; i++) {     
      /* Iteration over states */
      for( j = 0; j < halfSt; j++) {
         /* Find the distance between the input and codewords of '0' */
         distEven = 0; distOdd = 0;                
         for( k = 0; k < R; k++) {         
            /* Even states codewords when input is '0' */
            tmp = *(x + k + i*R)**(c0 + k + 2*j*R);
            if( tmp < 0) {
               tmp = -tmp;                       
            }
            else {
               tmp = 0;
            }
            distEven += tmp;
               
            /* Odd states codewords when input is '0' */
            tmp = *(x + k + i*R)**(c0 + k + (2*j + 1)*R);
            if( tmp < 0) {
               tmp = -tmp;                       
            }
            else {
               tmp = 0;
            }
            distOdd += tmp;
         } /* for k */
         /* Accumulate the distance with the previous branch metric */            
         distEven = *(metric + 2*j) + distEven;
         distOdd = *(metric + 2*j + 1) + distOdd;              
         /* Find the next metric and the surviving paths */
         if( distEven <= distOdd) {
            /* An even state branch is selected, update the next metric */
            *(nextMetric + j) = distEven;
            *(*(path + j) + i) = 2*j;
         }
         else {
            /* An odd state branch is selected, update the next metric */
            *(nextMetric + j) = distOdd;
            *(*(path + j) + i) = 2*j + 1;
         }
         /* For codewords of '1', any value does not matter */
         *(nextMetric + j + halfSt) = INF;
         *(*(path + j + halfSt) + i) = 0;                                  
      } /* for j */
      /* Update the branch metric vector */
      for( j = 0; j < nrSt; j++) {
         *(metric + j) = *(nextMetric + j);
      }         
   } /* for i */
   
   /* Trace the trellis back */ 
   /* The upper half states represent '0' and the opposite represent '1' */
   state = 0;
   for( i = L - 1; i >= 0; i--) {
      if( state < halfSt) {
         *(tmpOut + i) = 0;
      }
      else {
         *(tmpOut + i) = 1;
      }
      state = *(*(path + state) + i);
   }
   /* Remove the padding tail */
   for( i = 0; i < L - D; i++) {
      *(y + i) = *(tmpOut + i);
   }

} /* convViterbi */

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

/* input:  x = prhs[0], input matrix with size of R X L
 *         c0 = prhs[1], codewords '0' matrix (-1/1)
 *         c1 = prhs[2], codewords '1' matrix (-1/1)
 *         K = prhs[3], constraint length
 *     
 * output: y = plhs[0], output vector (0/1). */
{    
   double *x, *y, *c0, *c1, K;
   short R, L, nrSt;
   
   /* Catch the number of input argument's error if any */
   if( nrhs != 4) {
      mexErrMsgTxt( "Missing or wrong input argument(s)!");
   }
        
   /* Create pointer to the input matrix x */   
   x = mxGetPr( prhs[0]);
   
   /* Get the dimensions of the matrix x */
   R = mxGetM( prhs[0]);
   L = mxGetN( prhs[0]);   
   
   /* Create pointer to the codewords matrix c0 */
   c0 = mxGetPr( prhs[1]);
   
   /* Create pointer to the codewords matrix c1 */
   c1 = mxGetPr( prhs[2]);
   
   /* The size of c0 and c1 should be R X nrSt */   
   /* Get the length of c0 = c1 as the nrSt*/
   nrSt = mxGetN( prhs[1]);
   
   /* Create pointer to the constraint length K */
   K = mxGetScalar( prhs[3]);

   /* Set the output pointer to the output vector */   
   plhs[0] = mxCreateDoubleMatrix( 1, L - K + 1, mxREAL);
    
   /* Create a pointer to a copy ot the output vector y */
   y = mxGetPr( plhs[0]);
    
   /* Call the main subroutine */
   convViterbi( x, y, R, L, c0, c1, nrSt, (short) K - 1);            
}