turbo_map_punc.cpp

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

// -----------------------------------------------------------------------
// File: turbo_punc.cpp
// Date: May 5, 2002
// Authors: M. Walma and R.H. Morelos-Zaragoza
//
// Descirption: Binary turbo code with puncturing matrix for deriving a 
//              rate-1/2 code
//
// Based on a program written by Mathys Walma, who placed it in the public 
// domain in 1998. The user interface was modified by R. Morelos-Zaragoza.
// ------------------------------------------------------------------------
// This program is complementary material for the book:
//
// R.H. Morelos-Zaragoza, The Art of Error Correcting Coding, Wiley, 2002.
//
// ISBN 0471 49581 6
//
// This and other programs are available at http://the-art-of-ecc.com
//
// You may use this program for academic and personal purposes only. 
// If this program is used to perform simulations whose results are 
// published in a journal or book, please refer to the book above.
//
// The use of this program in a commercial product requires explicit
// written permission from the author. The author is not responsible or 
// liable for damage or loss that may be caused by the use of this program. 
//
// Copyright (c) 2002. Robert H. Morelos-Zaragoza. All rights reserved.
// ------------------------------------------------------------------------
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include "random.h"

#define MAXRAND       LONG_MAX

int          states;
long         seed;
unsigned int N;
unsigned int ITERATIONS;
long         maxsim;
double       No;
double       sigma;
int          temp;


double gaussian(double variance);
void encode(bool *mesg,bool *parity1,unsigned size, bool force);
void interleave(bool *mesg, unsigned size);
bool add(bool a,bool b);
void interleave(bool *mesg, unsigned size);
void deinterleave(bool *mesg, unsigned size);
void deinterleavedouble(double *mesg, unsigned size);
void interleavedouble(double *mesg, unsigned size);
void createinterleave(unsigned size);
void createencodetable(int states);
void inttobool(unsigned state, bool *array, unsigned size);
void booltoint(bool *array, unsigned size, unsigned *state);
unsigned decode (double *channelmesg, double *parity1, double *parity2, unsigned size, bool *mesg);

unsigned *interleavearray;
unsigned *deinterleavearray;

unsigned numstates;
unsigned memory;

// [2] = input, [16] = current state, tostate[2][16] = next state
unsigned *tostate[2];

// [2] = last input, [16] = current state, fromstate[2][16] = previous state
unsigned *fromstate[2];

// [2] = input, [16] = current state, output[2][16] = output of encoder
bool *output[2];

Random before;


   // Puncturing matrix
   int punc[3][2] ={ 1, 1, 
                     1, 0, 
                     0, 1 };


main(int argc, char *argv[])
{
   bool *mesg;
   bool *parity1;
   bool *parity2;
   double *channelmesg;
   double *channelparity1;
   double *channelparity2;
   register int i;
   char name1[40], name2[40];
   FILE *fp;

   // Command line processing
   if (argc != 10)
     {
       printf("Usage: %s Num_states Max_iter INT_file", argv[0]);
       printf(" BER_file EbNo_init EbNo_final EbNo_inc  max_sims  seed\n");
       exit(0);
     }

   double INIT, FINAL, INC;

   sscanf(argv[1],"%d",  &states);
   sscanf(argv[2],"%d",  &ITERATIONS);
   sscanf(argv[3],"%s",  name1);
   sscanf(argv[4],"%s",  name2);
   sscanf(argv[5],"%lf", &INIT);
   sscanf(argv[6],"%lf", &FINAL);
   sscanf(argv[7],"%lf", &INC);
   sscanf(argv[8],"%ld", &maxsim);
   sscanf(argv[9],"%ld", &seed);



   // READ INTERLEAVER
   fp = fopen(name1,"r");
   N = 0;
   while (!(fscanf(fp,"%d", &temp)==EOF))
     N++;
   fclose(fp);
   interleavearray = new unsigned[N];
   deinterleavearray = new unsigned[N];
   fp = fopen(name1,"r");
   i = 0;
   while (!(fscanf(fp,"%d", &interleavearray[i])==EOF))
     i++;
   fclose(fp);
   for (i=0; i<N; i++) deinterleavearray[interleavearray[i]]=i;

   double rate = 0.5;

   printf("\nParallel concatenated (turbo) code of rate = %lf\n", rate);
   printf("%2d-state constituent recursive convolutional codes \n\n", states);

   srandom(seed);

   mesg = new bool[N];
   parity1 = new bool[N];
   parity2 = new bool[N];
   channelmesg = new double[N];
   channelparity1 = new double[N];
   channelparity2 = new double[N];


   fp = fopen(name2,"w");

   createencodetable(states);

   rate = (double) N / ( 2.0*(double)N + log(states)/log(2.0) );


   for (double dB=INIT; dB<FINAL+1.0e-1; dB += INC)
   {
	unsigned totalN=0;
	unsigned totalerrors=0;
	unsigned numiter;
	unsigned totaliter=0;
	unsigned numloops=0;

      	// dB = 10*log(Eb/No) where Eb is 1
      	// No = 1.0/pow(10.0,dB/10.0);
      	sigma = 1.0/sqrt(2.0*0.3333333*pow(10.0,dB/10.0));
	
      	bool keepgoing=true;

	do {
		for (int x=0;x<N;x++)
		mesg[x] = before.boolrandom();

         	// create a random interleaver for each decode trial
		// createinterleave(N);

		encode(mesg,parity1,N,true);
		interleave(mesg,N);
		encode(mesg,parity2,N,false);
		deinterleave(mesg,N);

		for (int x=0;x<N;x++)
		   {
		   mesg[x]    ?    channelmesg[x]=1.0 :    channelmesg[x]=-1.0;
		   parity1[x] ? channelparity1[x]=1.0 : channelparity1[x]=-1.0;
		   parity2[x] ? channelparity2[x]=1.0 : channelparity2[x]=-1.0;
		   }

		// add gaussian noise, mean=0, variance=No/2
		for (int x=0;x<N;x++)
		   {
		    	//channelmesg[x] += gaussian(No/2);
			//channelparity1[x] += gaussian(No/2);
		      	//channelparity2[x] += gaussian(No/2);
		    	channelmesg[x] += gaussian(sigma);
			channelparity1[x] += gaussian(sigma);
		      	channelparity2[x] += gaussian(sigma);
		   }

		numiter = 
                    decode(channelmesg, channelparity1, channelparity2,N,mesg);

		unsigned numerrors=0;

		for (int x=0;x<N;x++)
		   {
		   bool temp = channelmesg[x] == 1 ? true : false;

		   if (mesg[x] != temp)
		   	numerrors++;
		   }

	      	totalerrors += numerrors;
	      	totalN += N;
         	totaliter += numiter;
         	numloops++;

	      	char name[30];

		// stopping criterion
		if (totalerrors> 1000)
         		keepgoing = false;
		else if (totalN>maxsim)
         		keepgoing=false;
	}
	while (keepgoing == true); 

	printf ("%5.2lf  %10.7e  %12ld %12ld     %lf     %d\n",  dB,
                (totalerrors/(double)totalN), totalerrors, totalN, 
                           (double)totaliter/(double)numloops, N);

        fprintf(fp, "%8.4lf %10.7e\n", dB, (totalerrors/(double)totalN));
        fflush(fp);
   }

   //delete interleavearray;
   //delete deinterleavearray;
   //delete mesg;
   //delete parity1;
   //delete parity2;
   //delete channelmesg;
   //delete channelparity1;
   //delete channelparity2;
}




void createencodetable(int states)
{
   bool *boolstate;
   bool *newstate;

if (states == 4) {         // Generators {7,5}
   numstates = 4;
   memory = 2;

   // create arrays used by encode and decode
   output[0] = new bool[numstates];
   output[1] = new bool[numstates];
   fromstate[0] = new unsigned[numstates];
   fromstate[1] = new unsigned[numstates];
   tostate[0] = new unsigned[numstates];
   tostate[1] = new unsigned[numstates];

   boolstate = new bool[memory];
   newstate = new bool[memory];

   for (unsigned input=0;input<2;input++)
      for (unsigned intstate=0;intstate<numstates;intstate++)
      {
         bool boolinput = (input == 0) ? false : true;

         inttobool(intstate,boolstate,memory);

         // calculate output due to the input
         output[input][intstate] = add(boolinput,boolstate[0]);
         output[input][intstate] = add(output[input][intstate],boolstate[1]);

         // calculate new states
         newstate[1] = boolstate[0];
         newstate[0] = add(add(boolinput,boolstate[0]),boolstate[1]);
         // from s' to s
         booltoint (newstate,memory,&tostate[input][intstate]);
         // from s to s'
         fromstate[input][tostate[input][intstate]] = intstate;
      }

}

else if (states == 8) {         // Generators {17,15}
   numstates = 8;
   memory = 3;

   // create arrays used by encode and decode
   output[0] = new bool[numstates];
   output[1] = new bool[numstates];
   fromstate[0] = new unsigned[numstates];
   fromstate[1] = new unsigned[numstates];
   tostate[0] = new unsigned[numstates];
   tostate[1] = new unsigned[numstates];

   boolstate = new bool[memory];
   newstate = new bool[memory];

   for (unsigned input=0;input<2;input++)
      for (unsigned intstate=0;intstate<numstates;intstate++)
      {
         bool boolinput = (input == 0) ? false : true;

         inttobool(intstate,boolstate,memory);

         // calculate output due to the input
         //output[input][intstate] = add(boolinput,boolstate[0]);
         output[input][intstate] = add(boolinput,boolstate[1]);
         output[input][intstate] = add(output[input][intstate],boolstate[2]);
         output[input][intstate] = add(output[input][intstate],boolstate[2]);

         // calculate new states
         newstate[2] = boolstate[1];
         newstate[1] = boolstate[0];
         newstate[0] = add(add(add(boolinput,boolstate[0]),boolstate[1]),
                                                               boolstate[2]);
         // from s' to s
         booltoint (newstate,memory,&tostate[input][intstate]);
         // from s to s'
         fromstate[input][tostate[input][intstate]] = intstate;
      }

}

else if (states == 16) {         // Generators {31,27}
   numstates = 16;
   memory = 4;

   // create arrays used by encode and decode
   output[0] = new bool[numstates];
   output[1] = new bool[numstates];
   fromstate[0] = new unsigned[numstates];
   fromstate[1] = new unsigned[numstates];
   tostate[0] = new unsigned[numstates];
   tostate[1] = new unsigned[numstates];

   boolstate = new bool[memory];
   newstate = new bool[memory];

   for (unsigned input=0;input<2;input++)
      for (unsigned intstate=0;intstate<numstates;intstate++)
      {
         bool boolinput = (input == 0) ? false : true;

         inttobool(intstate,boolstate,memory);

         // calculate output due to the input
         output[input][intstate] = add(boolinput,boolstate[0]);
         output[input][intstate] = add(output[input][intstate],boolstate[3]);
         output[input][intstate] = add(output[input][intstate],boolstate[1]);
         output[input][intstate] = add(output[input][intstate],boolstate[2]);
         output[input][intstate] = add(output[input][intstate],boolstate[3]);

         // calculate new states
         newstate[3] = boolstate[2];
         newstate[2] = boolstate[1];
         newstate[1] = boolstate[0];
         newstate[0] = add(add(boolinput,boolstate[0]),boolstate[3]);
         // from s' to s
         booltoint (newstate,memory,&tostate[input][intstate]);
         // from s to s'
         fromstate[input][tostate[input][intstate]] = intstate;
      }

}

   delete boolstate;
   delete newstate;
}





void inttobool(unsigned state, bool *array, unsigned size)
{
	for (unsigned x=0;x<size;x++)
   {
   	unsigned next = state >> 1;

		if ((next << 1) == state)
      	array[x] = false;
      else
      	array[x] = true;

      state = next;
   }
}




void booltoint(bool *array, unsigned size, unsigned *state)
{
	*state = 0;

   for (int x=0;x<size;x++)
   	if (array[x] == true)
      	(*state) |= (1 << x);
}





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