viterbi_binary_punct.c

error correction code

// ------------------------------------------------------------------------
// File: viterbi_binary_punct.c
// Date: April 1, 2002
// Description: Viterbi decoder, maximum likelihood decoding. For a
//              rate-1/n convolutional code. Correlation metric.
//
// Maximum of 1024 states; max. truncation length = 1024
//
// The puncturing rule must be in a separate file.
//
// CAUTION: The row order of the puncturing rule is reversed from that
//          in the book. Specify the rule starting from the last row
//          down to the first one. See the example files in this
//          directory: "rule_kn.data" specifies the puncturing rule for
//          a rate-k/n code derived from a rate-1/2 convolutional code.
//
//          If you give the wrong puncturing rule, the results will be
//          evident. Try it out. 
//
// ------------------------------------------------------------------------
// 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>
#include <time.h>

#define MAX_RANDOM LONG_MAX 

// Use the compiling directive -DSHOW_PROGRESS to see how the decoder
// converges to the decoded sequence by monitoring the survivor memory
#ifdef SHOW_PROGRESS
#define DELTA 1
#endif

int k2=1, n2, m2;
int length;
int NUM_STATES, OUT_SYM, NUM_TRANS;

long TRUNC_LENGTH;
int punct[10][10];                       // Puncturing rule

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

FILE *fp, *fp_ber, *fp_punct;

int g2[10][10];
unsigned int memory2, output;            /* Memory and output */
unsigned int data2;                      /* Data */

unsigned long seed;                      /* Seed for random generator */
unsigned int data_symbol[1024];          /* 1-bit data sequence */
long index;                              /* Index for memory management*/
double transmitted[10];                  /* Transmitted signals/branch */
double snr, amp;
double received[10];                     /* Received signals/branch */

int fflush();

// Data structures used for trellis sections and survivors
struct trel {
  int init;                /* initial state */
  int data;                /* data symbol */
  int final;               /* final state */
  double output[10];  /* output coded symbols (branch label) */
}; 
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 k_2-bit input data. The values read
//   are the final state and the output symbols 
struct trel trellis[1024][100];

/* 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 for C_{O2} */
int random_data(void);         /* Random data generator */
void transmit(void);           /* Encoder & BPSK modulator */
void awgn(void);               /* Add AWGN */
void viterbi(void);            /* Viterbi decoder */
double comp_metric(double rec, double ref); /* Metric calculator */
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;
char name[40], name1[40], name2[40];

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

  sscanf(argv[1],"%s", name1);
  sscanf(argv[2],"%s", name);
  sscanf(argv[3],"%s", name2);
  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);

printf("\nSimulation of Viterbi decoding 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 = 1.0 / (double) n2;

  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 */
	for (i=0; i < OUT_SYM; i++)
	  {
             fscanf(fp,"%d", &data);
             trellis[j][k].output[i] = (double) data;
	  }
      } /* end for states */
      /* end for branches */
  
  fclose(fp);
  
#ifdef PRINT
  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 */
      {
	printf("%3d %3d %3d  ", trellis[j][k].init,
	trellis[j][k].data, trellis[j][k].final);
	for (i=0; i < OUT_SYM; i++)
	  printf("%2.0lf ", trellis[j][k].output[i]);
	printf("\n");
       } /* end for states */
      /* end for branches */
#endif


  /* Open file with puncturing rule */
  if (!(fp_punct = fopen(name,"r")))
    {
    printf("Error opening file!\n");
    exit(1);
    }

  fscanf(fp_punct,"%d", &length);

  for (i=0; i<n2; i++)
    for (j=0; j<length; j++)
      fscanf(fp_punct, "%d", &punct[i][j]);

  // Compute rate of the code
  RATE = 0.0;
  for (i=0; i<n2; i++)
    for (j=0; j<length; j++) if (punct[i][j]) RATE += 1.0;
  RATE = (double)(length)/RATE;

  printf("This is a rate-%7.6lf punctured convolutional code\n", RATE);
  printf("Puncturing period = %d\n", length);
  printf("Puncturing rule:\n");
  for (i=0; i<n2; i++)
    {
    for (j=0; j<length; j++)
      printf("%d ", punct[i][j]);
    printf("\n");
    }
  printf("\n");

  snr = INIT_SNR;

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

  while ( snr < (FINAL_SNR+1.0e-6) )
    {
      amp = sqrt( 2.0 * RATE * pow(10.0, (snr/10.0)) );

      /* Random seed from current time */
      time(&seed);
      srandom(seed);
      
      /* Initialize transmitted data sequence */
      
      for (i=0; i<TRUNC_LENGTH; i++)
	data_symbol[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 */
      index = 0;
      
      /* Initialize encoder memory */
      memory2 = 0;

      /* Error counter */
      error_count = 0;


      /* ===================== SIMULATION LOOP ========================= */
      
      while (index < NUMSIM) 
	{ 
	
	  /* GENERATE a random bit */
	  data_symbol[index % TRUNC_LENGTH] = random_data(); /* */
	  
	  /* ENCODE AND MODULATE (BPSK) data bit */
	  transmit();
	  
	  /* ADD ADDITIVE WHITE GAUSSIAN NOISE */
	  awgn(); 

	  /* VITERBI DECODE */
	  viterbi();
	  
	  index += 1;           /* Increase simulation index */
	  
	  /* COMPUTE ERRORS */
	  error = survivor[0].data[index % TRUNC_LENGTH]
	                              ^ data_symbol[index % TRUNC_LENGTH]; 
	  if (error)
	    error_count++;
	  
#ifdef SHOW_PROGRESS
	  if ( (index % DELTA) == 0 )
	    {
              for (j=0; j<NUM_STATES; j++) 
                {
	          printf("%3d %2d   ", (index % TRUNC_LENGTH), j);
	          for (i=0; i<TRUNC_LENGTH; i++)
		    printf("%x", survivor[j].data[i]);
	          printf(" %ld\n", error_count);
                }
	    }
#endif
	}
      
      printf("%f  %10.4e\n",snr, ((double) error_count/(double) index) );
      
      fprintf(fp_ber, "%f %10.4e\n", snr, ((double)error_count /(double)index));
	      

      fflush(stdout);
      fflush(fp_ber);

      snr += SNR_INC;
      
    }
  
  fclose(fp_ber);
}





int random_data()
{
  /* Random bit generator */
  return( (random() >> 5) & 1 );
}







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

  encoder2(); /* */

  /* Modulate: {0,1} --> {+1,-1} */
  for (i=(OUT_SYM-1); i >=0; i--)
    if ( (output >> i) & 0x01 )
      transmitted[OUT_SYM-1-i] = -1.0;  /* if bit = 1 */
    else
      transmitted[OUT_SYM-1-i] =  1.0;  /* if bit = 0 */
}






void encoder2()
{
  // Binary convolutional encoder, rate 1/n2 
  register int i, j, result, temp;
  
  temp = memory2;
  output = 0;

  temp = (temp<<1) ^ data_symbol[index % 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 awgn()
{
  /* Add AWGN to transmitted sequence */
  double u1,u2,s,noise,randmum;
  int i;

  for (i=0; i<OUT_SYM; i++)
    {

      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 );
      received[i] = transmitted[i] + noise/amp;
    }
}








void viterbi()
{
  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 CORRELATION between received seq. and coded branch */
	  aux_metric = 0.0;
	  for (k=(OUT_SYM-1); k>=0; k--)
            {
            // Computation only for unpunctured positions:
            if (punct[k][(index%length)])
	      aux_metric += comp_metric(received[k],trellis[i][j].output[k]);
            }
	  aux_metric += survivor[i].metric;

	  /* compare with survivor metric at final state */
	  /* We compare CORRELATIONS */

          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[index%TRUNC_LENGTH] = j;
		  surv_temp[trellis[i][j].final].state[index%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];
	}
    }

}







double comp_metric(double rec, double ref)
{ /* CORRELATION between received and reference values, to reduce the
     computational effort */

/*  double aux; /* */
/*  aux = (rec-ref)*(rec-ref); /* */

  if (ref == 1.0)
    return(rec);
  else
    return(-rec);
}