turbo.h

turbo码的相关程序

#define CONSTRAINT 5			// memory length of the component convolution code 
#define POLY1 031				// polynomial of encoder1
#define POLY2 033				// polynomial of encoder2
#define BLOCKSIZE 1024			// the length of information blocks
#define ITERATION 15			// maximum iteration numbers

// puncture or not
#define PUNCTURE
#ifdef PUNCTURE
	#define CODERATE 1.0/2.0
#else
	#define CODERATE 1.0/3.0
#endif

// look up table or not
#define LOOKUPTABLE 1
#define MaxStar_scale 15.0
#define MaxStar_precision 256
#define MaxStar_trunc(x) ((x>MaxStar_precision-1) ? (MaxStar_precision-1) : x)
#define MaxStar_index(x) ((int)(floor(MaxStar_trunc(x*MaxStar_precision/MaxStar_scale))))

class TURBOCODE
{
public:
	int BlockSize;
	double RATE;
	double EbN0;
	double Noise;		// sigma	
	int * SourceBits, *DecodedBits;	
	int * EncodedBits;
	double * CodedBits;
	int ErrorBits;
	int ErrorBlocks;
	float AverageIterNum;
	int PunctureFlag;
	int LUTFlag;		
	
private:
	int	Constraint;
	int	PartabLen;
	int	State;	
	int G1;
	int G2;
	int TurboIteration;
	int mask;
	double INF;

	unsigned char * Partab;
	int * Interleaver, *Deinterleaver;		
	int * tempEncodedBits;
	int ** nextstate,**previousstate;
	int ** output;
	long looptimes;		
	double **Alpha1,**Beta1,***Gamma1,***GammaE1;
	double **Alpha2,**Beta2,***Gamma2,***GammaE2;	
	double *L12,*L21;	
	double *L1;
	
	double MaxStar_lookup[MaxStar_precision];

public:
	TURBOCODE::TURBOCODE(int blocksize,int constraint, int poly1, int poly2, double rate, int iteration, int punctureflag, int lutflag)
	{
		int i,j;

		BlockSize = blocksize;
		Constraint = constraint;
		PartabLen = (int) pow(2,5);
		State = (int) pow(2,Constraint - 1);		
		RATE = rate;
		G1 = poly1;
		G2 = poly2;
		TurboIteration = 15;
		mask = 0xF;
		INF = -1.0e200;
		
		PunctureFlag = punctureflag;		
		LUTFlag = lutflag;		

		Interleaver = new int [BlockSize + 1];
		Deinterleaver = new int [BlockSize + 1];
		
		SourceBits = new int[BlockSize+1];
		EncodedBits = new int[(int)(BlockSize/RATE)+1];
		tempEncodedBits = new int[BlockSize*3+1];
		DecodedBits = new int[BlockSize+1];
		CodedBits = new double[BlockSize*3+1];

		Partab = new unsigned char[PartabLen];
		nextstate = new int*[State];
		previousstate = new int*[State];
		output = new int*[State];
		for (i=0;i<State;i++)
		{
			nextstate[i] = new int[2];
			previousstate[i] = new int[2];
			output[i] = new int[2];
		}

		GeneratePar();
		GenerateTrellis();
		interleaver();
		deinterleaver();
		
		Alpha1 = new double*[BlockSize+1];
		Beta1 = new double*[BlockSize+1];
		Gamma1 = new double**[BlockSize+1];
		GammaE1 = new double**[BlockSize+1];

		Alpha2 = new double*[BlockSize+1];
		Beta2 = new double*[BlockSize+1];
		Gamma2 = new double**[BlockSize+1];
		GammaE2 = new double**[BlockSize+1];

		for (i=0;i<=BlockSize;i++)
		{
			Alpha1[i] = new double[State];
			Alpha2[i] = new double[State];
			Gamma1[i] = new double*[State];
			GammaE1[i] = new double*[State];
			Beta1[i] = new double[State];
			Beta2[i] = new double[State];
			Gamma2[i] = new double*[State];
			GammaE2[i] = new double*[State];
			for (j=0;j<State;j++)
			{
				Gamma1[i][j] =  new double[2];
				Gamma2[i][j] =  new double[2];
				GammaE1[i][j] =  new double[2];
				GammaE2[i][j] =  new double[2];
			}
		}

		L12 = new double[BlockSize+1];
		L21 = new double[BlockSize+1];
		L1 = new double[BlockSize+1];

		if (LUTFlag)
		{
			for (i=0;i<MaxStar_precision;i++)
				MaxStar_lookup[i] = log(1 + exp(-1.0*i*MaxStar_scale/MaxStar_precision));
		}
	}

	TURBOCODE::~TURBOCODE()
	{
		int i,j;

		delete [] Interleaver;
		delete [] Deinterleaver;
		
		delete [] Partab;

		for (i=0;i<State;i++)
		{
			delete [] nextstate[i];
			delete [] previousstate[i];
			delete [] output[i];
		}
		delete [] nextstate;
		delete [] previousstate;
		delete [] output;

		delete [] SourceBits;
		delete [] EncodedBits;
		delete [] DecodedBits;
		delete [] tempEncodedBits;
		delete [] CodedBits;

		for (i=0;i<=BlockSize;i++)
		{
			delete [] Alpha1[i];
			delete [] Alpha2[i];			
			delete [] Beta1[i];
			delete [] Beta2[i];

			for (j=0;j<State;j++)
			{
				delete Gamma1[i][j];
				delete Gamma2[i][j];
				delete GammaE1[i][j];
				delete GammaE2[i][j];
			}

			delete [] Gamma1[i];
			delete [] Gamma2[i];			
			delete [] GammaE1[i];
			delete [] GammaE2[i];
		}
		delete [] Gamma1;
		delete [] Gamma2;			
		delete [] GammaE1;
		delete [] GammaE2;

		delete [] L12;
		delete [] L21;
		delete [] L1;
	}

	int TURBOCODE::GeneratePar()
	{
		int i,j;
		int bits;
		int parity;		

		for (i=0;i<PartabLen;i++)
		{
			bits = i;
			parity = 0;
			for (j=0;j<Constraint;j++)
			{
				parity = parity ^ (bits & 1);
				bits = bits >> 1;

			}
			Partab[i] = parity;
		}
		return 1;
	}

	int TURBOCODE::GenerateTrellis()
	{
		int i,j;
		int RegState;
		int Feedback;
		int index;
		int parity;
		int RegState1,RegState2;

		RegState1 = 0;
		RegState2 = 0;


		for (i=0;i<State;i++)
		{
			RegState = i;
			for (j=0; j<2; j++)
			{
				RegState = i;
				RegState = RegState | (j<<(Constraint-1));
				index = RegState & G1;
				Feedback = Partab[index];
				RegState = (RegState & mask) | (Feedback<<(Constraint-1));
				index = RegState & G2;
				parity = Partab[index];
				nextstate[i][j] = RegState >> 1;
				previousstate[RegState>>1][j] = i;
				output[RegState>>1][j] = parity;
			}
		}
		return 1;
	}

	int TURBOCODE::interleaver()
	{
		FILE * F_interleaver;
		int i;
		int index;

		if ((F_interleaver = fopen("interleaver.txt","r")) == NULL)
		{
			printf("open interleaver.txt error");
			return -1;
		}

		for (i=1;i<BlockSize+1;i++)
		{
			fscanf(F_interleaver,"%d",&index);
			Interleaver[i] = index;
		}

		fclose(F_interleaver);
		return 1;
	}

	int TURBOCODE::deinterleaver()
	{
		FILE * F_deinterleaver;
		int i;
		int index;

		if ((F_deinterleaver = fopen("interleaver.txt","r")) == NULL)
		{
			printf("open interleaver.txt error");
			return -1;
		}

		for (i=1;i<BlockSize+1;i++)
		{
			fscanf(F_deinterleaver,"%d",&index);
			Deinterleaver[index] = i;
		}

		fclose(F_deinterleaver);
		return 1;
	}
	
	int TURBOCODE::forcezero(int regstate)
	{
		int input,feedback;
		int i;
		int Regstate;

		Regstate = regstate;

		for (i=Constraint -1; i > 0 ; i--)
		{
			Regstate = regstate;
			input = 0;
			Regstate = Regstate | (input << (Constraint-1));
			feedback = Partab[Regstate & G1];

			if (feedback == 0)
			{
				SourceBits[BlockSize + 1 - i] = input;
				Regstate = (Regstate & mask) | (feedback<<(Constraint-1));
				regstate = Regstate >> 1;

			}
			else
			{
				Regstate = regstate;
				input = 1;
				Regstate = Regstate | (input << (Constraint-1));
				feedback = Partab[Regstate & G1];

				if (feedback == 0)
				{
					SourceBits[BlockSize + 1 - i] = input;
					Regstate = (Regstate & mask) | (feedback<<(Constraint-1));
					regstate = Regstate >> 1;
				}
			}
		}
		return 1;
	}

	void TURBOCODE::encoder()
	{
		int RegState1,RegState2;
		int i;
		int input;
		int code1,code2; // code1 is produced by RSC1, code2 is produced by RSC2

		// RSC1
		RegState1 = 0;
		for (i=1; i<BlockSize+1-(Constraint-1) ; i++)
		{
			input = SourceBits[i];
			tempEncodedBits[3*i-2] = input;

			RegState1 = nextstate[RegState1][input];
			code1 = output[RegState1][input];

			tempEncodedBits[3*i-1] = code1;
		}

		forcezero(RegState1);

		for (i=Constraint-1; i>0; i--)
		{
			input = SourceBits[BlockSize+1-i];
			tempEncodedBits[3*(BlockSize+1-i)-2] = input;

			RegState1 = nextstate[RegState1][input];
			code1 = output[RegState1][input];

			tempEncodedBits[3*(BlockSize+1-i)-1] = code1;
		}

		//RSC2
		RegState2 = 0;
		for (i=1; i<BlockSize+1; i++)
		{
			input = SourceBits[Interleaver[i]];  // Interleaver
			RegState2 = nextstate[RegState2][input];
			code2 = output[RegState2][input];
			tempEncodedBits[3*i] =  code2;
		}		

		//puncturing
		if (PunctureFlag)
		{
			int k = 0;
			for (i=1;i<BlockSize+1;i++)
			{
				EncodedBits[i*2-1] = tempEncodedBits[i*3-2] > 0 ? 1 : -1;
				k++;
				if (k%2 != 0)			
					EncodedBits[i*2] = tempEncodedBits[i*3-1] > 0 ? 1 : -1;
				else			
					EncodedBits[i*2] = tempEncodedBits[i*3] > 0 ? 1 : -1;
			}
		}
		else
			for (i=1;i<BlockSize+1;i++)
			{
					EncodedBits[i*3-2] = tempEncodedBits[i*3-2] > 0 ? 1 : -1;
					EncodedBits[i*3-1] = tempEncodedBits[i*3-1] > 0 ? 1 : -1;
					EncodedBits[i*3] = tempEncodedBits[i*3] > 0 ? 1 : -1;
			}
	}


	void TURBOCODE::log_BCJR_decoder()
	{
		int i,j,k;
		int times;		
		double L12_d,L21_n,L21_d,L12_n;
		double Lc;				
		double x,y,tempL,absx_y;
		int ErrorFlag;

		double y1,y1p,y2,y2p;			

		//Initialization
		Lc = 2 / (Noise * Noise);
		Alpha1[0][0]=0;
		Alpha2[0][0]=0;

		for (i=1;i<State;i++)
		{
			Alpha1[0][i] = INF;
			Alpha2[0][i] = INF;
		}

		for (i=1;i<BlockSize+1;i++)
		{
			L21[i] = 0;
		}

		Beta1[BlockSize][0] = 0;

		for (i=1; i<State; i++)
		{
			Beta1[BlockSize][i] = INF;
		}


		ErrorFlag = 1;
		for (times=0;(times<TurboIteration) && (ErrorFlag == 1);times++)
		{

			for (k=1; k<BlockSize+1; k++)
			{
				y1 = CodedBits[3*k-2];
				y1p = CodedBits[3*k-1];

				for (i=0;i<State;i++)
					for (j=0;j<2;j++)
						Gamma1[k][i][j] = (j?1.0:-1.0)*Lc*y1 + 0.5*(j?1.0:-1.0) * L21[Deinterleaver[k]]
											+ (output[i][j]?1.0:-1.0)*y1p*Lc;						
			}

			for (k=1; k<BlockSize+1; k++)
			{
				for (i=0;i<State;i++)
				{
					x = Alpha1[k-1][previousstate[i][0]] + Gamma1[k][i][0];
					y = Alpha1[k-1][previousstate[i][1]] + Gamma1[k][i][1];
					absx_y = (x>y)?(x-y):(y-x);
					if (LUTFlag)						
						Alpha1[k][i] = (__max(x,y) + MaxStar_lookup[MaxStar_index(absx_y)]);
					else
						Alpha1[k][i] = (__max(x,y) + log(1 + exp(-1.0 * absx_y)));
				}
			}

			for (k=BlockSize;k>1;k--)
			{

				for (i=0;i<State;i++)
				{
					x = Beta1[k][nextstate[i][0]] + Gamma1[k][nextstate[i][0]][0];
					y = Beta1[k][nextstate[i][1]] + Gamma1[k][nextstate[i][1]][1];
					absx_y = (x>y)?(x-y):(y-x);					
					if (LUTFlag)
						Beta1[k-1][i] = (__max(x,y) + MaxStar_lookup[MaxStar_index(absx_y)]);
					else
						Beta1[k-1][i] = (__max(x,y) + log(1 + exp(-1.0 * absx_y)));
				}
			}

			for (k=1;k<BlockSize+1;k++ )
			{
				L12_n = INF;
				L12_d = INF;
				y1p =  CodedBits[3*k-1];

				for (i=0;i<State;i++)
				{
					y =Alpha1[k-1][previousstate[i][1]] + y1p*(output[i][1]?1.0:-1.0)*Lc + Beta1[k][i];					
					tempL = __max(L12_n,y);
					absx_y = (L12_n>y)?(L12_n-y):(y-L12_n);					
					if (LUTFlag)
						L12_n = tempL + MaxStar_lookup[MaxStar_index(absx_y)];
					else
						L12_n = tempL + log(1 + exp(-1.0 * absx_y));

					y = Alpha1[k-1][previousstate[i][0]] + y1p*(output[i][0]?1.0:-1.0)*Lc + Beta1[k][i];					
					tempL = __max(L12_d,y);
					absx_y = (L12_d>y)?(L12_d-y):(y-L12_d);
					if (LUTFlag)
						L12_d = tempL + MaxStar_lookup[MaxStar_index(absx_y)];
					else
						L12_d = tempL + log(1 + exp(-1.0 * absx_y));
				}
				L12[k] = L12_n - L12_d;
			}

			//  D2
			for (k=1; k<BlockSize+1; k++)
			{
				y2 = CodedBits[Interleaver[k] * 3 - 2];
				y2p = CodedBits[k * 3];

				for (i=0;i<State;i++)
					for (j=0;j<2;j++)
						Gamma2[k][i][j] = (j?1.0:-1.0)*Lc*y2 + 0.5*(j?1.0:-1.0) * L12[Interleaver[k]]
											+ (output[i][j]?1.0:-1.0)*y2p*Lc;						
			}

			for (k=1; k<BlockSize+1; k++)
			{
				for (i=0;i<State;i++)
				{
					x = Alpha2[k-1][previousstate[i][0]] + Gamma2[k][i][0];
					y = Alpha2[k-1][previousstate[i][1]] + Gamma2[k][i][1];
					absx_y = (x>y)?(x-y):(y-x);
					if (LUTFlag)						
						Alpha2[k][i] = (__max(x,y) + MaxStar_lookup[MaxStar_index(absx_y)]);
					else
						Alpha2[k][i] = (__max(x,y) + log(1 + exp(-1.0 * absx_y)));
				}
			}

			for (i=0;i<State;i++)
				Beta2[BlockSize][i] = Alpha2[BlockSize][i];

			for (k=BlockSize;k>1;k--)
			{
				for (i=0;i<State;i++)
				{
					x = Beta2[k][nextstate[i][0]] + Gamma2[k][nextstate[i][0]][0];
					y = Beta2[k][nextstate[i][1]] + Gamma2[k][nextstate[i][1]][1];
					absx_y = (x>y)?(x-y):(y-x);
					if (LUTFlag)
						Beta2[k-1][i] = (__max(x,y) + MaxStar_lookup[MaxStar_index(absx_y)]);
					else
						Beta2[k-1][i] = (__max(x,y) + log(1 + exp(-1.0 * absx_y)));
				}
			}

			for (k=1;k<BlockSize+1;k++ )
			{
				L21_n = INF;
				L21_d = INF;
				y2p = CodedBits[k * 3];

				for (i=0;i<State;i++)
				{
					y =Alpha2[k-1][previousstate[i][1]] + y2p*(output[i][1]?1.0:-1.0)*Lc + Beta2[k][i];					
					tempL = __max(L21_n,y);
					absx_y = (L21_n>y)?(L21_n-y):(y-L21_n);					
					if (LUTFlag)
						L21_n = tempL + MaxStar_lookup[MaxStar_index(absx_y)];
					else
						L21_n = tempL + log(1 + exp(-1.0 * absx_y));

					y = Alpha2[k-1][previousstate[i][0]] + y2p*(output[i][0]?1.0:-1.0)*Lc + Beta2[k][i];					
					tempL = __max(L21_d,y);
					absx_y = (L21_d>y)?(L21_d-y):(y-L21_d);
					if (LUTFlag)
						L21_d = tempL + MaxStar_lookup[MaxStar_index(absx_y)];
					else
						L21_d = tempL + log(1 + exp(-1.0 * absx_y));
				}
				L21[k] = L21_n - L21_d;
			}

			for (k=1;k<BlockSize+1;k++)
			{
				L1[k] = Lc*CodedBits[3*k-2] + L21[Deinterleaver[k]] + L12[k];				
				
				if (L1[k] > 0)
					DecodedBits[k]=1;
				else
					DecodedBits[k]=0;
			}
			
			ErrorFlag = 0;
			for (i=1; i<BlockSize+1; i++)
				if (SourceBits[i] != DecodedBits[i])			
				{
					ErrorFlag = 1;			
					break;
				}	
		}				

		//printf("%d\n",times);

		AverageIterNum = AverageIterNum + times;
		
		return;
	}
};