/* Non-binary Turbo code Test Driver: Max_log_map *//* Creation date : Jan.17 2001   *//* Programmer: Gao yingzi        *//* Last modified date: Nov.16 2001 *//* This program simulates the DVB-RCS encoding-decoding system on PC. *//* It uses the CRSC codes in DVB-RCS system model described in Figure 3.12 in Chapter 3. *//* Two component CRSC ( Circular Recursive Systematic Convolutional) encoders are used.  *//* Random information bits are modulated into +1/-1, and transmitted through an AWGN channel. *//* Max-Log-MAP algorithm without quantization or approximation is used.      *//* By making use of ln(e^x+e^y+e^z+e^w) = max (x,y,z,w), we have MAX-Log-MAP.*//* To set the number of iterations, the globle variable "DECITER" should be changed. *//* To set the frame length, the globle variable "N" should be changed.    *//* To set the code rate, the globle variable "R" should be changed.       *//* To simulate a different range of Eb/N0, change the globle variable LOEBNO and HIEBNO, *//* as well as the step size of Eb/No increment. */#include <malloc.h>#include <math.h>#include <stdio.h>#include <stdlib.h>#include <time.h>#define K 4              /* constraint length */#define NUMSTATES 8      /* 2^(K - 1) -- change as required */ #define PI 3.141592654   /* circumference of circle divided by diameter */#define N 212            /* Block sizes (53 bytes) */#define R 3              /* code rate number: R=1, code rate r=1/3; R=2, code rate r=2/5;     */                         /* R=3, code rate r=1/3; R=4, code rate r=2/3; R=5, code rate r=3/4; */                         /* R=6, code rate r=4/5; R=7, code rate r=6/7; */#define DECITER 8        /* the number of decoding iterations */#define LOEBN0 1.5       /* minimum Eb/No at which to test */#define HIEBN0 3.5       /* maximum Eb/No at which to test */#define EBN0STEP 0.5     /* Eb/No increment for test driver */void gen01dat(int data_len, int *out_array);void permutation(int *alpha);void deci2bin(int d, int size, int *b);void turbo_encd(int g[3][K], int M, int *alpha, int *in_array, int *turbo_en_out);void addnoise(float eb_ovr_n0, int encd_block_len, float rate, int *in_array, float *out_array);void de_multiplex(int msg_length, int M, int *alpha, float *in_array, float *sys_array,		  float *pinfo_array, float *parity1_array, float *parity2_array);void maxlogmap(int g[3][K], float *rec_s, float *rec_p, float *L_a, float *L_u, float *Alpha0, float *BetaN);float max3(float x, float y, float z);main(){  int i, j, ll, iter, t, start;   /* loop variables */  int encd_block_len,total_len;   /* the length of the encoded block depending on the code rate */  int M,msg_length,channel_len;   /* parameter of encoded block length */  int a,dec_iter;                 /* turbo decoding iteration */  int binary_input[2];  int *alpha;                     /* permutation table */  int *onezer;  int *encoded;                   /* original, encoded array */  float *splusn;                  /* noisy data array */  float *sys_array;               /* received systematic array for decoder one */  float *pinfo_array;             /* received systematic array, interleaved for decoder two */  float *parity1_array;           /* received parity1 array for decoder one */  float *parity2_array;           /* received parity2 array for decoder two */  float *L_a;                     /* deinterleaved or interleaved extrinsic info */  float *L_e_1;                   /* decoder two output, then deinterleaved as input of decoder one */  float *L_e_2;                   /* decoder one output, as decoder two input */  float *L_u;                     /* log_likelihood ratios */  float *L_all;                   /* finial soft output after iterations */  int *turbout;                   /* estimated bit of decoder */  float yx[N][K],L_a_p,L_e_2_p;   /* Medium values for caculating extrinsic info, interleaving, deinterleaving */  float eb_ovr_n0,rate,L_c,bit_errors,temp_bit_errors,frame_errors,e_threshold,e_ber,f_ber;   /* various statistics */  float Ad1[NUMSTATES],Bd1[NUMSTATES],Ad2[NUMSTATES],Bd2[NUMSTATES];   /* circular trellis state values */  float *Alpha0,*BetaN;  int g[3][K] = {{1, 1,0,1},       /* 15 */	         {1, 0,1,1},       /* 13 */		 {1, 0,0,1}};      /* 11 */  FILE *fp;  fp = fopen("TCC_DVB-RCS_PC_data","a");  if ((fp = fopen("TCC_DVB-RCS_PC_data","a"))==NULL){    printf("cannot open file\n");    exit (0);  }  printf("\ng1 = %d%d%d%d", g[0][0], g[0][1], g[0][2], g[0][3]);  printf("\ng2 = %d%d%d%d", g[1][0], g[1][1], g[1][2], g[1][3]);  printf("\ng3 = %d%d%d%d\n", g[2][0], g[2][1], g[2][2], g[2][3]);  fprintf(fp,"\ng1 = %d%d%d%d", g[0][0], g[0][1], g[0][2], g[0][3]);  fprintf(fp,"\ng2 = %d%d%d%d", g[1][0], g[1][1], g[1][2], g[1][3]);  fprintf(fp,"\ng3 = %d%d%d%d\n", g[2][0], g[2][1], g[2][2], g[2][3]);  #if R == 1  rate = 1.0/3.0;  M = N;  encd_block_len = 2*N + M;  #endif  #if R == 2  rate = 2.0/5.0;  M = N/2;  encd_block_len = 2*N + M;  #endif  #if R == 3  rate = 1.0/2.0;  M = N;  encd_block_len = N + M;  #endif  #if R == 4  rate = 2.0/3.0;  M = N/2;  encd_block_len = N + M;  #endif  #if R == 5  rate = 3.0/4.0;  if (N%3 == 0) M = N/3;  else if (N%3 == 1) M = (N-4)/3 + 2;  else if (N%3 == 2) M = (N-8)/3 + 3;  encd_block_len = N + M;  #endif  #if R == 6  rate = 4.0/5.0;  M = N/4;  encd_block_len = N + M;  #endif  #if R == 7  rate = 6.0/7.0;  if (N%3 == 0) M = N/6;  else if (N%3 == 1) M = (N-4)/6 + 1;  else if (N%3 == 2) M = (N-8)/6 + 2;  encd_block_len = N + M;  #endif  msg_length = 2*N;  total_len = 3*N;  channel_len = 2*encd_block_len;  printf("\nK = %d   Frame size N = %d  M = %d  code rate = %f\n", K, N, M, rate);  fprintf(fp,"\nK = %d   Frame size N = %d  M = %d  code rate = %f\n", K, N, M, rate);  fclose(fp);  for (eb_ovr_n0 = LOEBN0; eb_ovr_n0 <= HIEBN0; eb_ovr_n0 += EBN0STEP) {    fp = fopen("TCC_DVB-RCS_PC_data","a");    if ((fp = fopen("TCC_DVB-RCS_PC_data","a"))==NULL){      printf("cannot open file\n");      exit (0);    }    start = time(NULL);    a = 1;    L_c = 4 * a * rate * pow(10,eb_ovr_n0/10);     bit_errors = 0.0;    frame_errors = 0.0;    e_ber = 0.0;    f_ber = 0.0;    iter = 0;    if (eb_ovr_n0 <= 6.5) e_threshold = 100; /* +/- 20% */    else e_threshold = 20;                   /* +/- 100% */    printf("\nreliability = %f   BER threshold = %f", L_c, e_threshold);    fprintf(fp,"\nReliability = %f   BER threshold = %f", L_c, e_threshold);    while (bit_errors < e_threshold) {      alpha = malloc( N * sizeof(int) );      if (alpha == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating interleaver table, aborting!");        exit(1);      }       onezer = malloc( msg_length * sizeof(int) );      if (onezer == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating onezer array, aborting!");        exit(1);      }      encoded = malloc( channel_len * sizeof(int) );      if (encoded == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating encoded array, aborting!");        exit(1);      }      splusn = malloc( channel_len * sizeof(float) );      if (splusn == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating splusn array, aborting!");        exit(1);      }      sys_array = malloc( msg_length * sizeof(float) );      if (sys_array == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating received symbol sys array, aborting!");        exit(1);      }      pinfo_array = malloc( msg_length * sizeof(float) );      if (pinfo_array == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating received symbol interleaved sys array, aborting!");        exit(1);      }      parity1_array = malloc( msg_length * sizeof(float) );      if (parity1_array == NULL) {        printf("\n TCC_DVB-RCS_UNIX.c:  error allocating received symbol parity1 array, aborting!");        exit(1);      }      parity2_array = malloc( msg_length * sizeof(float) );      if (parity2_array == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating received symbol parity2 array, aborting!");        exit(1);      }      Alpha0 = malloc( NUMSTATES * sizeof(float) );      if (Alpha0 == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating the values of circular trellis state, aborting!");        exit(1);      }      BetaN = malloc( NUMSTATES * sizeof(float) );      if (BetaN == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating the values of circular trellis state, aborting!");        exit(1);      }      L_a = malloc( total_len * sizeof(float) );      if (L_a == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating L_a array, aborting!");        exit(1);      }      L_e_1 = malloc( total_len * sizeof(float) );      if (L_e_1 == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating L_e_1 array, aborting!");        exit(1);      }      L_e_2 = malloc( total_len * sizeof(float) );      if (L_e_2 == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating L_e_2 array, aborting!");        exit(1);      }      L_u = malloc( total_len * sizeof(float) );      if (L_u == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating L_u array, aborting!");        exit(1);      }      L_all = malloc( total_len * sizeof(float) );      if (L_all == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating L_all array, aborting!");        exit(1);      }      turbout = malloc( msg_length * sizeof(int) );      if (turbout == NULL) {        printf("\n TCC_DVB-RCS_PC.c:  error allocating turbo decoding out array, aborting!");        exit(1);      }      iter += 1;      gen01dat(msg_length, onezer);      permutation(alpha);      turbo_encd(g, M, alpha, onezer, encoded);      addnoise(eb_ovr_n0, encd_block_len, rate, encoded, splusn);      de_multiplex(msg_length, M, alpha, splusn, sys_array, pinfo_array, parity1_array, parity2_array);      for (j=0; j<msg_length; j++){ 	*(sys_array+j) = 0.5*L_c**(sys_array+j);	*(pinfo_array+j) = 0.5*L_c**(pinfo_array+j);	*(parity1_array+j) = 0.5*L_c**(parity1_array+j);	*(parity2_array+j) = 0.5*L_c**(parity2_array+j);      }      for (j=0; j<total_len; j++){           /* initial extrinsic info and log_likelihood ratio */	*(L_a+j) = 0.0;	*(L_e_1+j) = 0.0;	*(L_e_2+j) = 0.0;	*(L_u+j) = 0.0;	*(L_all+j) = 0.0;      }      for (i=0; i<NUMSTATES; i++){           /* initial circular trellis values */	Ad1[i] = 0.0;	Bd1[i] = 0.0;	Ad2[i] = 0.0;	Bd2[i] = 0.0;      }      for (dec_iter=0; dec_iter<DECITER; dec_iter++){	for (i=0; i<3; i++){                 /* deinterleave the extrinsic info for Decoder one */          for (j=0; j<N; j++)                /* deinterleave the level two */	    *(L_a + *(alpha+j) + i*N) = *(L_e_1 + j + i*N);	}        for (j=0; j<N; j++){                 /* deinterleave the level one */          if (j%2==0){	    L_a_p = *(L_a+j);	    *(L_a+j) = *(L_a+j+N);	    *(L_a+j+N) = L_a_p;	  }	}        for (i=0; i<NUMSTATES; i++){	  *(Alpha0 + i) = Ad1[i];	  *(BetaN + i) = Bd1[i];	}        maxlogmap(g, sys_array, parity1_array, L_a, L_u, Alpha0, BetaN);   /* Decoder one */        for (i=0; i<NUMSTATES; i++){	  Ad1[i] = *(Alpha0 + i);	  Bd1[i] = *(BetaN + i);	}        for (j=0; j<N; j++){                 /* get extrinsic info for Decoder two */          for (i=0; i<4; i++){	    deci2bin(i, 2, binary_input);	    for (ll=0; ll<2; ll++)	      binary_input[ll] = 1 - 2*binary_input[ll];	    yx[j][i] = *(sys_array + 2*j)*binary_input[1] + *(sys_array + 2*j + 1)*binary_input[0];	  }          for (i=0; i<3; i++)            *(L_e_2+j+i*N) = *(L_u+j+i*N) - yx[j][i+1] + yx[j][0] - *(L_a+j+i*N);	}        for (j=0; j<N; j++){                 /* Interleave the extrinsic info for Decoder two, level one */          if (j%2==0){	    L_e_2_p = *(L_e_2+j);	    *(L_e_2+j) = *(L_e_2+j+N);	    *(L_e_2+j+N) = L_e_2_p;	  }	}	for (i=0; i<3; i++){                 /* level two */          for (j=0; j<N; j++)	    *(L_a + j + i*N) = *(L_e_2 + i*N + *(alpha+j));	}        for (i=0; i<NUMSTATES; i++){	  *(Alpha0 + i) = Ad2[i];	  *(BetaN + i) = Bd2[i];	}	maxlogmap(g, pinfo_array, parity2_array, L_a, L_u, Alpha0, BetaN);  /* Decoder two */        for (i=0; i<NUMSTATES; i++){	  Ad2[i] = *(Alpha0 + i);	  Bd2[i] = *(BetaN + i);	}        for (j=0; j<N; j++){                 /* get extrinsic info for Decoder one */          for (i=0; i<4; i++){	    deci2bin(i, 2, binary_input);	    for (ll=0; ll<2; ll++)	      binary_input[ll] = 1 - 2*binary_input[ll];	    yx[j][i] = *(pinfo_array + 2*j)*binary_input[1] + *(pinfo_array + 2*j + 1)*binary_input[0];	  }          for (i=0; i<3; i++)            *(L_e_1+j+i*N) = *(L_u+j+i*N) - yx[j][i+1] + yx[j][0] - *(L_a+j+i*N);	}      }  /* end dec_iter-loop */      for (i=0; i<3; i++){                   /* deinterleave the log_likelihood ratios */        for (j=0; j<N; j++)	  *(L_all + *(alpha+j) + i*N) = *(L_u + j + i*N);      }      for (j=0; j<N; j++){                   /* estimate the decoded output */	if ( max3(*(L_all+j),*(L_all+j+N),*(L_all+j+2*N)) > 0.0 ){	  if ( *(L_all+j)==max3(*(L_all+j),*(L_all+j+N),*(L_all+j+2*N)) ){	    if(j%2 == 0){                    /* deinterleaver for j mod 2 = 0, invert A and B */	      *(turbout+2*j) = 0;	      *(turbout+2*j+1) = 1;	    }	    else {	      *(turbout+2*j) = 1;	      *(turbout+2*j+1) = 0;	    }	  }	  else if ( *(L_all+j+N)==max3(*(L_all+j),*(L_all+j+N),*(L_all+j+2*N)) ){	    if(j%2 == 0){                    /* deinterleaver for j mod 2 = 0, invert A and B */	      *(turbout+2*j) = 1;	      *(turbout+2*j+1) = 0;	    }	    else {	      *(turbout+2*j) = 0;	      *(turbout+2*j+1) = 1;	    }	  }	  else if ( *(L_all+j+2*N)==max3(*(L_all+j),*(L_all+j+N),*(L_all+j+2*N)) ){	    *(turbout+2*j) = 1;	    *(turbout+2*j+1) = 1;	  }	}	else {	  *(turbout+2*j) = 0;	  *(turbout+2*j+1) = 0;	}      } /* end j for-estimating loop */      temp_bit_errors = bit_errors;      for (t = 0; t < msg_length; t++){     /* count the bit errors */        if (*(onezer + t) != *(turbout + t))          bit_errors = bit_errors + 1;      } /* end t count bit errors-loop */      if ((bit_errors - temp_bit_errors) > 0) frame_errors = frame_errors + 1;   /* count the frame errors */          free(alpha);      free(onezer);      free(encoded);      free(splusn);      free(sys_array);      free(pinfo_array);      free(parity1_array);      free(parity2_array);      free(Alpha0);      free(BetaN);      free(L_a);      free(L_e_1);      free(L_e_2);      free(L_u);      free(L_all);      free(turbout);    } /* end of while loop */    e_ber = bit_errors / (msg_length * iter);    f_ber = frame_errors / iter;    printf("\nThe elapsed time was %d seconds for %d frames and %d decoding iterations",                   time(NULL) - start, iter, dec_iter);    fprintf(fp,"\nThe elapsed time was %d seconds for %d frames and %d decoding iterations",                   time(NULL) - start, iter, dec_iter);    printf("\nAt %1.1fdB Eb/No, ", eb_ovr_n0);    fprintf(fp,"\nAt %1.1fdB Eb/No, ", eb_ovr_n0);