conv_dec_soft.cpp

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

/*
 *-----------------------------------------------------------------------------
 * This is a C++ mex-file for MATLAB
 *
 * Decode received 1-D data vector with trellis (soft inputs).
 *
 *
 * Format:
 * -------
 *
 * outputs = ...
 *   tcm_dec_CC(Data_received, NextState, Output_modulated, fading_gains, m)
 *
 *
 * Author: MVe
 * Date:   06.09.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;


/***************************************************
 * C++ function that implements the trellis search *
 ***************************************************/

static void conv_dec_CC(double *in_r, double *in_i, int Sinput1, int Sinput2, 
			double *nstate, int Nstates, int Ninputs, 
			double *TRout_r, double *TRout_i, 
			double *fad_r, double *fad_i, bool complex_ch, int m,
			double *output)

{

  /* Local variables */ 

  int i, j, t, state, input, Linput, decode_ptr, next_state;
  double faded_output_r=0.0, faded_output_i=0.0;
  double tmp_dist, euclid_metric, output_i, output_r;//, fad_tmp_r, fad_tmp_i;
  const double INF = 1E+300;

  // For handling matrix dimensions
  firstIndex i_ind;
  secondIndex j_ind;

  /* Allocate memory for arrays - note that indexing starts from one */

  // local matrices
  TinyVector<int,1> vec_tmp;
  IArray1D survivor_ptr(Range(1,Nstates));
  IArray1D ptr_tmp(Range(1,Nstates));
  DArray1D Metrics_dist(Range(1,Nstates));
  IArray2D TRstatepath(Range(1,Nstates), Range(1,Sinput2+1));
  DArray2D Metrics_dist_tmp(Range(1,Nstates), Range(1,Nstates));
  IArray3D TRstatepath_tmp(Range(1,Nstates), Range(1,Nstates), 
			   Range(1,Sinput2+1));

  // input matrices 
  DArray1D fading_r(fad_r, shape(Sinput2), neverDeleteData, 
		    fortranArray);
  DArray1D fading_i(fad_i, shape(Sinput2), neverDeleteData, 
		    fortranArray);
  DArray2D input_r(in_r, shape(Sinput1,Sinput2) , neverDeleteData, 
		   fortranArray);
  DArray2D input_i(in_i, shape(Sinput1,Sinput2), neverDeleteData, 
		   fortranArray);
  DArray2D TRoutput_r(TRout_r, shape(Nstates,Ninputs), neverDeleteData, 
			fortranArray);
  DArray2D TRoutput_i(TRout_i, shape(Nstates,Ninputs), neverDeleteData, 
		      fortranArray);
  DArray2D TRnextstate(nstate, shape(Nstates,Ninputs), neverDeleteData, 
		       fortranArray);


  /* Initialise local matrices */ 

  euclid_metric = 0.0;
  TRstatepath = 1;
  Metrics_dist_tmp = INF;
  Metrics_dist = INF;
  Metrics_dist(1) = 0.0;
  survivor_ptr = -1;
  ptr_tmp = -1;


  /*
   * -- Trellis calculations begin (soft inputs) --
   *
   */

  Linput = Sinput2;		// length of the input sequence

  // Loop over received symbols
  for (t = 1; t <= Linput; t++)
    {

      // Loop over allowed states
      for (state = 1; state <= Nstates; state++)
	{
	
	  // Calculate metrics if transition is allowed, otherwise skip to
	  // next state
	  if (Metrics_dist(state)<INF)
	    {
	      // Loop over all transitions
	      for (input = 1; input <= Ninputs; input++)
		{
		  
		  // Next state of the trellis 
		  next_state =(int) TRnextstate(state,input);
		  
		  // Modulated branch output (real part)
		  output_r = TRoutput_r(state,input);
		  // Modulated branch output (imaginary part)
		  output_i = TRoutput_i(state,input);

		  if (complex_ch)
		    {
		      // Modify branch output with channel gains (real)
		      faded_output_r = 
			fading_r(t)*output_r - fading_i(t)*output_i;
		  
		      // Modify branch output with channel gains (imag)
		      faded_output_i = 
			fading_i(t)*output_r + fading_r(t)*output_i;
		  
		      // Euclidian metric
		      euclid_metric = 
			((faded_output_r-input_r(t))
			 *(faded_output_r-input_r(t)))
			+((faded_output_i-input_i(t))
			  *(faded_output_i-input_i(t)));
		    }
		  else
		    {
		      // Modify branch output with channel gains (real)
		      faded_output_r = fading_r(t)*output_r;			
		      // Modify branch output with channel gains (imag)
		      faded_output_i = fading_r(t)*output_i;
			
		      // Euclidian metric
		      euclid_metric = 
			((faded_output_r-input_r(t))
			 *(faded_output_r-input_r(t)))
			+((faded_output_i-input_i(t))
			  *(faded_output_i-input_i(t)));		      
		    }
		  // Branch euclidian metric
		  Metrics_dist_tmp(next_state,state) = 
		    Metrics_dist(state) + euclid_metric;
		  
		  // Keep track of the previous paths that have lead to
		  // current states.
 		  TRstatepath_tmp(next_state,state,Range(1,t)) = 
 		    TRstatepath(state,Range(1,t));
		  
		  // Current state
		  TRstatepath_tmp(next_state,state,t+1) = next_state;
		  
		} // for (input = 1; input <= Ninputs; input++)
	      
	    } // if (Metrics_dist(state)<INF)
	  
	} // for (state = 1; state <= Nstates; state++)
      

      /* Find the survivors */ 
      
      // Initialise survivor pointer
      survivor_ptr = -1;

      // Find shortest metrices 
      Metrics_dist = min(Metrics_dist_tmp(i_ind, j_ind),j_ind);
      // Find indexes of the shortest metrices       
      ptr_tmp = minIndex(Metrics_dist_tmp(i_ind, j_ind),j_ind);

      // Have we gone past the trellis initialisation state?
      if (t>m+1)
 	survivor_ptr = ptr_tmp;
      else
	{
	  for (i = 1; i <= Nstates; i++)
	    {
	      if (any(Metrics_dist_tmp(i,Range(1,Nstates))<INF))
		survivor_ptr(i) = ptr_tmp(i);
	    }
 	}

      /* Keep track of the surviving paths */

      // Loop over next states
      for (i = 1; i <= Nstates; i++)
	{
	  // Has this state been reached yet?
	  if (survivor_ptr(i) > -1)
	    {
 	      TRstatepath(i,Range(1,t+1)) = 
 		TRstatepath_tmp(i,survivor_ptr(i),Range(1,t+1));
	    }
	} // for (i = 1; i <= Nstates; i++)
        

    } // for (t = 1; t <= Linput; t++)
  

  /* Final surviving state path */
  
  // Find the index of the minimum metric
  vec_tmp = minIndex(Metrics_dist);
  decode_ptr = vec_tmp(0);

  // Copy the final state path to output
  for (i = 1; i <= Linput+1; i++)
    *(output + (i-1)) =(double) TRstatepath(decode_ptr,i);
  
  return;

} // st_dec_CC



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

/*
 * conv_dec_soft(Data_received(C) (1xL), NextState, Output_mod(C),
 *               fading_gains(C) (1xL), m);
 *           
 */

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

  // Inputs (complex)
  double *input_r, *input_i, *fading_r, *fading_i;
  double *TRoutput_r, *TRoutput_i;
  // Inputs (real)
  double *TRnextstate;
  int m;
  // Outputs
  double *output;
  // Local variables
  int Sinput1, Sinput2, Nstates, Ninputs, Sfading1;
  bool complex_ch, complex_mod, complex_in;
  
  /* Check for proper number of arguments */  
  if (nrhs != 5)
    mexErrMsgTxt("'conv_dec_soft' requires five input arguments.");
  

  /**************************************
   * Assign pointers to input variables *
   * ************************************/
  
  /* Received / estimated symbols */
  complex_in = mxIsComplex(prhs[0]);
  input_r = mxGetPr(prhs[0]);		/* real part of input  */
  if (complex_in)  
    input_i = mxGetPi(prhs[0]);		/* imaginary part of input */
  else
    mexErrMsgTxt("'conv_dec_soft' works only with complex constellations.");
  //    input_i = mxGetPr(prhs[0]);		/* DUMMY */

  /* Next states of the trellis transitions */
  TRnextstate = mxGetPr(prhs[1]);

  /* Modulated symbol outputs corresponding to trellis transitions */
  TRoutput_r = mxGetPr(prhs[2]);	 /* real part of trellis output */
  TRoutput_i = mxGetPi(prhs[2]);	 /* imaginary part of trellis output */

  /* Channel fading information */
  complex_ch = mxIsComplex(prhs[3]);
  fading_r = mxGetPr(prhs[3]);		/* real part of fading coefficients */
  if (complex_ch)
    fading_i = mxGetPi(prhs[3]);	/* imaginary part of fading coeff's */
  else
    fading_i = mxGetPr(prhs[3]); 	/* DUMMY */
  
  m =(int) mxGetScalar(prhs[4]);	/* length of encoder memory */

  /* Size of the input matrix */
  Sinput1 = mxGetM(prhs[0]);
  Sinput2 = mxGetN(prhs[0]);
  
  /* Size of the look-up tables */
  Nstates = mxGetM(prhs[1]);
  Ninputs = mxGetN(prhs[1]);

  /* Length of the fading vector */
  Sfading1 = mxGetM(prhs[3]);

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

  /* Function call */
  conv_dec_CC(input_r, input_i, Sinput1, Sinput2, TRnextstate, Nstates, 
	      Ninputs, TRoutput_r, TRoutput_i, 
	      fading_r, fading_i, complex_ch, m, output);

  return;

}