turbo.cpp

Turbo Code simulator

unsigned decode (double *mesg, double *parity1, double *parity2, unsigned size, bool *boolmesg)
{
	static double **a[2];
   static double **b[2];
   static double *L[2];
   static double **gamma[2][2];
   static double **gammaE[2][2];
   static bool initialized=false;
   static unsigned lastsize=0;
   unsigned returnvalue=ITERATIONS;

   // minimize new's and delete's to only when needed
   if (size != lastsize && initialized == true)
   {
	   // delete all the arrays and rebuild
	   for (int y=0;y<2;y++)
	   	for (int x=0;x<2;x++)
	      {
	         for (int z=0;z<lastsize;z++)
	         {
	  	      	delete gamma[y][x][z];
	  	      	delete gammaE[y][x][z];
	         }

	      	delete gamma[y][x];
	      	delete gammaE[y][x];
	  		}

	   // create L[encoder #]
	   for (int y=0;y<2;y++)
	      delete L[y];

	   // create alpha[encoder #][k][state]
	   for (int x=0;x<2;x++)
	   {
	   	for (int y=0;y<lastsize;y++)
	      {
	         delete a[x][y];
	         delete b[x][y];
		   }

	      delete a[x];
	      delete b[x];
		}
   }

   if (initialized == false || size != lastsize)
   {
   	initialized = true;
      lastsize = size;

	   // create the arrays dynamically at runtime, delete at end of routine
	   // create gamma[encoder #][uk][k][state]
	   for (int y=0;y<2;y++)
	   	for (int x=0;x<2;x++)
	      {
	      	gamma[y][x] = new double*[size];
	         gammaE[y][x] = new double*[size];
	         for (int z=0;z<size;z++)
	         {
	  	      	gamma[y][x][z] = new double[numstates];
	            gammaE[y][x][z] = new double[numstates];
	         }
	  		}

	   // create L[encoder #]
	   for (int y=0;y<2;y++)
	      L[y] = new double[size];

	   // create alpha[encoder #][k][state]
	   for (int x=0;x<2;x++)
	   {
	      a[x] = new double*[size];
	      b[x] = new double*[size];
	      // each Yk has 'numstates' values of gamma
	   	for (int y=0;y<size;y++)
	      {
	         a[x][y] = new double[numstates];
	         b[x][y] = new double[numstates];
		   }
		}
   }

   // initialization of iteration arrays
   for (int x=0;x<numstates;x++)
   {
   	a[0][0][x] = b[0][size-1][x] = a[1][0][x] = (x==0) ? 1.0 : 0.0;
      // extrinsic information from 2-1
   }

	// initialization of extrinsic information array from decoder 2, used in decoder 1
   for (int x=0;x<size;x++)
   	L[1][x] = 0.0;

	// 4*Eb/No
   double Lc = (4.0*1.0)/No;

   for (int c=0;c<ITERATIONS;c++)
   {
		// k from 0 to N-1 instead of 1 to N
	   for (int k=0;k<size;k++)
	   {
	   	// calculate the gamma(s',s);
	   	for (int input=0;input<2;input++)
			{
				double uk = (input == 0) ? -1.0 : 1.0;

	      	for (int s=0;s<numstates;s++)
	         {
	         	double xk = (output[input][s] == 0) ? -1.0 : 1.0;

		      	gammaE[0][input][k][s]=exp(0.5*Lc*parity1[k]*xk);
		      	gamma[0][input][k][s]=exp(0.5*uk*(L[1][k]+Lc*mesg[k]))*gammaE[0][input][k][s];
	         }
	      }
	   }

	   // calculate the alpha terms
	   // from 1 to N-1, 0 is precalculated, N is never used
		for (int k=1;k<size;k++)
	   {
			double temp=0;

	   	// calculate denominator
	      for (int state=0;state<numstates;state++)
		     	temp += a[0][k-1][fromstate[0][state]]*gamma[0][0][k-1][fromstate[0][state]] + a[0][k-1][fromstate[1][state]]*gamma[0][1][k-1][fromstate[1][state]];

	     	for (int state=0;state<numstates;state++)
				a[0][k][state] = (a[0][k-1][fromstate[0][state]]*gamma[0][0][k-1][fromstate[0][state]] + a[0][k-1][fromstate[1][state]]*gamma[0][1][k-1][fromstate[1][state]])/temp;
		}

		// from N-1 to
		for (int k=size-1;k>=1;k--)
	   {
	   	double temp=0;

	   	// calculate denominator
	      for (int state=0;state<numstates;state++)
		     	temp += a[0][k][fromstate[0][state]]*gamma[0][0][k][fromstate[0][state]] + a[0][k][fromstate[1][state]]*gamma[0][1][k][fromstate[1][state]];

	   	for (int state=0;state<numstates;state++)
	      	b[0][k-1][state] = (b[0][k][tostate[0][state]]*gamma[0][0][k][state] + b[0][k][tostate[1][state]]*gamma[0][1][k][state])/temp;
	   }

	   for (int k=0;k<size;k++)
	   {
	   	double min=0;

	      // find the minimum product of alpha, gamma, beta
	      for (int u=0;u<2;u++)
	      	for (int state=0;state<numstates;state++)
	         {
	         	double temp=a[0][k][state]*gammaE[0][u][k][state]*b[0][k][tostate[u][state]];

	            if ((temp < min && temp != 0)|| min == 0)
	            	min = temp;
	         }

	      // if all else fails, make min real small
	      if (min == 0 || min > 1)
	      	min = 1E-100;

	      double topbottom[2];

		   for (int u=0;u<2;u++)
			{
		   	topbottom[u]=0.0;

	   	   for(int state=0;state<numstates;state++)
	         	topbottom[u] += (a[0][k][state]*gammaE[0][u][k][state]*b[0][k][tostate[u][state]]);
			}

	      if (topbottom[0] == 0)
	      	topbottom[0] = min;
	      else if (topbottom[1] == 0)
	      	topbottom[1] = min;

	      L[0][k] = (log(topbottom[1]/topbottom[0]));
	   }

	   interleavedouble(L[0],size);
	   // remember to deinterleave for next iteration
	   interleavedouble(mesg,size);

	   // start decoder 2
	   // code almost same as decoder 1, could combine code into one but too lazy
	   for (int k=0;k<size;k++)
	   {
	   	// calculate the gamma(s',s);
	   	for (int input=0;input<2;input++)
			{
				double uk = (input == 0) ? -1.0 : 1.0;

	      	for (int s=0;s<numstates;s++)
	         {
	         	double xk = (output[input][s] == 0) ? -1.0 : 1.0;

		      	gammaE[1][input][k][s]=exp(0.5*Lc*parity2[k]*xk);
		      	gamma[1][input][k][s]=exp(0.5*uk*(L[0][k]+Lc*mesg[k]))*gammaE[1][input][k][s];
	         }
	      }
	   }

	   // calculate the alpha terms
		for (int k=1;k<size;k++)
	   {
			double temp=0;

	   	// calculate denominator
	      for (int state=0;state<numstates;state++)
		     	temp += a[1][k-1][fromstate[0][state]]*gamma[1][0][k-1][fromstate[0][state]] + a[1][k-1][fromstate[1][state]]*gamma[1][1][k-1][fromstate[1][state]];

	     	for (int state=0;state<numstates;state++)
				a[1][k][state] = (a[1][k-1][fromstate[0][state]]*gamma[1][0][k-1][fromstate[0][state]] + a[1][k-1][fromstate[1][state]]*gamma[1][1][k-1][fromstate[1][state]])/temp;
		}

	   // in the first iteration, set b[1][N-1] = a[1][N-1] for decoder 2
	   // this decoder can't be terminated to state 0 because of the interleaver
	   // the performance loss is supposedly neglible.
      if (c==0)
	   {
			double temp=0;

		  	// calculate denominator
		   for (int state=0;state<numstates;state++)
			  	temp += a[1][size-1][fromstate[0][state]]*gamma[1][0][size-1][fromstate[0][state]] + a[1][size-1][fromstate[1][state]]*gamma[1][1][size-1][fromstate[1][state]];

		  	for (int state=0;state<numstates;state++)
				b[1][size-1][state] = (a[1][size-1][fromstate[0][state]]*gamma[1][0][size-1][fromstate[0][state]] + a[1][size-1][fromstate[1][state]]*gamma[1][1][size-1][fromstate[1][state]])/temp;
	   }

		for (int k=size-1;k>=1;k--)
	   {
	   	double temp=0;

	   	// calculate denominator
	      for (int state=0;state<numstates;state++)
		     	temp += a[1][k][fromstate[0][state]]*gamma[1][0][k][fromstate[0][state]] + a[1][k][fromstate[1][state]]*gamma[1][1][k][fromstate[1][state]];

	   	for (int state=0;state<numstates;state++)
	      	b[1][k-1][state] = (b[1][k][tostate[0][state]]*gamma[1][0][k][state] + b[1][k][tostate[1][state]]*gamma[1][1][k][state])/temp;
	   }

	   for (int k=0;k<size;k++)
	   {
	   	double min = 0;

	      // find the minimum product of alpha, gamma, beta
	      for (int u=0;u<2;u++)
	      	for (int state=0;state<numstates;state++)
	         {
	         	double temp=a[1][k][state]*gammaE[1][u][k][state]*b[1][k][tostate[u][state]];

	            if ((temp < min && temp != 0)|| min == 0)
	            	min = temp;
	         }
	      // if all else fails, make min real small
	      if (min == 0 || min > 1)
	      	min = 1E-100;

			double topbottom[2];

		   for (int u=0;u<2;u++)
			{
		   	topbottom[u]=0.0;

	   	   for(int state=0;state<numstates;state++)
					topbottom[u] += (a[1][k][state]*gammaE[1][u][k][state]*b[1][k][tostate[u][state]]);
			}

	      if (topbottom[0] == 0)
	      	topbottom[0] = min;
	      else if (topbottom[1] == 0)
	      	topbottom[1] = min;

	      L[1][k] = (log(topbottom[1]/topbottom[0]));
	   }

		deinterleavedouble(mesg,size);
	   deinterleavedouble(L[1],size);
		// get L[0] back to normal after decoder 2
		deinterleavedouble(L[0],size);

      bool temp=true;
      for (int k=0;k<size;k++)
			if (boolmesg[k] != ((Lc*mesg[k] + L[0][k] + L[1][k] > 0.0) ? true : false))
				temp = false;

      // we can quit prematurely since it has been decoded
      if (temp==true)
      {
      	returnvalue = c;
			c=ITERATIONS;
      }
	}
	// end decoder 2

   // make decisions
   for (int k=0;k<size;k++)
		if ((Lc*mesg[k] + L[0][k] + L[1][k]) > 0)
      	mesg[k] = 1.0;
      else
      	mesg[k] = -1.0;

      return returnvalue;
}

void deinterleavedouble(double *mesg, unsigned size)
{
	double *temp;

   temp = new double[size];

   for (int x=0;x<size;x++)
   	temp[x] = mesg[x];

	for (int x=0;x<size;x++)
		mesg[deinterleavearray[x]] = temp[x];

   delete temp;
}

void interleavedouble(double *mesg, unsigned size)
{
	double *temp;

   temp = new double[size];

   for (int x=0;x<size;x++)
   	temp[x] = mesg[x];

   for (int x=0;x<size;x++)
   	mesg[interleavearray[x]] = temp[x];

   delete temp;
}

void interleave(bool *mesg,unsigned size)
{
	bool *temp;

   temp = new bool[size];

   for (int x=0;x<size;x++)
   	temp[x] = mesg[x];

   for (int x=0;x<size;x++)
   	mesg[interleavearray[x]] = temp[x];

   delete temp;
}

void deinterleave(bool *mesg,unsigned size)
{
	bool *temp;

   temp = new bool[size];

   for (int x=0;x<size;x++)
   	temp[x] = mesg[x];

	for (int x=0;x<size;x++)
		mesg[deinterleavearray[x]] = temp[x];

   delete temp;
}

void createinterleave(unsigned size)
{
   bool *yesno;

   yesno = new bool[size];

   for (int x=0;x<N;x++)
   	yesno[x]=false;

	// create an interleave array
   for (int x=0;x<N;x++)
   {
   	unsigned  val;

      do
      {
      	val=before.longrandom(N);
      }
		while(yesno[val] == true);

      yesno[val] = true;
		interleavearray[x] = val;
      deinterleavearray[val] = x;
	}

   delete yesno;
}

void encode(bool *mesg,bool *parity,unsigned size, bool force)
{

	unsigned state=0;

  	for (int x=0;x<size;x++)
	{
		// force the encoder to zero state at the end
     	if (x>=size-memory && force)
      {
			if (tostate[0][state]&1)
         	mesg[x] = true;
         else
         	mesg[x] = false;
      }

		// can't assume the bool type has an intrinsic value of 0 or 1
		// may differ from platform to platform
      int uk = mesg[x] ? 1 : 0;

		// calculate output due to new mesg bit
      parity[x] = output[uk][state];
      // calculate the new state
		state = tostate[uk][state];
	}
}

bool add(bool a, bool b)
{
   return a==b ? false : true;
}

double gaussian(double variance)
{
	// static becuase we don't want to have it initialized each time we go in
   double returnvalue=0;
   double k;

	k = sqrt(variance/2.0);

   // add 24 uniform RV to obtain a simulation of normality
   for (int x=0;x<24;x++)
		returnvalue += before.doublerandom();

	return k*(returnvalue-0.5*24);
}