tcc_dvb-rcs_pc.c

MAP 算法 源程序 第三代 数字通信 turbo码相关

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