/*
 * File: tc-8psk-iq.c
 * Date: January 8, 1998
 * Description: Two-stage decoding of TC-8PSK
 *
 * The I and Q components of the received 8-PSK symbols are mapped into
 * a two-dimensional space with a QPSK constellation whose points represent
 * cosets of BPSK sets in the original 8-PSK constellation.
 * 
 */
// ------------------------------------------------------------------------
// This program is complementary material for the book:
//
// R.H. Morelos-Zaragoza, The Art of Error Correcting Coding, Wiley, 2002.
//
// ISBN 0471 49581 6
//
// This and other programs are available at http://the-art-of-ecc.com
//
// You may use this program for academic and personal purposes only. 
// If this program is used to perform simulations whose results are 
// published in a journal or book, please refer to the book above.
//
// The use of this program in a commercial product requires explicit
// written permission from the author. The author is not responsible or 
// liable for damage or loss that may be caused by the use of this program. 
//
// Copyright (c) 2002. Robert H. Morelos-Zaragoza. All rights reserved.
// ------------------------------------------------------------------------
 
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <limits.h>
 

int k2=1, n2, m2;
int NUM_STATES, OUT_SYM, NUM_TRANS;
long TRUNC_LENGTH;

double RATE;
double INIT_SNR, FINAL_SNR, SNR_INC;
long NUMSIM;

char name1[40], name2[40];

#define MAX_RANDOM LONG_MAX 

//#define DELTA 1000
//#define SHOW_PROGRESS       /* Print simulation progress, */

FILE  *fp;                 /* Pointer for trellis data file */

//unsigned int g2[n2][k2] =  { /* rate-1/2 memory=6 */
//        0x4f, /* REVERSE ORDER */
//        0x6d }; /* */

unsigned int g2[10][10];
unsigned int memory2, output;            /* Memory,output transmit encoder */
unsigned int memory, output2;            /* Memory,output receive encoder */
unsigned int data2;                      /* Data */

unsigned long seed;                      /* Seed for random generator */
unsigned int data_symbol[1024];          /* 1-bit data sequence */
unsigned int data_symbol2[1024];         /* 1-bit data sequence */
unsigned long indxx;                     /* Simulation index */


double psk_I[8], psk_Q[8];        /* Coordinated of 8-PSK signals */
double qpsk_I[4], qpsk_Q[4];      /* Coordinated of 4-PSK signals */

int s[4];                         /* Coset-sector look-up table */

int transmitted;                  /* index of transmitted signal */
double transmitted_I;             /* Transmitted signals/branch */
double transmitted_Q;             /* Transmitted signals/branch */
int reconstructed;                /* Decoded coset index */
int estimate_data2;               /* Estimate of uncoded info. bit */
double snr, amp;

double received_I;                       /* Received signals/branch */
double received_Q;                       /* Received signals/branch */

int sector;

double transf_I;                       /* I-coordinate of sector  */
double transf_Q;                       /* Q-coordinate of sector  */

double buffer_I[1024];
double buffer_Q[1024];
int buffer_sector[1024];

/* int fflush(); /* */

/* Data structures used for trellis sections and survivors */

struct trel {
  int init;                /* initial state */
  int data;                /* data symbol */
  int final;               /* final state */
  int output;              /* index of output coded symbols */
}; 

struct surv {
  double metric;           /* metric */
  int data[1024];  /* estimated data symbols */
  int state[1024];  /* state sequence */
};

/* A trellis section is an array of branches, indexed by an initial
   state and a k2-bit input data. The values that can be read
   are the final state and the index of the output symbols */

struct trel trellis[1024][64];

/* A survivor is a sequence of states and estimated data, of length
   equal to TRUNC_LENGTH, together with its corresponding metric.
   A total of NUM_STATES survivors are needed */

struct surv survivor[1024], surv_temp[1024];

/* Function prototypes */

void build_sector_table(void); /* Sector look-up table */
void encoder2(void);           /* Encoder in transmitter */
void encoder(void);            /* Encoder in receiver */
int random_data(void);         /* Random data generator */
void transmit(void);           /* Encoder & BPSK modulator */
void awgn(void);               /* AWGN generator */
void find_sect(void);          /* Cartesian to polar conversion */
void transform(void);          /* Find coordinates QPSK sector */
void viterbi(void);            /* Viterbi decoder */
void re_encode(void);          /* Reencoder of decoded information bit */
void sector_table(void);       /* Construct branch metric look-up table */
void open_files(void);         /* Open files for output */
void close_files(void);        /* Close files */

main(int argc, char *argv[])
{

  register int i, j, k;
  int init, data, final, output;
  register int error;
  unsigned long error_count, error_count_1, error_count_2;
  FILE  *fp_ber;             /* Pointer for overall BER data file */
  double ii;

  // Command line processing
  if (argc != 9)
    {
      printf("Usage %s file_input file_output truncation snr_init snr_final snr_inc num_sim seed\n", argv[0]);
      exit(0);
    }

  sscanf(argv[1],"%s", name1);
  sscanf(argv[2],"%s", name2);
  sscanf(argv[3],"%ld", &TRUNC_LENGTH);
  sscanf(argv[4],"%lf", &INIT_SNR);
  sscanf(argv[5],"%lf", &FINAL_SNR);
  sscanf(argv[6],"%lf", &SNR_INC);
  sscanf(argv[7],"%ld", &NUMSIM);
  sscanf(argv[8],"%ld", &seed);

  printf("\nSimulation of TCM decoding with 8-PSK modulation over an AWGN channel\n");
  printf("%ld simulations per Eb/No (dB) point\n", NUMSIM);

  fp_ber = fopen(name2,"w");

  /* Open file with trellis data */
  if (!(fp = fopen(name1,"r"))) 
    {
      printf("Error opening file!\n");
      exit(1);
    }

  fscanf(fp, "%d %d", &n2, &m2);

  RATE = (2.0 / (double) n2)*3.0;  // 2 bits per symbol if encoder is rate 1/2

  fscanf(fp, "%d %d %d", &NUM_STATES, &OUT_SYM, &NUM_TRANS);
  for (j=0; j<n2; j++)
    fscanf(fp, "%x", &g2[j][0]);

  printf("\n%d-state rate-1/%d binary convolutional encoder\n",
          NUM_STATES, n2);
  printf("with generator polynomials ");
  for (j=0; j<n2; j++) printf("%x ", g2[j][0]); printf("\n");
  printf("\nDecoding depth = %ld\n\n", TRUNC_LENGTH);


    /* =================== READ TRELLIS STRUCTURE ==================== */

      for (j=0; j<NUM_STATES; j++) /* For each state in the section */
	    for (k=0; k<NUM_TRANS; k++) /* and for each outgoing branch */
	      {
	        /* Read initial state, input data and final state */
	        fscanf(fp,"%d %d %d", &trellis[j][k].init, &trellis[j][k].data,
		       &trellis[j][k].final);
	        /* Read the output symbols of the branch */
		    fscanf(fp,"%d", &data);
		    trellis[j][k].output = data;
	      } /* end for states */
        /* end for branches */
      
      fclose(fp);



  snr = INIT_SNR;

  /* Pragmatic bits-to-8PSK signal mapping


                                |
                     [3] 011 #  |  * 001 [1]
                                |
                                |
                [2] 010 .       |       + 000 [0]
                                |
                   ---------------------------
                                |
                [4] 100 +       |       . 110 [6]
                                |
                                |
                     [5] 101 *  |  # 111 [7]
                                |

   */

  /* Coordinates of the 8-PSK symbols */

  psk_I[0] = cos(M_PI/8.0);        psk_Q[0] = sin(M_PI/8.0);
  psk_I[1] = cos(3.0*M_PI/8.0);    psk_Q[1] = sin(3.0*M_PI/8.0);
  psk_I[3] = -cos(3.0*M_PI/8.0);   psk_Q[3] = sin(3.0*M_PI/8.0);
  psk_I[2] = -cos(M_PI/8.0);       psk_Q[2] = sin(M_PI/8.0);
  psk_I[4] = -cos(M_PI/8.0);       psk_Q[4] = -sin(M_PI/8.0);
  psk_I[5] = -cos(3.0*M_PI/8.0);   psk_Q[5] = -sin(3.0*M_PI/8.0);
  psk_I[7] = cos(3.0*M_PI/8.0);    psk_Q[7] = -sin(3.0*M_PI/8.0);
  psk_I[6] = cos(M_PI/8.0);        psk_Q[6] = -sin(M_PI/8.0);

  /* Coordinates of the QPSK coset constellation */

  qpsk_I[0] = cos(M_PI/4.0);        qpsk_Q[0] = sin(M_PI/4.0);
  qpsk_I[1] = -cos(M_PI/4.0);       qpsk_Q[1] = sin(M_PI/4.0);
  qpsk_I[3] = -cos(M_PI/4.0);       qpsk_Q[3] = -sin(M_PI/4.0);
  qpsk_I[2] = cos(M_PI/4.0);        qpsk_Q[2] = -sin(M_PI/4.0);

  /* Build the sector look-up table for estimating uncoded bit */

  build_sector_table();

#ifdef SECTOR
  printf("s[0] = %x, s[1] = %x, s[3] = %x, s[2] = %x\n",s[0],s[1],s[3],s[2]);
#endif

  /* ======================== SNR LOOP ============================= */

  while ( snr < (FINAL_SNR+0.001) )

    {

  srandom(seed);

  amp = sqrt(2.0*RATE*pow(10.0,(snr/10.0)));

  /* Initialize transmitted data sequence */
      
  for (i=0; i<TRUNC_LENGTH; i++)
    {
      data_symbol[i]=0;
      data_symbol2[i]=0;
      buffer_sector[i] = 0;
    }

  /* Initialize survivor sequences and metrics */
      
  for (i=0; i<NUM_STATES; i++)
	{
	  survivor[i].metric = 0.0;             /* Metric = 0 */

	  for (j=0; j<TRUNC_LENGTH; j++)
	    {
	      survivor[i].data[j] = 0;        /* Estimated data = 0 */
	      survivor[i].state[j] = 0;       /* Estimated state = 0 */
	    }
	}
      
  /* Index used in simulation loop */

  indxx = 0;
      
  /* Initialize encoder memories */

  memory2 = 0;
  memory = 0;

  /* Error counters */

  error_count = 0;
  error_count_1 = 0;
  error_count_2 = 0;

  /* ===================== SIMULATION LOOP ========================= */
      
  while (indxx < NUMSIM) 

	{ 
	
	  /* GENERATE RANDOM TWO-BIT INFORMATION SYMBOLS */
	  
      i = random_data();
	  data_symbol[indxx % TRUNC_LENGTH] = (i & 1);          /* BIT 1 */
	  data_symbol2[indxx % TRUNC_LENGTH] = ( (i>>1) & 1 );  /* BIT 2 */
	  
#ifdef PRINT
	  printf("Transmitted data sequence:\n");
	  printf("%2d %x   ",(indxx % TRUNC_LENGTH),
		 data_symbol[indxx % TRUNC_LENGTH]);
	  for (i=0; i<TRUNC_LENGTH; i++)
	    printf("%x", data_symbol[i]);
	  printf("\n");
#endif

	  
	  /* ENCODE AND MODULATE data bits */
	  
	  transmit();
	  
#ifdef PRI
	  printf("Transmitted = ");
	  printf("%d ", transmitted);
	  printf("\n");
#endif
	  
	
	  /* ADD ADDITIVE WHITE GAUSSIAN NOISE */
	  
	  awgn(); /* */
	  


      /* DETERMINE THE PHASE OF THE RECEIVED 8-PSK SYMBOL */

      find_sect();

      /* Store sector in buffer */

      buffer_sector[indxx % TRUNC_LENGTH] = sector;



      /* TRANSFORMATION (I_in,Q_in) to (I,Q) */

      transform();



	  /* VITERBI DECODE (FOR THE ENCODED INFORMATION BIT) */
	  
	  viterbi();



	  indxx += 1;           /* Increment simulation index */



      /* ESTIMATE UNCODED INFORMATION BIT */

      /* Re-encode estimated encoded information bit, BIT 1 */

      re_encode();


      /* "reconstructed" is the estimated coset (output of re-encoder) */

      estimate_data2 = (s[reconstructed]>>buffer_sector[indxx % TRUNC_LENGTH])
                       & 0x01;

#ifdef PRINT
printf("buffer_I:\n");
for (i=0; i<TRUNC_LENGTH; i++)
printf("%f", buffer_I[i]);
printf("\n");
printf("coset = %d, estimate_data2 = %d\toriginal data2 = %d\n", 
reconstructed, estimate_data2, data_symbol2[indxx % TRUNC_LENGTH]);
#endif


	  /* COMPUTE ERRORS */
	  
	  error = (survivor[0].data[indxx % TRUNC_LENGTH]
	          ^ data_symbol[indxx % TRUNC_LENGTH]) & 0x01; 
  
      error_count_1 += error;

	  error += (estimate_data2 ^ data_symbol2[indxx % TRUNC_LENGTH]); 

      error_count_2 += (estimate_data2 ^ data_symbol2[indxx % TRUNC_LENGTH]);

	  if (error)
	    error_count+=error;
	  
#ifdef SHOW_PROGRESS
	  if ( (indxx % DELTA) == 0 )
	    {
        printf("index = %ld\n", indxx);
	      for (i=0; i<TRUNC_LENGTH; i++)
		printf("%x", survivor[0].data[i]);
	      printf("\n--- errors = %ld %ld %ld\n", error_count,
          error_count_1, error_count_2);
	    }
#endif
	  
	}
      
    /* PRINT SIMULATION RESULTS */

      printf("%10.5f  %10.4e  %10ld %10ld %10ld %10ld\n",
             snr, ( (double) error_count / ((double) indxx * 2.0) ),
             indxx, error_count, error_count_1, error_count_2);
      
    fprintf(fp_ber, "%f %20.15f %20.15f %20.15f\n", snr,
	      ( (double) error_count / ((double) indxx * 2.0) ),
	      ( (double) error_count_1 / ((double) indxx) ),
	      ( (double) error_count_2 / ((double) indxx) ));

    fflush(stdout);
    fflush(fp_ber);

    snr += SNR_INC;

  }
  
  fclose(fp_ber);

  return 0;

}



void build_sector_table()
{

  /*
   * ---------------------------------------------------------------------
   *            Look-up table for decoding BIT 2
   * ---------------------------------------------------------------------
   */

                /*                Sector      */
                /* -------------------------  */
                /* COSET     7 6 5 4 3 2 1 0  */
                /* -------------------------  */
  s[0] = 0x78;  /*  00       0 1 1 1 1 0 0 0  */
  s[1] = 0xf0;  /*  01       1 1 1 1 0 0 0 0  */
  s[3] = 0xe1;  /*  11       1 1 1 0 0 0 0 1  */
  s[2] = 0xc3;  /*  10       1 1 0 0 0 0 1 1  */
}



int random_data()
{
  /* 
   * ---------------------------------------------------------------------
   *            Two-bit random number generator 
   * ---------------------------------------------------------------------
   */
  return( (random() >> 5) & 0x03 );
}



void encoder2()
{
  /* 
   * ---------------------------------------------------------------------
   *      Convolutional encoder, rate 1/n2 (fixed k_2=1) 
   * ---------------------------------------------------------------------
   */
  register int i, j, result, temp;
  
  temp = memory2;
  output = 0;

  temp = (temp<<1) ^ data_symbol[indxx % TRUNC_LENGTH];

  for (i=0; i<n2; i++)
   {
     result = 0;
     for (j=m2; j>=0; j--)
       result ^= ( ( temp & g2[i][0] ) >> j ) & 1;
     output = ( output<<1 ) ^ result;
   }
  memory2 = temp ;
}



void encoder()
{
  /* 
   * ---------------------------------------------------------------------
   *       Convolutional encoder, rate 1/n2 (fixed k_2=1) 
   * NOTE: Used in re-encoding coset indexes
   * ---------------------------------------------------------------------
   */
  register int i, j, result, temp;
  
  temp = memory;
  output2 = 0;

  temp = (temp<<1) ^ survivor[0].data[indxx % TRUNC_LENGTH];

  for (i=0; i<n2; i++)
   {
     result = 0;
     for (j=m2; j>=0; j--)
       result ^= ( ( temp & g2[i][0] ) >> j ) & 1;
     output2 = ( output2<<1 ) ^ result;
   }
  memory = temp ;
}



void transmit()
{
  /* 
   * ---------------------------------------------------------------------
   *        Encode and modulate 2-bit data sequence
   * ---------------------------------------------------------------------
   */
  int i;

  /* Encode BIT 1 */

  encoder2();

  /* Transmitted is an index. BIT 2 determines opposite 8-PSK symbols */

  transmitted = output + 4*data_symbol2[indxx % TRUNC_LENGTH];

}


void re_encode()
{
  /* 
   * ----------------------------------------------------------------------
   * Re-encode information bit (BIT 1) from Viterbi decoder to obtain coset
   * ----------------------------------------------------------------------
   */

  encoder();
  reconstructed = output2;
}


void awgn()
{
  /* 
   * ---------------------------------------------------------------------
   * Add AWGN to transmitted sequence
   * ---------------------------------------------------------------------
   */
  double u1,u2,s,noise,randmum;
  int i;

#ifdef PRINT
      printf("Received    = ");
#endif

  do
	{	
	  randmum = (double)(random())/MAX_RANDOM;
	  u1 = randmum*2.0 - 1.0;
	  randmum = (double)(random())/MAX_RANDOM;
	  u2 = randmum*2.0 - 1.0;
	  s = u1*u1 + u2*u2;
	} while( s >= 1);
  noise = u1 * sqrt( (-2.0*log(s))/s ) / amp; /* */
#ifdef NO_NOISE
  noise = 0.0;
#endif
  received_I = psk_I[transmitted] + noise;

  do
    {
      randmum = (double)(random())/MAX_RANDOM;
      u1 = randmum*2.0 - 1.0;
      randmum = (double)(random())/MAX_RANDOM;
      u2 = randmum*2.0 - 1.0;
      s = u1*u1 + u2*u2;
    } while( s >= 1);
  noise = u1 * sqrt( (-2.0*log(s))/s ) / amp; /* */
#ifdef NO_NOISE
  noise = 0.0;
#endif
  received_Q = psk_Q[transmitted] + noise;

#ifdef PRINT
  printf("%4.2lf %4.2lf ", received_I, received_Q);
  printf("\n");
#endif
}




void viterbi()
{
  /* 
   * ---------------------------------------------------------------------
   *                         Viterbi decoder
   * ---------------------------------------------------------------------
   */

  double aux_metric, surv_metric[NUM_STATES];
  register int i,j,k;

  for (i=0; i<NUM_STATES; i++) /* Initialize survivor branch metric */
    surv_metric[i] = DBL_MAX;

  for (i=0; i<NUM_STATES; i++) /* Loop over inital states */
    {
    for (j=0; j<NUM_TRANS; j++) /* Loop over data */
      {

	  /* Compute metric between received and coded symbols */

      aux_metric = ( (qpsk_I[trellis[i][j].output]-transf_I)
                     * (qpsk_I[trellis[i][j].output]-transf_I) )
                   +
                   ( (qpsk_Q[trellis[i][j].output]-transf_Q)
                       * (qpsk_Q[trellis[i][j].output]-transf_Q) );

	  aux_metric += survivor[i].metric;

#ifdef PRINT
	  printf("rec_I   rec_Q  index_coded\n");
	  printf("%6.2lf %6.2lf %d\n", received_I,received_Q,trellis[i][j].output);
	  printf("init, data, final, metric = %2d %2d %2d %lf\n", i, j,
		 trellis[i][j].final, aux_metric);
#endif

	  /* compare with survivor metric at final state */
      if ( aux_metric < surv_metric[trellis[i][j].final] )

	    { /* Good candidate found */

	      surv_metric[trellis[i][j].final] = aux_metric;

	      /* Update data and state paths */
	      for (k=0; k<TRUNC_LENGTH; k++)
		    {
		      surv_temp[trellis[i][j].final].data[k] =
		        survivor[i].data[k];
		      surv_temp[trellis[i][j].final].state[k] =
		        survivor[i].state[k];
		    }

          surv_temp[trellis[i][j].final].data[indxx%TRUNC_LENGTH] = j;
          surv_temp[trellis[i][j].final].state[indxx%TRUNC_LENGTH] = 
		          trellis[i][j].final;
	    }
      }
    }
  
  for (i=0; i<NUM_STATES; i++)
    {
      /* Update survivor metrics */
      survivor[i].metric = surv_metric[i];

      for (k=0; k<TRUNC_LENGTH; k++)
        {
          survivor[i].data[k] = surv_temp[i].data[k];
          survivor[i].state[k] = surv_temp[i].state[k];
        }
    }
  
#ifdef PRINT
  for (i=0; i<NUM_STATES; i++)
    {
      printf("survivor %2d, initial state = %2d metric = %10.5lf, ", i,
	     survivor[i].state[(indxx-1) % TRUNC_LENGTH], survivor[i].metric);
      printf("\t estimate = %d\n", survivor[i].data[indxx % TRUNC_LENGTH]);
    }
  printf("\n");
#endif

#ifdef DEBUG
  printf("Surviving paths (data):\n");
   for (i=0; i<NUM_STATES; i++)
     {
       printf("%2d --> ",i);
       for (k=0; k<TRUNC_LENGTH; k++)
	 printf("%1x", survivor[i].data[k]);
       printf("\n");
     }
  printf("\n");
#endif

}

void find_sect()
{
  /* ---------------------------------------------------------------------
   * Determine the sector in which the received point lies
   * ---------------------------------------------------------------------
   */

  register int i, flag;
  double ii, received_angle, phase1, phase2;

  flag = 0;
  received_angle = atan2(received_Q, received_I);

  if (received_angle < 0.0)
    received_angle += (2.0*M_PI);

  if ( ( received_angle > 2.0*M_PI-M_PI/8 ) || (received_angle < M_PI/8.0) )
    sector = 0;
  else
  for (i=1; ( (i<8) && (!flag)); i++)
    {
      ii = (double) i;
      phase1 = (ii*M_PI/4.0) - (M_PI/8.0);
      phase2 = (ii*M_PI/4.0) + (M_PI/8.0);
#ifdef DEBUG
printf("phase1, phase2 = %lf %lf\n", (180./M_PI)*phase1,(180./M_PI)*phase2);
#endif
      if ( (received_angle > phase1) && (received_angle <= phase2) )
        {
          flag = 1;
          sector = i;
        }
    }

#ifdef SECTOR
printf("---\nrec_I, rec_Q = %lf, %lf\nangle = %lf, sector = %d\n",
    received_I, received_Q, (180./M_PI)*received_angle, sector);
#endif
}


void transform()
{

  double amplitude, phase;

  amplitude = sqrt (received_I*received_I + received_Q*received_Q);
  phase = 2.0*atan(received_Q/received_I);

  transf_I = amplitude * cos(phase);
  transf_Q = amplitude * sin(phase);

#ifdef IQ
printf("rec_I, rec_Q = (%lf, %lf) \ntransf_I, transf_Q = (%lf, %lf)\n",
  received_I,received_Q,transf_I,transf_Q);
#endif
}