#include <stdio.h>
#include <math.h>
#include <stdlib.h>

/*
Convolutional binary rate 1/3  nonsystematic code
Dfree=16
K=7  (trellis length = 8)

Connection vectors (from K. J. Larsen):

D0 = x7 + x4 + x2 + 1                   (0x95)
D1 = x7 + x6 + x4 + x3 + 1              (0xd9)
D2 = x7 + x6 + x5 + x4 + x2 + x1 + 1    (0xf7)
*/

#define WINDOW_BYTES        8           //  Parameter - decoding window length in bytes


typedef struct {
               float cost[256];
               unsigned char path[256][WINDOW_BYTES];
               } DECODER_STATUS;

typedef struct {
               unsigned char id;
               unsigned char path_no[2];
               unsigned char path_continuation[2];
               float path_cost[2];
               } NODE;



unsigned char const trellis_to_tribit[256] =
{
 0,7,1,6,5,2,4,3,2,5,3,4,7,0,6,1,
 7,0,6,1,2,5,3,4,5,2,4,3,0,7,1,6,
 1,6,0,7,4,3,5,2,3,4,2,5,6,1,7,0,
 6,1,7,0,3,4,2,5,4,3,5,2,1,6,0,7,
 3,4,2,5,6,1,7,0,1,6,0,7,4,3,5,2,
 4,3,5,2,1,6,0,7,6,1,7,0,3,4,2,5,
 2,5,3,4,7,0,6,1,0,7,1,6,5,2,4,3,
 5,2,4,3,0,7,1,6,7,0,6,1,2,5,3,4,
 7,0,6,1,2,5,3,4,5,2,4,3,0,7,1,6,
 0,7,1,6,5,2,4,3,2,5,3,4,7,0,6,1,
 6,1,7,0,3,4,2,5,4,3,5,2,1,6,0,7,
 1,6,0,7,4,3,5,2,3,4,2,5,6,1,7,0,
 4,3,5,2,1,6,0,7,6,1,7,0,3,4,2,5,
 3,4,2,5,6,1,7,0,1,6,0,7,4,3,5,2,
 5,2,4,3,0,7,1,6,7,0,6,1,2,5,3,4,
 2,5,3,4,7,0,6,1,0,7,1,6,5,2,4,3
};



/*  ***********************************************************

Encode one bit to tribit

***********************************************************  */

unsigned char EncodeBit(unsigned char bit,unsigned char *trellis)
{
// Put bit into trellis

*trellis = ((*trellis)<<1)|bit;

// Encode trellis state to corresponding tribit

return trellis_to_tribit[*trellis];
}



/*  ***********************************************************

Initialize the paths and the path costs

************************************************************ */

void InitDecoderStatus(DECODER_STATUS *ds)
{
unsigned int i,j;


for(i=0;i<256;i++)
 {
 ds->cost[i]=0.0;
 for(j=0;j<(WINDOW_BYTES-1);j++) ds->path[i][j]=0;
 ds->path[i][WINDOW_BYTES-1] = i;
 }
}

/*  *****************************************************************

Soft Decision Viterbi Decoding

Input - the array of three float values (normaly, -1.0 corresponds to 0, +1.0 to 1)
Output - decoded bit

The delay between the input and the output  bitstreams is 8*WINDOW_BYTES - 1

**************************************************************** */


unsigned char DecodeBit(float *triplet,DECODER_STATUS *ds)
{
unsigned char decoded_bit,min_no,tmp;
unsigned int i,j,k;
float min_cost,x,ftmp;
float tribit_cost[8];
NODE node[256];
unsigned char new_path[256][WINDOW_BYTES];
float new_cost[256];


// Init nodes

for(i=0;i<256;i++) node[i].id=0;


// Calculate the current costs of tribits

for(i=0;i<8;i++)
 {
 tmp=i;
 ftmp=0.0;

 for(j=0;j<3;j++)
  {
  x = (tmp&0x04) ? 1.0 : -1.0;
  tmp<<=1;
  x-=triplet[j];
  ftmp+=x*x;  // Minimum Square Error metrics
  }
 tribit_cost[i]=ftmp;
 }

// Check the all paths

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

 // Shift the path

 for(j=0;j<(WINDOW_BYTES-1);j++)
  {
  ds->path[i][j]=(ds->path[i][j]<<1)|((ds->path[i][j+1]&0x80) ? 1:0);
  }

 ds->path[i][WINDOW_BYTES-1]<<=1;

 // Find the nodes in continuation of the path

 for(j=0;j<2;j++)
  {
  // Trellis continuation

  tmp =(ds->path[i][WINDOW_BYTES-1])|j;

  // Cost of this continuation

  ftmp = tribit_cost[trellis_to_tribit[tmp]] + ds->cost[i];

  // Index

  k=node[tmp].id;

  // Remember continuation parameters

  node[tmp].path_no[k]  = i;
  node[tmp].path_continuation[k] = j;
  node[tmp].path_cost[k] = ftmp;
  node[tmp].id++;
  }
 }

 // Check each node for the path with minimum cost
 // Save the survived paths and costs
 // Find the minimum cost

min_cost = node[0].path_cost[0];
min_no = 0;

for(i=0;i<256;i++)
 {
 k = (node[i].path_cost[0] < node[i].path_cost[1]) ? 0 : 1;


 for(j=0;j<WINDOW_BYTES;j++)
  {
  new_path[i][j]=ds->path[node[i].path_no[k]][j];
  }
 new_path[i][WINDOW_BYTES-1]|=node[i].path_continuation[k];
 new_cost[i]=node[i].path_cost[k];

 if(new_cost[i]<min_cost)
  {
  min_cost=new_cost[i];
  min_no=i;
  }
 }

// Substract the minimum cost
// Store the paths

for(i=0;i<256;i++)
 {
 ds->cost[i]=new_cost[i]-min_cost;

 for(j=0;j<WINDOW_BYTES;j++)
  {
  ds->path[i][j] = new_path[i][j];
  }
 }

decoded_bit = (new_path[min_no][0] & 0x80) ? 1:0;

return decoded_bit;
}


//==========================================================
void main(void)
{
int i,j;
unsigned char x,tmp;
unsigned char trellis = 0;
unsigned char input[1024];
unsigned char output[1024];
float triplet[3];

float frnd;

DECODER_STATUS ds;


InitDecoderStatus(&ds);

for(i=0;i<1024;i++)
 {
 input[i]= x = rand()&1;

 x=EncodeBit(x,&trellis);

 tmp=1;
 for(j=0;j<3;j++)
  {
  frnd = ((float)rand())/((float)RAND_MAX);
  if(frnd<0.15) x^=tmp;
  tmp<<=1;
  }

 for(j=0;j<3;j++)
  {
  triplet[j]=(x&0x04) ? 1.0 : -1.0;
  x<<=1;
  }
 output[i]=DecodeBit(triplet,&ds);
 }

for(i=0;i<1024;i+=32)
 {
 printf("\n\n");
 for(j=0;j<32;j++) printf("%d ",input[i+j]);
 printf("\n");
 for(j=0;j<32;j++) printf("%d ",output[i+j+63]);
 }
}