tc-8psk-sect.c

The same two-stage decoder as above. However, when transforming the symbols prior to Viterbi decodin

/*
 * File: tc-8psk-sect.c
 * Date: January 8, 1998
 * Description: Suboptimal trellis decoder using sector information. 
 */
// ------------------------------------------------------------------------
// 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];
FILE  *fp;                 /* Pointer for trellis data file */

#define MAX_RANDOM LONG_MAX /* Maximum value of random() */

long NUM_SECT;

//#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 */

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

double received_angle;
double angle_inc;
int sector;
int buffer_sector[1024];
double branch_metric[512][8];

/* 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 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 viterbi(void);            /* Viterbi decoder */
void re_encode(void);          /* Reencoder of decoded information bit */
void sector_table();           /* 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;
  FILE  *fp_ber;             /* Pointer for overall BER data file */
  double RATE;
  double ii, reference_I[512], reference_Q[512];

  // Command line processing
  if (argc != 10)
    {
      printf("Usage %s file_input file_output NUM_SECT 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", &NUM_SECT);
  sscanf(argv[4],"%ld", &TRUNC_LENGTH);
  sscanf(argv[5],"%lf", &INIT_SNR);
  sscanf(argv[6],"%lf", &FINAL_SNR);
  sscanf(argv[7],"%lf", &SNR_INC);
  sscanf(argv[8],"%ld", &NUMSIM);
  sscanf(argv[9],"%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", TRUNC_LENGTH);
  printf("Number of sectors = %ld\n\n", NUM_SECT);


  /* =================== 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;
  angle_inc = 2.0*M_PI/(double)NUM_SECT;

  /* Bits-to-8PSK signal mapping as used in Japanese Satellite ISDB */

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

  sector_table(); /* Branch metric lookup table */

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

  while ( snr < (FINAL_SNR+0.001) )

    {

  /* ==================== INITIALIZE SEQUENCES ===================== */
      
  /* Random seed from current time */
  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;

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

	{ 
	
	  /* GENERATE random two-bit symbols */
	  
      i = random_data();
	  data_symbol[indxx % TRUNC_LENGTH] = (i & 1); /* */
	  data_symbol2[indxx % TRUNC_LENGTH] = ( (i>>1) & 1 ); /* */
	  
#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 (BPSK) data bit */
	  
	  transmit();
	  
#ifdef PRI
	  printf("Transmitted = ");
	  printf("%d ", transmitted);
	  printf("\n");
#endif
	  
	
	  /* ADD ADDITIVE WHITE GAUSSIAN NOISE */
	  
	  awgn(); /* */
	  

      /* CARTESIAN TO POLAR COORDINATE */

      find_sect();

      /* Store sector in buffer */

      buffer_sector[indxx % TRUNC_LENGTH] = sector;


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

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

      /* ESTIMATE UNCODED INFORMATION BIT */

      re_encode();

      estimate_data2 = 0;

#ifdef HARD_DECISION
      /* The half-plane defines uncoded bit */
      if (buffer_sector[indxx%TRUNC_LENGTH] > NUM_SECT/2)
        estimate_data2 = 1;
#else
      if ( ( pow((cos(buffer_sector[indxx%TRUNC_LENGTH]*angle_inc)
                      - psk_I[reconstructed]),2.0)
           + pow((sin(buffer_sector[indxx%TRUNC_LENGTH]*angle_inc)
                      - psk_Q[reconstructed]),2.0))
           > ( pow((cos(buffer_sector[indxx%TRUNC_LENGTH]*angle_inc)
                      - psk_I[reconstructed+4]),2.0)
            + pow((sin(buffer_sector[indxx%TRUNC_LENGTH]*angle_inc)
                      - psk_Q[reconstructed+4]),2.0)) )
        estimate_data2 = 1;
#endif

#ifdef ABS
      estimate_data2 = 0;
      if ( ( fabs(cos(buffer_sector[indxx%TRUNC_LENGTH]*angle_inc)
                      - psk_I[reconstructed])
           + fabs(sin(buffer_sector[indxx%TRUNC_LENGTH]*angle_inc)
                      - psk_Q[reconstructed]))
           > ( fabs(cos(buffer_sector[indxx%TRUNC_LENGTH]*angle_inc)
                      - psk_I[reconstructed+4])
            + fabs(sin(buffer_sector[indxx%TRUNC_LENGTH]*angle_inc)
                      - psk_Q[reconstructed+4])) )
        estimate_data2 = 1;
#endif

#ifdef PRINT
printf("Sector buffer:\n");
for (i=0; i<TRUNC_LENGTH; i++)
printf("%d", buffer_sector[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];