tcc_binary_unix.c

This CD-ROM is distributed by Kluwer Academic Publishers with ABSOLUTELY NO SUPPORT and NO WARRANTY

  free(output1);
  free(rsc2_in_array);
  free(unuse);  
  free(output2);
}
 
void rsc_encd(int g[2][K], int input_len,int R, int *in_array, int *out_array_0, int *out_array)
{
  int m;                 /* K-1 */
  int i,j;               /* loop variables */
  int t,tt;              /* bit time, symble time */
  int *unencoded_data;   /* pointer to data array */
  int shift_reg[K];      /* the 1st RSC encoder shift register */
  int sr_head;           /* index to the first element in the shift register */
  int x,p;               /* feedback bit and output */
  int y;                 /* test bits */

  m = K-1;
  unencoded_data = malloc((input_len + m +1)*sizeof(int));
                         /* allocate space for the zero-padded input data array */

  for (t=0; t<input_len; t++)     /* read in the data and store it in the array */
    *(unencoded_data+t) = *(in_array+t);

  for (j=0; j<K; j++)
    shift_reg[j] = 0;    /* Initialize the shift register */

  sr_head = 0;           /* Initialize the index of the shift register */

  tt = 0;                /* Initialize the channel symbol output index */

  for (t=0; t<input_len+m; t++)          /* Start the encoding process */
    {
      shift_reg[sr_head] = *(unencoded_data+t);
      x = 0;
      p = 0;

      for (j=0; j<K; j++){
	i = (j+sr_head) % K;
	x ^=shift_reg[i] & g[0][j];
      }

      shift_reg[sr_head] = x;
      for (j=0; j<K; j++){
	i = (j+sr_head) % K;
	p ^=shift_reg[i] & g[1][j];
      }

      *(out_array+tt) = p;
      tt = tt + 1;
      
      sr_head -= 1;
      if ( sr_head<0 ) sr_head = m;   /* make sure we adjust pointer modulo K */
      
      if (R==1){                      /* get the 1st RSC encoded data */
        if (t>=input_len-1){          /* zero-pad the end of the data */     
	  shift_reg[sr_head] = 0;
	  y = 0;
          for (j=0; j<K; j++){
	  i = (j+sr_head) % K;
	  y ^=shift_reg[i] & g[0][j];
	  }
	  if (y==0) *(unencoded_data+t+1) = 0;
          else *(unencoded_data+t+1) = 1;
	}
      }
    }

  if (R==1){                          /* get the input data for 2nd RSC encoder */
  for (t=0; t<input_len+m; t++)
    *(out_array_0 + t) = *(unencoded_data + t);
  }

  free (unencoded_data);              /* free the dynamically allocated array */
}

/*******************************************/
/* BPSK Binary symmetric channel simulator */
/*******************************************/
float gngauss(float mean, float sigma);

void addnoise(float eb_ovr_n0, int channel_len, int *in_array, float *out_array)
{
  int t;
  float mean, eb, sn_ratio, sigma, signal,rate;

  mean = 0;
  eb = 1;
  rate = (float) 1.0/(2.0 + PUNCTURE);
  sn_ratio = (float) pow(10, ( eb_ovr_n0 / 10) );
  sigma =  (float) sqrt (eb / (2 * sn_ratio * rate) );

  for (t = 0; t < channel_len; t++) { /* transform the data from 0/1 to -1/+1 and add noise */
    signal = 2 * *( in_array + t ) - 1;
    *( out_array + t ) = signal + gngauss(mean, sigma);
  }
}

float gngauss(float mean, float sigma) 
{    
  double u, r;            /* uniform and Rayleigh random variables */

      /* generate a uniformly distributed random number u between 0 and 1 - 1E-6*/
  u = (double)rand() / RAND_MAX;
  if (u == 1.0) u = 0.999999999;

      /* generate a Rayleigh-distributed random number r using u */
  r = sigma * sqrt( 2.0 * log( 1.0 / (1.0 - u) ) );

      /* generate another uniformly-distributed random number u as before*/
  u = (double)rand() / RAND_MAX;
  if (u == 1.0) u = 0.999999999;

      /* generate and return a Gaussian-distributed random number using r and u */
  return( (float) ( mean + r * cos(2 * PI * u) ) );
}

/**************************************************/
/* Serial to parallel demultiplex at the receiver */
/**************************************************/
//void f_inter_block(int frame_size, int m, int n, float *read_in_array, float *read_out_array);

void de_multiplex(int alpha[FRAMESIZE], int puncture, float *in_array, 
                                   float *decd1_in_array, float *decd2_in_array)
{
  int i,j;                           /* loop variables */
  float x_sys[FRAMESIZE];            /* make a matrix to store the systematic data */
  float *info_array,*output_array;   /* pointer to systematic data array */
  int len_total;

  len_total = FRAMESIZE;

  info_array = malloc((len_total)*sizeof(int));
                                /* allocate space for systematic data array */
  output_array = malloc((len_total)*sizeof(int));
                                /* allocate space for interleaved data array */

  if (puncture==1){             /* unpunctured: rate = 1/3 */
    for (i=0; i<len_total; i++){
     x_sys[i] = *(in_array + 3*i);
     *(decd1_in_array + 2*i + 1) = *(in_array + 3*i + 1);
     *(decd2_in_array + 2*i + 1) = *(in_array + 3*i + 2);    
    }
  }
  else {                        /* punctured: rate = 1/2 */
    for (i=0; i<len_total; i++){
      x_sys[i] = *(in_array + 2*i);
      *(decd1_in_array + 2*i + 1) = 0;
      *(decd2_in_array + 2*i + 1) = 0;

      if (i%2==0) *(decd1_in_array + 2*i + 1) = *(in_array + 2*i + 1);
      else *(decd2_in_array + 2*i + 1) = *(in_array + 2*i + 1);
    }
  }

  for (i=0; i<len_total; i++){
    *(decd1_in_array + 2*i) = x_sys[i];      /* Extract the systematic bits for both decoders */
    *(info_array+i) = x_sys[i];
  }

  for (i=0; i<len_total; i++){
    *(output_array + i) = *(info_array + alpha[i]);
    *(decd2_in_array + 2*i) = *(output_array+i); /*Interleave the systematic bits for decoder 2*/
  }

  free(info_array);
  free(output_array);
}

/****************************************/
/* Turbo code decoding with log_map */
/****************************************/
void trellis(int g[2][K], int nextstate[NUMSTATES][2], int nextoutput[NUMSTATES][NUMSTATES],
		    int prestate[NUMSTATES][2], int preoutput[NUMSTATES][NUMSTATES]);
double max(double x, double y, double z, double w);
void logmap(int g[2][K], int len_total, int num_decd, float *rec_s, float *L_a, float *L_u)
{
  int i,j,l;                              /* loop variables */
  int k;                                  /* times */  
  int nextoutput[NUMSTATES][NUMSTATES];   /* maps current/nxt sts to next output */
  int nextstate[NUMSTATES][2];            /* for current state, gives next for given input */
  int preoutput[NUMSTATES][NUMSTATES];    /* maps previous/next states to previous output */
  int prestate[NUMSTATES][2];             /* gives conv. encoder output */

  double Alpha[len_total][NUMSTATES];
  double Beta[len_total][NUMSTATES];

  double gamma[NUMSTATES] = {-INF};  
  double tempmax[FRAMESIZE] = {0};
  double temp0[NUMSTATES] = {0};
  double temp1[NUMSTATES] = {0};
  double gamma0 = 0;
  double gamma1 = 0;

  for (i = 1; i < NUMSTATES; i++)      /* initialize Alpha data array */
    Alpha[0][i] = -INF;

  Alpha[0][0] = 0;

  for (i = 1; i < NUMSTATES; i++)      /* initialize Beta data array */
    Beta[len_total-1][i] = -INF;

  if (num_decd == 1) Beta[len_total-1][0] = 0;
  else if (num_decd == 2)
    {
      for (i = 0; i < NUMSTATES; i++)
	Beta[len_total-1][i] = 0;
    }
  else printf("\n decoder number is limited to 1 and 2!");

  trellis(g, nextstate, nextoutput, prestate, preoutput);

    #ifdef DEBUG
    printf("\nNext output:");
    for (j = 0; j < NUMSTATES; j++) {
        printf("\n");
        for (l = 0; l < NUMSTATES; l++)
            printf("%2d ", nextoutput[j][l]);
    } /* end j for-loop */

    printf("\nPrevious output:");
    for (j = 0; j < NUMSTATES; j++) {
        printf("\n");
        for (l = 0; l < NUMSTATES; l++)
            printf("%2d ", preoutput[j][l]);
    } /* end j for-loop */

    printf("\nNext State:");
    for (j = 0; j < NUMSTATES; j++) {
        printf("\n");
        for (l = 0; l < 2; l++)
            printf("%2d ", nextstate[j][l]);
    } /* end j for-loop */

    printf("\nPrevious State:");
    for (j = 0; j < NUMSTATES; j++) {
        printf("\n");
        for (l = 0; l < 2; l++)
            printf("%2d", prestate[j][l]);
    } /* end j for-loop */
    #endif

  for (k=1; k<len_total; k++) {        /* trace forward, compute Alpha */
    for (i=0; i<NUMSTATES; i++) {
      for (j=0; j<NUMSTATES; j++)
	gamma[j] = -INF;                 /* initialize gamma */

      gamma[prestate[i][0]] = -*(rec_s+2*k-2) + *(rec_s+2*k-1)*preoutput[i][prestate[i][0]] - log(1+exp(*(L_a+k-1)));
      gamma[prestate[i][1]] = *(rec_s+2*k-2) + *(rec_s+2*k-1)*preoutput[i][prestate[i][1]] + *(L_a+k-1) - log(1+exp(*(L_a+k-1)));

      if ( ( exp(gamma[prestate[i][0]]+Alpha[k-1][prestate[i][0]])
            +exp(gamma[prestate[i][1]]+Alpha[k-1][prestate[i][1]]) )<1e-300)
	Alpha[k][i] = -INF;
      else Alpha[k][i] = log(exp(gamma[prestate[i][0]] + Alpha[k-1][prestate[i][0]])
                            +exp(gamma[prestate[i][1]] + Alpha[k-1][prestate[i][1]]));
    }				

    tempmax[k] = max(Alpha[k][0], Alpha[k][1],Alpha[k][2], Alpha[k][3]);
 
    for (j=0; j<NUMSTATES; j++){    /* to be sure that Alpha do not overflow */ 

      Alpha[k][j] = Alpha[k][j] - tempmax[k];
    }
  }
   
  for (k=len_total-2; k>=0; k--) {        /* trace forward, compute Beta */
    for (i=0; i<NUMSTATES; i++) {
      for (j=0; j<NUMSTATES; j++)
	gamma[j] = -INF;                 /* initialize Beta */

      gamma[nextstate[i][0]] = -*(rec_s+2*k+2) + *(rec_s+2*k+3)*nextoutput[i][nextstate[i][0]] - log(1+exp(*(L_a+k+1)));
      gamma[nextstate[i][1]] = *(rec_s+2*k+2) + *(rec_s+2*k+3)*nextoutput[i][nextstate[i][1]]+ *(L_a+k+1)-log(1+exp(*(L_a+k+1)));

      if ((exp(gamma[nextstate[i][0]]+Beta[k+1][nextstate[i][0]])+exp(gamma[nextstate[i][1]]+Beta[k+1][nextstate[i][1]]))<1e-300)
	Beta[k][i] = -INF;
      else Beta[k][i] =log(exp(gamma[nextstate[i][0]] + Beta[k+1][nextstate[i][0]])
                          +exp(gamma[nextstate[i][1]] + Beta[k+1][nextstate[i][1]]));
    }

    for (j=0; j<NUMSTATES; j++)    /* to be sure that Beta do not overflow */ 
      Beta[k][j] = Beta[k][j] - tempmax[k+1];
  }

  for (k=0; k<len_total; k++) {        /* compute the soft output, log_likelihood ratio */
    for (i=0; i<NUMSTATES; i++) {
      gamma0 = -*(rec_s+2*k) + *(rec_s+2*k+1)*preoutput[i][prestate[i][0]] - log(1 + exp(*(L_a+k)));
      gamma1 = *(rec_s+2*k) + *(rec_s+2*k+1)*preoutput[i][prestate[i][1]] + *(L_a+k) - log(1 + exp(*(L_a+k)));

      temp0[i] = exp(gamma0 + Alpha[k][prestate[i][0]] + Beta[k][i]);
      temp1[i] = exp(gamma1 + Alpha[k][prestate[i][1]] + Beta[k][i]);
    }
    *(L_u+k) = log(temp1[0]+temp1[1]+temp1[2]+temp1[3]) - log(temp0[0]+temp0[1]+temp0[2]+temp0[3]);
  }
}

double max(double x, double y, double z, double w)
{
  double maxmum;

  maxmum = x>y?x:y;        
  maxmum = maxmum>z?maxmum:z;
  maxmum = maxmum>w?maxmum:w;
  
  return(maxmum);
}
       
/*************************************************/
/* Set up the trellis for a given code generator */
/*************************************************/
int bin2deci(int *b, int size);
void deci2bin(int d, int size, int *b);
int nxt_stat(int g[2][K], int current_state, int input, int *memory_contents);

void trellis(int g[2][K], int nextstate[NUMSTATES][2], int nextoutput[NUMSTATES][NUMSTATES],
	     int prestate[NUMSTATES][2], int preoutput[NUMSTATES][NUMSTATES])
{
  int i,j,l;                            /* loop variables */
  int memory_contents[K];               /* feedback bit + conv. encoder sr */
  int encd_output;                      /* store trial encoder output */
  int n, next_state;                    

  n = 2;                                /* n is 2^1 = 2 for rate 1/2 */ 
  for (i = 0; i < NUMSTATES; i++) {     /* initialize data structures */
    for (j = 0; j < NUMSTATES; j++){
      preoutput[i][j] = 0;
      nextoutput[i][j] = 0;
    }

    for (j = 0; j < n; j++){
      prestate[i][j] = 0; 
      nextstate[i][j] = 0;
    }
  }

  for (j = 0; j < NUMSTATES; j++) {     /* generate the matrices */
    for (l = 0; l < n; l++) {
      next_state = nxt_stat(g,j, l, memory_contents);

      encd_output = 0;
      for (i = 0; i < K; i++) 
        encd_output ^= memory_contents[i] & g[1][i];

      nextstate[j][l] = next_state;              /* next state, given current state and input */
      nextoutput[j][next_state] = 2*encd_output - 1;   /* generate the next output */

      prestate[next_state][l] = j;
      preoutput[next_state][j] = nextoutput[j][next_state];
    }  /* end of l for loop */
  }  /* end of j for loop */
}
/* ************************************************************************** */
int nxt_stat(int g[2][K], int current_state, int input, int *memory_contents) {

  int binary_state[K - 1];              /* binary value of current state */
  int next_state_binary[K - 1];         /* binary value of next state */
  int next_state;                       /* decimal value of next state */
  int i;                                /* loop variable */

  deci2bin(current_state, K - 1, binary_state);
                     /* convert the decimal value of the current state number to binary */

  for (i=0; i<K-1; i++){ 
    input ^= binary_state[i] & g[0][i+1];      /* get the feedback bit */
    next_state_binary[i] = 0;
  }

  next_state_binary[0] = input;
                     /* given the input and current state number, compute the next state number */
  for (i = 1; i < K - 1; i++)
    next_state_binary[i] = binary_state[i - 1];

  next_state = bin2deci(next_state_binary, K - 1);
                     /* convert the binary value of the next state number to decimal */

  memory_contents[0] = input;   /* memory_contents are feedback bits in the encoder */
  for (i = 1; i < K; i++)
    memory_contents[i] = binary_state[i - 1];
  return(next_state);
}

void deci2bin(int d, int size, int *b) {  /* converts a decimal number to a binary number */
   int i;

   for(i = 0; i < size; i++)
     b[i] = 0;

   b[size - 1] = d & 0x01;

   for (i = size - 2; i >= 0; i--) {
     d = d >> 1;
     b[i] = d & 0x01;
   }
}

int bin2deci(int *b, int size) {          /* converts a binary number to decimal number */
   int i, d;

   d = 0;

   for (i = 0; i < size; i++)
     d = d + (int)pow(2, size - i - 1) * b[i];

   return(d);
}