#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <memory.h>
#include "turbo.h"

// for rand2()
#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
#define PI 3.1415926
long id2=-16000;
long id3=-17000;

double rand2(long *idnum);
double AWGN(double mean, double std_dev);

int main()
{
	int Max_Error;
	long i,j;	
	double EbN0;	// dB
	double EbN0Pow;
	double noise;	// sigma
	
	int looptimes;	

	double Pb;
	int ErrorFlag;
	int ErrorBits,ErrorBlocks;	

	int SimPoints;
	float StartEbN0;

	printf("Please input the Max Errors\n");
	scanf("%d",&Max_Error);

	printf("Please input the Start EbN0\n");
	scanf("%f",&StartEbN0);

	printf("Please input the number of EbN0 points\n");
	scanf("%d",&SimPoints);
	
	#ifdef PUNCTURE
		TURBOCODE TurboCode = TURBOCODE(BLOCKSIZE, CONSTRAINT, POLY1, POLY2, CODERATE, 
									ITERATION, 1, LOOKUPTABLE);
	#else
		TURBOCODE TurboCode = TURBOCODE(BLOCKSIZE, CONSTRAINT, POLY1, POLY2, CODERATE, 
									ITERATION, 0, LOOKUPTABLE);
	#endif
	

	FILE * errorfile;
	errorfile = fopen("error.dat", "w");
	fclose(errorfile);

	FILE * result;
	result=fopen("result.dat", "w");
	fclose(result);	

	for (j=0;j<SimPoints;j++)
	{
		ErrorBits = 0;
		ErrorBlocks = 0;
		EbN0 = StartEbN0 + 0.2*j;
		EbN0Pow = pow(10.0,EbN0/10.0);		
		noise = 0.5/(EbN0Pow*TurboCode.RATE);
		noise = sqrt(noise);   //sigma
		
		TurboCode.EbN0 = EbN0;	
		TurboCode.Noise = noise;
		
		srand(time(NULL));		
		looptimes = 0;
		ErrorBits = 0;
		ErrorBlocks = 0;	
		
		TurboCode.AverageIterNum = 0;		
		
		while ( (ErrorBlocks < Max_Error) )		
		{
			looptimes++;
			if (looptimes % 500 == 0)
			{
				printf("%d blocks decoded, Average %f Iterationv\n",looptimes,TurboCode.AverageIterNum/looptimes);
			}			
			
			for (i=1; i<=TurboCode.BlockSize; i++)
				TurboCode.SourceBits[i] = rand() & 1;				
			
			TurboCode.encoder();
			
			if (TurboCode.PunctureFlag)
			{
				int k = 0;
				for (i=1;i<TurboCode.BlockSize+1;i++)
				{
					TurboCode.CodedBits[i*3-2] = AWGN(TurboCode.EncodedBits[i*2-1],noise);					
					k++;
					if (k%2 != 0)			
					{
						TurboCode.CodedBits[i*3-1] = AWGN(TurboCode.EncodedBits[i*2],noise);
						TurboCode.CodedBits[i*3] = 0.0;
					}
					else
					{
						TurboCode.CodedBits[i*3-1] = 0.0;
						TurboCode.CodedBits[i*3] = AWGN(TurboCode.EncodedBits[i*2],noise);						
					}
				}
			}
			else			
				for (i=1;i<TurboCode.BlockSize+1;i++)
				{
					TurboCode.CodedBits[i*3-2] = AWGN(TurboCode.EncodedBits[i*3-2],noise);
					TurboCode.CodedBits[i*3-1] = AWGN(TurboCode.EncodedBits[i*3-1],noise);
					TurboCode.CodedBits[i*3] = AWGN(TurboCode.EncodedBits[i*3],noise);
				}
			
			TurboCode.log_BCJR_decoder();

			ErrorFlag = 0;
			
			FILE* errorfile;
			errorfile=fopen("error.dat","a");
				
			for (i=1;i<=TurboCode.BlockSize;i++)
			{
				if ( !(TurboCode.SourceBits[i] == TurboCode.DecodedBits[i]) )
				{
					ErrorBits++;
					ErrorFlag = 1;
					fprintf(errorfile,"EbN0 = %f, No %d Block, No. %d ErrorBits at No. %d Bits of No. %d ErrorBlocks\n",EbN0,looptimes,ErrorBits,i,ErrorBlocks+1);
					printf("EbN0 = %f, No %d Block, No. %d ErrorBits at No. %d Bits of No. %d ErrorBlocks\n",EbN0,looptimes,ErrorBits,i,ErrorBlocks+1);
				}
			}

			if (ErrorFlag == 1)
			{				
				ErrorBlocks ++;				
				Pb = (double) (ErrorBits) / (double) (looptimes*TurboCode.BlockSize);
			}					

			fclose(errorfile);			
		}
		
		TurboCode.AverageIterNum = TurboCode.AverageIterNum / looptimes;

		FILE*  result;
		result = fopen("result.dat","a");
		fprintf(result,"%f	%f	%f	%f\n",EbN0,Pb,Max_Error/1.0/looptimes,TurboCode.AverageIterNum);
		fclose(result);		
	}	

	return 1;
}

// function rand2 follows
double rand2(long *idnum)
{
	int j;
	long k;
	static long idnum2=123456789;
	static long iy=0;
	static long iv[NTAB];
	double temp;
	if (*idnum<=0)
	{
		if (-(*idnum)<1) *idnum=1;
		else *idnum=-(*idnum);
		idnum2=(*idnum);
		for (j=NTAB;j>=0;j--)
		{
			k=(*idnum)/IQ1;
			*idnum=IA1*(*idnum-k*IQ1)-k*IR1;
			if (*idnum<0) *idnum+=IM1;
			if (j<NTAB) iv[j]=*idnum;
		}
		iy=iv[0];
	}
	k=(*idnum)/IQ1;
	*idnum=IA1*(*idnum-k*IQ1)-k*IR1;
	if (*idnum<0) *idnum+=IM1;
	k=idnum2/IQ2;
	idnum2=IA2*(idnum2-k*IQ2)-k*IR2;
	if (idnum2<0) idnum2+=IM2;
	j=iy/NDIV;
	iy=iv[j]-idnum2;
	iv[j]=*idnum;
	if (iy<1) iy+=IMM1;
	if ((temp=AM*iy)>RNMX) return RNMX;
	else return temp;
}

double AWGN(double mean, double std_dev)
{
	double rand2(long *idnum);	
	double u1,u2,a,b;
	u1=(double) (rand2(&id2));
	u2=(double) (rand2(&id3)); 
	a=2*PI*u1;
	b=sqrt(-2.0*log(u2));
	
	return mean + std_dev * b * cos(a);
}