//
------------------------------------------------------------------------
// File: viterbi_binary_nsc.c
// Date: April 1, 2002
// Description: Viterbi decoder, maximum likelihood
decoding. For a
// rate-1/n
convolutional code. Correlation metric.
//
// NOTE: The generator polynomials are given in
hexadecimal and in
// reversed
order. For example, the rate-1/2 standard code
// has
generator polynomials (octal) (g1,g2)=(133,171). Then
//
// Base
8 Reverse Base 16
// 133 =
1011011 ---> 1101101 = 6d
// 171 =
1111001 ---> 1001111 = 4f
//
// These
values can be found in file "trellis_64states.data"
//
// Also, the
order in the file is (g2,g1). That is, the last
// generator
in the first row, down to the first generator in
// the last
row.
//
// With this
format, the software implementation of the
// algorithm
is simplified.
//
// NOTE: 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;
unsigned int data2;
unsigned long seed;
unsigned int data_symbol[1024];
long index;
double transmitted[10];
double snr, amp;
double received[10];
int fflush();
// Data structures used for trellis sections and
survivors
struct trel {
int init;
int data;
int final;
double
output[10];
};
struct surv {
double
metric;
int
data[1024];
int
state[1024];
};
// 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];
struct surv survivor[1024], surv_temp[1024];
void encoder2(void);
int random_data(void);
void transmit(void);
void awgn(void);
void viterbi(void);
double comp_metric(double rec, double ref);
void open_files(void);
void close_files(void);
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\n");
printf("%ld simulations per Eb/No (dB)
point\n", NUMSIM);
fp_ber =
fopen(name2,"w");
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
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
printf("\n\nDecoding depth = %ld\n\n", TRUNC_LENGTH);
for (j=0;
j
for (k=0;
k
{
fscanf(fp,"%d %d %d", &trellis[j][k].init,
&trellis[j][k].data,
&trellis[j][k].final);
for (i=0; i
< OUT_SYM; i++)
{
fscanf(fp,"%d", &data);
trellis[j][k].output[i] = (double) data;
}
}
fclose(fp);
#ifdef PRINT
for (j=0;
j
for (k=0;
k
{
printf("= = = ",
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");
}
#endif
snr = INIT_SNR;
while ( snr <
(FINAL_SNR+1.0e-6) )
{
amp = sqrt(
2.0 * RATE * pow(10.0, (snr/10.0)) );
time(&seed);
srandom(seed);
for (i=0;
i
data_symbol[i]=0;
for (i=0;
i
{
survivor[i].metric = 0.0;
for (j=0;
j
{
survivor[i].data[j] = 0;
survivor[i].state[j] = 0;
}
}
index = 0;
memory2 = 0;
error_count =
0;
while (index
< NUMSIM)
{
data_symbol[index % TRUNC_LENGTH] = random_data();
transmit();
awgn();
viterbi();
index +=
1;
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
{
printf("= - ",
(index % TRUNC_LENGTH), j);
for (i=0; i
printf("%x", survivor[j].data[i]);
printf(" %ld\n", error_count);
}
}
#endif
}
printf("%f .4e\n",snr,
((double) error_count/(double) index) );
fprintf(fp_ber, "%f .4e\n", snr, ((double)error_count
/(double)index));
fflush(stdout);
fflush(fp_ber);
snr +=
SNR_INC;
}
fclose(fp_ber);
}
int random_data()
{
return( (random()
>> 5) & 1 );
}
void transmit()
{
int i;
encoder2();
for
(i=(OUT_SYM-1); i >=0; i--)
if ( (output
>> i) & 0x01 )
transmitted[OUT_SYM-1-i] = -1.0;
else
transmitted[OUT_SYM-1-i] =
1.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
{
result = 0;
for (j=m2;
j>=0; j--)
result ^= (
( temp & g2[i][0] ) >> j ) & 1;
output = ( output<<1 ) ^ result;
}
memory2 = temp ;
}
void awgn()
{
double
u1,u2,s,noise,randmum;
int i;
for (i=0;
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
surv_metric[i]
= -DBL_MAX;
for (i=0;
i
{
for (j=0;
j
{
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;
if ( aux_metric >
surv_metric[trellis[i][j].final] )
{
surv_metric[trellis[i][j].final] = aux_metric;
for
(k=0; 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
{
survivor[i].metric = surv_metric[i];
for (k=0;
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)
{
if (ref == 1.0)
return(rec);
else
return(-rec);
}