turbo_map_punc.cpp

Turbo码BCJR传统经典算法的C语言实现,经典,对于学习MAP算法的人很有帮助.

	// 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;

   // 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/No;
   double Lc = 2.0/sigma;

   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;

                 if (punc[k%2][1])
		   gammaE[0][input][k][s]=exp(0.5*Lc*parity1[k]*xk);
                 else
		   gammaE[0][input][k][s]=1.0;
                 if (punc[k%2][0])
		   gamma[0][input][k][s]
                      = exp(0.5*uk*(L[1][k]+Lc*mesg[k]))*gammaE[0][input][k][s];
                 else
		   gamma[0][input][k][s]
                      = exp(0.5*uk*L[1][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;
	   }

	// calculate the alpha terms
	// from N-1 to 1
	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;
	   }

	// compute the extrinsic LLRs
	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;

                 if (punc[k%2][2])
		    gammaE[1][input][k][s] = exp(0.5*Lc*parity2[k]*xk);
                 else
		    gammaE[1][input][k][s] = 1.0;
                 if (punc[k%2][0])
		    gamma[1][input][k][s] = exp(0.5*uk*(L[0][k]+Lc*mesg[k]))
                                                     * gammaE[1][input][k][s];
                 else
		    gamma[1][input][k][s] = exp(0.5*uk*L[0][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+1;
		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 sigma)
{
  double rndm, u1, u2, s, noise;
    do {
        rndm = (double)(random())/MAXRAND;
        u1 = rndm * 2 - 1.0;
        rndm = (double)(random())/MAXRAND;
        u2 = rndm * 2 - 1.0;
        s = u1 * u1 + u2 * u2;
      } while( s >= 1 );
    noise = u1 * sqrt( (-2.0*log(s))/s );
  return(noise*sigma);
}