st_dec_cc.cpp

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

/*
 *-----------------------------------------------------------------------------
 * This is a C++ mex-file for MATLAB
 *
 * Decode received 1-D data vector with trellis (soft inputs).
 *
 *
 * Format:
 * -------
 *
 * outputs = ...
 *   st_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 st_dec_CC(double *in_r, double *in_i, int Sinput1, int Sinput2, 
		      double *nstate, int Nstates, int Ninputs, 
		      double *TRout_1_r, double *TRout_1_i, 
		      double *TRout_2_r, double *TRout_2_i, 
		      double *fad_r, double *fad_i, int m, int SF, int NRx,
		      double *output)

{

  /* Local variables */ 

  int i, j, t, state, input, Linput, decode_ptr, next_state, rx_ant, c;
  double faded_output_r=0.0, faded_output_i=0.0;
  double tmp_dist;
  const double INF = 1E+100;

  // 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;
  DArray1D fad_tmp_r(Range(1,2));
  DArray1D fad_tmp_i(Range(1,2));
  IArray1D survivor_ptr(Range(1,Nstates));
  IArray1D ptr_tmp(Range(1,Nstates));
  DArray1D output_r(Range(1,2));
  DArray1D output_i(Range(1,2));
  DArray1D euclid_metric(Range(1,NRx));
  DArray1D Metrics_dist(Range(1,Nstates));
  IArray2D TRstatepath(Range(1,Nstates), Range(1,Sinput2/SF+1));
  DArray2D Metrics_dist_tmp(Range(1,Nstates), Range(1,Nstates));
  IArray3D TRstatepath_tmp(Range(1,Nstates), Range(1,Nstates), 
                           Range(1,Sinput2/SF+1));

  // input matrices 
  DArray2D input_r(in_r, shape(Sinput1,Sinput2) , neverDeleteData, 
		   fortranArray);
  DArray2D input_i(in_i, shape(Sinput1,Sinput2), neverDeleteData, 
		   fortranArray);
  DArray2D fading_r(fad_r, shape(2*NRx,Sinput2), neverDeleteData, 
		    fortranArray);
  DArray2D fading_i(fad_i, shape(2*NRx,Sinput2), neverDeleteData, 
		    fortranArray);
  DArray2D TRoutput_1_r(TRout_1_r, shape(Nstates,Ninputs), neverDeleteData, 
			fortranArray);
  DArray2D TRoutput_1_i(TRout_1_i, shape(Nstates,Ninputs), neverDeleteData, 
		      fortranArray);
  DArray2D TRoutput_2_r(TRout_2_r, shape(Nstates,Ninputs), neverDeleteData, 
		      fortranArray);
  DArray2D TRoutput_2_i(TRout_2_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 decisions) --
   *
   */

  Linput = Sinput2/SF;		// 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);
		  
		  /* TX Antenna 1 */
		  // Modulated branch output (real part)
		  output_r(1) = TRoutput_1_r(state,input);
		  // Modulated branch output (imaginary part)
		  output_i(1) = TRoutput_1_i(state,input);

		  /* TX Antenna 2 */
		  // Modulated branch output (real part)
		  output_r(2) = TRoutput_2_r(state,input);
		  // Modulated branch output (imaginary part)
		  output_i(2) = TRoutput_2_i(state,input);

          euclid_metric = 0.0;

		  /* Loop over receive antennas */
		  for (rx_ant = 1; rx_ant <= NRx; rx_ant++)
		    {
              for (c = 1; c <= SF; c++)
                {
                  /* THIS IS KLUDGE!!! FIND BETTER WAY TO DO THIS!!! */
                  if (rx_ant == 1)
                    {
                      // Modify branch output with channel gains (real)
                      fad_tmp_r(Range(1,2)) = 
                        fading_r(Range(1,2),SF*(t-1)+c)*output_r(Range(1,2))-
                        fading_i((Range(1,2)),SF*(t-1)+c)*output_i(Range(1,2));
                      
                      // Modify branch output with channel gains (imag)
                      fad_tmp_i(Range(1,2)) = 
                        fading_i((Range(1,2)),SF*(t-1)+c)*output_r(Range(1,2))+
                        fading_r((Range(1,2)),SF*(t-1)+c)*output_i(Range(1,2));
                    }
                  else
                    {
                      // Modify branch output with channel gains (real)
                      fad_tmp_r(Range(1,2)) = 
                        fading_r(Range(3,4),SF*(t-1)+c)*output_r(Range(1,2)) -
                        fading_i((Range(3,4)),SF*(t-1)+c)*output_i(Range(1,2));
                      
                      // Modify branch output with channel gains (imag.)
                      fad_tmp_i(Range(1,2)) = 
                        fading_i((Range(3,4)),SF*(t-1)+c)*output_r(Range(1,2))+
                        fading_r((Range(3,4)),SF*(t-1)+c)*output_i(Range(1,2));
                    } // (rx_ant == 1)
                  
                  faded_output_r = sum(fad_tmp_r);
                  faded_output_i = sum(fad_tmp_i);

                  // Euclidian metric
                  euclid_metric(rx_ant) = euclid_metric(rx_ant)+
                    (faded_output_r-input_r(rx_ant,SF*(t-1)+c))*
                    (faded_output_r-input_r(rx_ant,SF*(t-1)+c))+
                    (faded_output_i-input_i(rx_ant,SF*(t-1)+c))*
                    (faded_output_i-input_i(rx_ant,SF*(t-1)+c));

                } // for (c = 1; c <= SF; c++)
		    } // for (rx_ant = 1; rx_ant <= NRx; rx_ant++)
		 
		  // Branch euclidian metric
		  Metrics_dist_tmp(next_state,state) = 
		    Metrics_dist(state) + sum(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 */
/***********/

/*
 * st_dec_CC(Data_received(C) (2xL), NextState, Output_mod_1(C),
 *           Output_mod_2(C), fading_gains(C) (4xL), m, SF);
 *           
 */

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

  /* Inputs (complex) */
  double *input_r, *input_i, *fading_r, *fading_i;
  double *TRoutput_1_r, *TRoutput_1_i,*TRoutput_2_r, *TRoutput_2_i;
  /* Inputs (real) */
  double *TRnextstate;
  int m, SF;
  /* Outputs */
  double *output;
  /* Local variables */
  int Sinput1, Sinput2, Nstates, Ninputs, Sfading1, NRx;
  
  /* Check for proper number of arguments */  
  if (nrhs != 7)
    mexErrMsgTxt("'st_dec_CC' requires seven input arguments.");
  
  /* Assign pointers to input variables */
  input_r = mxGetPr(prhs[0]);		 /* real part of input  */
  input_i = mxGetPi(prhs[0]);		 /* imaginary part of input */
  TRnextstate = mxGetPr(prhs[1]);

  /* Antenna 1 */
  TRoutput_1_r = mxGetPr(prhs[2]);	 /* real part of trellis output */
  TRoutput_1_i = mxGetPi(prhs[2]);	 /* imaginary part of trellis output */
  /* Antenna 2 */
  TRoutput_2_r = mxGetPr(prhs[3]);	 /* real part of trellis output */
  TRoutput_2_i = mxGetPi(prhs[3]);	 /* imaginary part of trellis output */

  fading_r = mxGetPr(prhs[4]);		 /* real part of fading coefficients */
  fading_i = mxGetPi(prhs[4]);		 /* imaginary part of fading coeff's */

  m =(int) mxGetScalar(prhs[5]);	 /* encoder memory */
  SF =(int) mxGetScalar(prhs[6]);	 /* spreading factor */

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

  /* Size of the fading vector and number of Rx antennas */
  Sfading1 = mxGetM(prhs[4]);
  NRx = Sfading1/2;

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

  /* Function call */
  st_dec_CC(input_r, input_i, Sinput1, Sinput2, TRnextstate, Nstates, 
            Ninputs, TRoutput_1_r, TRoutput_1_i, TRoutput_2_r, TRoutput_2_i, 
            fading_r, fading_i, m, SF, NRx, output);

  return;

}