siso_app_cc.cpp

通信中常用的卷积码信道译码源码程序

/*
 *-----------------------------------------------------------------------------
 * Maximum a posteriori decoder with soft-inputs and soft-outputs. 
 *
 * Calculates Log-likelihood ratios of received bits. A priori information 
 * of information and coded bits from the other decoder and soft democulator 
 * accepted. Soft inputs, soft outputs. Log-domain calculations for numerical 
 * stability. No approximations used in the calculation of logarithm
 * of sum of exponents.
 *
 * Note: PCCC-schemes shouldn't use a priori LLRs: log(P(c;O)). See references.
 * 
 * Algorithm reference:
 *   Benedetto, S., Divsalar, D., Montorsi, G. and Pollara, F.,
 *   "A Soft-Input Soft-Output APP Module for Iterative Decoding of 
 *    Concatenated Codes". IEEE Communications Letters, Vol.1, No.1, 1997.
 *   Benedetto, S., Divsalar, D., Montorsi, G. and Pollara, F.,
 *   "A Soft-Input Soft-Output Maximum A Posterioiri (MAP) Module to Decode
 *    Parallel and Serial Concatenated Codes". TDA Progress Report 42-127.
 * 
 * Restrictions: Two 1/2-rate RSCs as CCs assumed.
 *
 * 
 * Format:
 * -------
 * 
 * [LuO,LcO] = siso_app(LuI, LcI, Trellis, BitsLut, dec, LcOFlag);
 *
 * Inputs:
 * -------
 *  LuI - log(P(u;I)), a priori log-likelihoods of the information bits 
 *                     from other component decoder
 *  LcI - log(P(c;I)), a priori log-likelihoods of the coded bits
 *                     from demodulator (PCCC) or other component decoder(SCCC)
 *  Trellis - Structure containing necessary information of CC trellises
 *  BitsLut - Bit look-up table for edge output symbols
 *  dec - Decoder 1 or 2
 *  LcOFlag - Calculate likelihoods log(P(c;O)) or not
 *
 * Outputs:
 * --------
 *  LuO - log(P(u;O)), extrinsic bit information of the information bits
 *  LcO - log(P(c;O)), extrinsic bit information of the coded bits
 * 
 * Author: MVe, [e-mail]
 * Date:   14.10.2002
 *-----------------------------------------------------------------------------
 */

#include <cmath>            // basic mathematical functions
#include <blitz/array.h>    // matrix operations
#include "mex.h"            // matlab interface

using namespace std;
using namespace blitz;


/* Define different types of arrays */

// vector with type 'int' elements
typedef Array<int,1> IArray1D;
// matrix with type 'int' elements
typedef Array<int,2> IArray2D;
// 3-D Matrix with type 'int' elements 
typedef Array<int,3> IArray3D;

// vector with type 'double' elements
typedef Array<double,1> DArray1D;
// matrix with type 'double' elements
typedef Array<double,2> DArray2D;
// 3-D Matrix with type 'double' elements 
typedef Array<double,3> DArray3D;

double maxx(double,double);	// function, maxx(x,y) = log(exp(x)+exp(y))


/**************************************************
 * C++ function that implements the MAP algorithm *
 **************************************************/

static void siso_app_f(double *LuI_in, int L_total, double *LcI_in, 
		       double *edge_in, int Nstates, int Ninputs, 
		       double *BitsLut_in, int Noutputs, int Nbitsout, 
		       int dec, 
		       double *output1, double *output2)


{

  /* Local variables */ 
  int i, j, k, e, b, M, N;
  double log_norma_alpha, log_norma_beta, log_norma_LuO, scale;
  const double INF = 1E+100;
  const int e_S=1, e_E=2, e_U=3, e_C1=4, e_C2=5,E_ENTRIES=5;

  // local matrices 
  IArray2D binC(Range(1,1),Range(1,Nbitsout));
  IArray2D edge(Range(1,Nstates*Ninputs),Range(1,E_ENTRIES));
  IArray2D BitsLut(Range(1,Noutputs),Range(1,Nbitsout));
  DArray1D LuO(Range(1,L_total));
  DArray1D LuO_1(Range(1,L_total));
  DArray1D LuO_0(Range(1,L_total));
  DArray1D LcO1(Range(1,L_total));
  DArray1D LcO1_0(Range(1,L_total));
  DArray1D LcO1_1(Range(1,L_total));
  DArray1D LcO2(Range(1,L_total));
  DArray1D LcO2_0(Range(1,L_total));
  DArray1D LcO2_1(Range(1,L_total));
  DArray2D LcO(Range(1,Nbitsout),Range(1,L_total));
  DArray2D alpha(Range(1,Nstates),Range(1,L_total+1));
  DArray2D beta(Range(1,Nstates),Range(1,L_total+1));


  // input matrices 
  DArray2D LuI(LuI_in, shape(1,L_total), neverDeleteData, 
               fortranArray);
  DArray2D LcI(LcI_in, shape(Nbitsout,L_total), neverDeleteData, 
               fortranArray);
  DArray2D edge_tmp(edge_in, shape(Nstates*Ninputs,E_ENTRIES), 
		    neverDeleteData, fortranArray);
  DArray2D BitsLut_tmp(BitsLut_in, shape(Noutputs,Nbitsout), neverDeleteData, 
                       fortranArray);


  /* Initialise local matrices and variables */

  alpha = -INF;
  beta = -INF;
  LuO = -INF;
  LuO_0 = -INF;
  LuO_1 = -INF;
  LcO = -INF;
  LcO1 = -INF;
  LcO1_0 = -INF;
  LcO1_1 = -INF;
  LcO2 = -INF;
  LcO2_0 = -INF;
  LcO2_1 = -INF;

  for (i = 1; i <= Nstates*Ninputs; i++){
    for (j = 1; j <= E_ENTRIES; j++){
      edge(i,j) =(int) edge_tmp(i,j);
    }
  }

  for (i = 1; i <= Noutputs; i++){
    for (j = 1; j <= Nbitsout; j++)
      BitsLut(i,j) =(int) BitsLut_tmp(i,j);
  }


  /*************************************
   ** -- alpha (forward recursion) -- ** 
   *************************************
   *
   * Note:
   *
   * L_e(u(n)) = log(P(u(n)==1) / P(u(n)==0)) = L(u(n)==1) - L(u(n)==0)
   * 
   *  ==> L(u(n)==1) = L_e(u(n)==1) - log(1 + exp(L_e(u(n)==1)))
   *  ==> L(u(n)==0) = - log(1 + exp(L_e(u(n)==0)))
   *    ==> L(u(n)) = u(n)*L_e(u(n)) - log(1 + exp(L_e(u(n))))  | u(n) = {0,1}
   *
   * L_e(u(n)) - extrinsic information from other component decoder
   *
   */
  
  /* Format of 'alpha':
   *   row = state
   *   column = time
   */

  alpha(1,1)=0.0;

  for (k=2; k<=L_total+1;k++){

    for (e=1; e<=Nstates*Ninputs;e++){
      
      alpha(edge(e,e_E),k) = 
	maxx(alpha(edge(e,e_E),k), 
	     alpha(edge(e,e_S),k-1)+ 
	     edge(e,e_U)*LuI(1,k-1)-maxx(0.0,LuI(1,k-1))+ 
	     edge(e,e_C1)*LcI(1,k-1)-maxx(0.0,LcI(1,k-1))+ 
	     edge(e,e_C2)*LcI(2,k-1)-maxx(0.0,LcI(2,k-1)));
      
    }
  }



  /*************************************
   ** -- beta (backward recursion) -- **
   *************************************/
  
  /* Format of 'beta':
   *  row = state
   *  column = time
   */
  
  
  if (dec==1)
    beta(1,L_total+1)=0.0;
  else
    beta(Range(1,Nstates),L_total+1)=alpha(Range(1,Nstates),L_total+1);
 
  for (k=L_total;k>=1;k--){

    for (e=1; e<=Nstates*Ninputs;e++){

      beta(edge(e,e_S),k) = 
	  maxx(beta(edge(e,e_S),k), 
	       beta(edge(e,e_E),k+1)+ 
	       edge(e,e_U)*LuI(1,k)-maxx(0.0,LuI(1,k))+ 
	       edge(e,e_C1)*LcI(1,k)-maxx(0.0,LcI(1,k))+ 
	       edge(e,e_C2)*LcI(2,k)-maxx(0.0,LcI(2,k)));

    }
  }


  // ******************* LuO(k,l)=log(P(u(k,l);O)) ************************

  for (k=2;k<=L_total+1;k++){

    for (e=1; e<=Nstates*Ninputs;e++){  

      if (edge(e,e_U))
	{
	  LuO_1(k-1) = 
	    maxx(LuO_1(k-1), 
		 alpha(edge(e,e_S),k-1)+ 
		 edge(e,e_C1)*LcI(1,k-1)-maxx(0.0,LcI(1,k-1))+ 
		 edge(e,e_C2)*LcI(2,k-1)-maxx(0.0,LcI(2,k-1))+ 
		 beta(edge(e,e_E),k));
	}        
      else
	{
	  LuO_0(k-1) = 
	    maxx(LuO_0(k-1), 
		 alpha(edge(e,e_S),k-1)+ 
		 edge(e,e_C1)*LcI(1,k-1)-maxx(0.0,LcI(1,k-1))+ 
		 edge(e,e_C2)*LcI(2,k-1)-maxx(0.0,LcI(2,k-1))+ 
		 beta(edge(e,e_E),k));        
	}
      LuO(k-1)=LuO_1(k-1)-LuO_0(k-1);  
    }
  }

  // *******************    LcO   ************************

//   for (t=2;t<=L_total+1;t++){

//     log_norma_LuO=-INF;

//     for (m=1;m<=Nstates;m++){
//       for (b=1;b<=Nbitsout;b++){

//         binC(1,Range(1,Nbitsout))=BitsLut(Xfi(m,b),Range(1,Nbitsout));

//         if (binC(1,1)==1)
//           {
//             LcO1_1(t-1)=maxx(LcO1_1(t-1),
//                              alpha(m,t-1)+beta(Sfi(m,b),t)
//                              +(b-1)*LuI(1,t-1)-maxx(0.0,LuI(1,t-1)) 
//                              +binC(1,2)*LcI(2,t-1)-maxx(0.0,LcI(2,t-1)));  
//           }

//         if (binC(1,1)==0)
//           {
//             LcO1_0(t-1)=maxx(LcO1_0(t-1),
//                              alpha(m,t-1)+beta(Sfi(m,b),t)
//                              +(b-1)*LuI(1,t-1)-maxx(0.0,LuI(1,t-1)) 
//                              +binC(1,2)*LcI(2,t-1)-maxx(0.0,LcI(2,t-1)));
//           }

//         if (binC(1,2)==1)
//           {
//             LcO2_1(t-1)=maxx(LcO2_1(t-1),
//                              alpha(m,t-1)+beta(Sfi(m,b),t)
//                              +(b-1)*LuI(1,t-1)-maxx(0.0,LuI(1,t-1)) 
//                              +binC(1,1)*LcI(1,t-1)-maxx(0.0,LcI(1,t-1)));  
//           }

//         if (binC(1,2)==0)
//           {
//             LcO2_0(t-1)=maxx(LcO2_0(t-1),
//                              alpha(m,t-1)+beta(Sfi(m,b),t)
//                              +(b-1)*LuI(1,t-1)-maxx(0.0,LuI(1,t-1)) 
//                              +binC(1,1)*LcI(1,t-1)-maxx(0.0,LcI(1,t-1)));
//           }
//       }
//     }
//     LcO(1,t-1)=LcO1_1(t-1)-LcO1_0(t-1);
//     LcO(2,t-1)=LcO2_1(t-1)-LcO2_0(t-1);
  
//   }


  // Write matrices to outputs

  N = L_total;

  for (j=1; j<=N; j++)
    *(output1 + (j-1)) = LuO(j);


  M = Nbitsout;

  for (i=1; i<=M; i++){
    for (j=1; j<=N; j++)
      *(output2 + (j-1)*M+(i-1)) = LcO(i,j);
  }



}

/*
 * maxx - 
 *
 * Function that takes two scalars x and y, and calculates:
 *
 * maxx_out = max(x,y) + ln(1 + exp(-abs(x-y)));
 *
 */

inline double maxx(double x, double y)
{

  double maxx_out,max_xy,exp_xy = 0.0;
  
  /* <cmath> provides functions:
   *     double exp(double x) -- Compute exponential of x 
   *     double log(double x) -- Compute (base-e logarithm) log(x)
   *     double fabs(double x) -- Compute absolute value of x
   */
  
  /* max_xy = max(x,y) */
  max_xy = (x>y) ? x : y;
  
  /* exp_xy = exp(-abs(x-y)) */
  exp_xy = exp(-fabs(x-y));
  
  /* Final result */
  maxx_out = max_xy + log(1+exp_xy);

  return(maxx_out);
}


/***********/
/* GATEWAY */
/***********/

/*
 * [LuO, LcO]=siso_app(LuI, LcI, edge, BitsLut, k0, dec, LcOFlag);
 */

void mexFunction(
		 int          nlhs,
		 mxArray      *plhs[],
		 int          nrhs,
		 const mxArray *prhs[]
		 )
{

  /* Inputs */
  double *LuI, *LcI, *edge, *BinLut;
  int dec, k0, LcOFlag;
  /* Outputs */
  double *output1,*output2;
  /* Local variables */
  int  n0,Nstates, Ninputs, Noutputs, L_total;
  double Nedges,Nstates_tmp;
  
  /* Check for proper number of arguments */  
  if (nrhs != 7)
    mexErrMsgTxt("'siso_app_CC' requires seven input arguments.");
  
  /* Assign pointers to input variables */
  LuI = mxGetPr(prhs[0]);            /* external information */
  LcI = mxGetPr(prhs[1]);            /* log-likelihoods */

  /* Trellis */
  edge = mxGetPr(prhs[2]);           /* trellis edge */
  BinLut = mxGetPr(prhs[3]);         /* lut for edge outputs*/
  k0 =(int) mxGetScalar(prhs[4]);    /* # of input bits / output */

  dec =(int) mxGetScalar(prhs[5]);
  LcOFlag =(int) mxGetScalar(prhs[6]);

  /* Size of the input matrix */
  n0 = mxGetM(prhs[1]);
  L_total = mxGetN(prhs[1]);

  /* Size of the look-up tables */
  Nedges =(double) mxGetM(prhs[2]);
  Nstates_tmp = Nedges/pow(2.0, (double) k0);
  Nstates = (int) Nstates_tmp;
  Ninputs =(int) pow(2.0, (double) k0);
  Noutputs = (int) pow(2.0, (double) n0);

  /* Assign pointers to output variables */
  plhs[0] = mxCreateDoubleMatrix(1,L_total,mxREAL);
  output1 = mxGetPr(plhs[0]);

  plhs[1] = mxCreateDoubleMatrix(n0,L_total,mxREAL);
  output2 = mxGetPr(plhs[1]);

  /* Function call */
  siso_app_f(LuI,L_total,LcI,edge,Nstates,Ninputs,BinLut,Noutputs,
	     n0,dec,output1,output2);


  return;

}