// Description: Viterbi decoder, hard decision decoding. For a
// rate-1/n convolutional code. Hamming metric.
//
// Maximum of 1024 states; max. truncation length = 1024
//
#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 NUM_STATES, OUT_SYM, NUM_TRANS;
long TRUNC_LENGTH;
double RATE;
double INIT_SNR, FINAL_SNR, SNR_INC;
long NUMSIM;
FILE *fp, *fp_ber;
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 name1[40], name2[40];
// Command line processing
if (argc != 8)
{
printf("Usage %s file_input file_output truncation snr_init snr_final snr_inc num_sim\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);
printf("\nSimulation of Viterbi decoding over an AWGN channel with hard decisions\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);
fprintf(stderr,"SCAN n2=%d m2=%d\n",n2,m2);
RATE = 1.0 / (double) n2;
fprintf(stderr,"RATE: %f\n",RATE);
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("\n\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
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;
}
fprintf(stderr,"mem2: %d\n",memory2);
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--)
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 HAMMING DISTANCE*/
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)
{ // HAMMING DISTANCE between received and reference values
// This is equivalent to hard decision + Hamming distance computation:
if ( (rec<0.0) && (ref<0.0) ) return(0.0);
if ( (rec<0.0) && (ref>0.0) ) return(1.0);
if ( (rec>0.0) && (ref<0.0) ) return(1.0);
if ( (rec>0.0) && (ref>0.0) ) return(0.0);
}
19.04.2012, 10:53 