#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "random.h"
// make the complex structure
typedef struct complex1
{
	double real;
	double imag;
}complex;
// define the number of the iteration
#define ITERATION pow(10, 4)
#define THENUMBEROFDATA 1000

void convolution(int *data, int *codeddata);
double modulation(int data);
void demodulation(double * r_data);
int bercheck(int data, double d_data);
double comparison(double com1, double com2);
void viterbi_decode(int *d_data, int *decode_data);
void quantization(double r_data, int *d_data);
double dis(int data, int data1, int data2, int value1, int value2);

void main()
{
	FILE *fp1;
	double Noisepower;
	complex Noise;
	int data[THENUMBEROFDATA];
	int codeddata[2*(THENUMBEROFDATA+2)];
	int cr_data[2];
	int ci_data[2];
	complex m_data;
	complex r_data;
//	double d_datareal[2];
//	double d_dataimag[2];
	int d_data[3];
	int d_data1[3];
	int dd_data[6*(THENUMBEROFDATA+2)];
	int decode_data[THENUMBEROFDATA+2];
	int error;
//	int error1;
//	int error2;
//	int error3;
	double avgerr;
	int EbofN0;
	
	// open and make the text file 
	fp1=fopen("result_16qam_ViterbiSoft.txt","w");
	srand(time(0));
	for(EbofN0=11; EbofN0<=15; EbofN0+=1)
	{
		double sumerr=0.0;
		int run;		
		int i;
		int l;
	
		for(run=0;run<ITERATION;run++)
		{
			// generate the real part data 		
			generate_random(&data[0], THENUMBEROFDATA);
			
			convolution(data, codeddata);
			// generate the imaginary part data
            //generate_random(&data1[0], THENUMBEROFDATA/2);
			
            for (i=0; i<=(THENUMBEROFDATA/2);i++)
			{
				
				cr_data[0]=codeddata[4*i];
				cr_data[1]=codeddata[4*i+1];
                ci_data[0]=codeddata[4*i+2];
				ci_data[1]=codeddata[4*i+3];
				//printf("%d",cr_data[0]);
				//printf("%d",cr_data[1]);
				//printf("%d",ci_data[0]);
				//printf("%d\n",ci_data[1]);
				// modulated the real part signal
				m_data.real = modulation(cr_data);
				// modulated the imaginary part signal 
				m_data.imag = modulation(ci_data);

				// make the noisepower	
				Noisepower=1/(2*pow(10,(double)EbofN0/10));
				// make the real part of noise
				Noise.real = sqrt(Noisepower/2)*GaussRand();
				// make the imaginary part of noise
				Noise.imag = sqrt(Noisepower/2)*GaussRand();
			
				// receive the real part of the signal 
				r_data.real = sqrt(1.0/10.0)*m_data.real + Noise.real;
				// receive the imaginary part of the signal 
				r_data.imag = sqrt(1.0/10.0)*m_data.imag + Noise.imag;
						
				quantization(r_data.real,d_data);
				quantization(r_data.imag,d_data1);
                
				dd_data[6*i]=d_data[0];
				dd_data[6*i+1]=d_data[1];
				dd_data[6*i+2]=d_data[2];
				dd_data[6*i+3]=d_data1[0];
				dd_data[6*i+4]=d_data1[1];
				dd_data[6*i+5]=d_data1[2];
				//printf("%d",dd_data[6*i]);
				//printf("%d",dd_data[6*i+1]);
				//printf("%d",dd_data[6*i+2]);
				//printf("%d",dd_data[6*i+3]);
				//printf("%d",dd_data[6*i+4]);
				//printf("%d\n",dd_data[6*i+5]);
				//d_datareal[0] = r_data.real; 
				//d_dataimag[0] = r_data.imag;
				// demodulated the real part of the signal 
				//demodulation(&d_datareal);
				// demodulated the imaginary part of the signal 
				//demodulation(&d_dataimag);
				//d_data[4*i]=(int)d_datareal[0];
				//d_data[4*i+1]=(int)d_datareal[1];
				//d_data[4*i+2]=(int)d_dataimag[0];
				//d_data[4*i+3]=(int)d_dataimag[1];

			}

			viterbi_decode(dd_data,decode_data);
			// check the error 

			for (l=0;l<=(THENUMBEROFDATA-1);l++)
			{
				//printf("%d\n",decode_data[l]);
				//printf("%d\n",data[l]);
				error= bercheck(data[l], decode_data[l]);
				sumerr = sumerr + error;
			}
			//error= bercheck(data[0], d_datareal[0]);
			//error1= bercheck(data[1], d_datareal[1]);
			//error2= bercheck(data1[0], d_dataimag[0]);
			//error3= bercheck(data1[1], d_dataimag[1]);

			//sumerr = sumerr + error + error1 + error2 + error3;
		}
		//calculated the BER
		avgerr = sumerr / (ITERATION*THENUMBEROFDATA);
		fprintf(fp1,"%.7f\n",avgerr);		

		printf("%.7f\n", avgerr);
	}
	fclose(fp1);
}

void quantization(double r_data, int *d_data)
{
	if (r_data<=-sqrt(1.0/10.0)*3*3/4)
	{
		d_data[0]=0;
		d_data[1]=0;
		d_data[2]=0;
	}else if (r_data<=-sqrt(1.0/10.0)*3*2/4)
	{
		d_data[0]=0;
		d_data[1]=0;
		d_data[2]=1;
	}else if (r_data<=-sqrt(1.0/10.0)*3*1/4)
	{
		d_data[0]=0;
		d_data[1]=1;
		d_data[2]=0;
	}else if (r_data <= 0 )
	{
		d_data[0]=0;
		d_data[1]=1;
		d_data[2]=1;
	}else if (r_data<=sqrt(1.0/10.0)*3*1/4)
	{
		d_data[0]=1;
		d_data[1]=0;
		d_data[2]=0;
	}else if (r_data<=sqrt(1.0/10.0)*3*2/4)
	{
		d_data[0]=1;
		d_data[1]=0;
		d_data[2]=1;
	}else if (r_data<=sqrt(1.0/10.0)*3*3/4)
	{
		d_data[0]=1;
		d_data[1]=1;
		d_data[2]=0;
	}else 
	{
		d_data[0]=1;
		d_data[1]=1;
		d_data[2]=1;
	}
}
// modulation function
double modulation(int *data)
{
	int m_data;
	if (*data==1)
	{
		if(*(data+1)==1)
		{
			m_data=1;
		}
		else
		{
			m_data=3;
		}
	}
	else
	{
		if(*(data+1)==1)
		{
			m_data=-1;
		}
		else
		{
			m_data=-3;
		}		
	}
	return (m_data);
}
// demodulation function
void demodulation(double *r_data)
{
	if (r_data[0]>=0.0)
	{
		if (r_data[0]>=sqrt(1.0/10.0)*2.0)
		{
			r_data[0]=1;
     		r_data[1]=0;
		}
		else
		{
			r_data[0]=1;
			r_data[1]=1;
		}
	}
	else
	{
		if (r_data[0]<=sqrt(1.0/10.0)*(-2.0))
		{
			r_data[0]=0;
     		r_data[1]=0;
		}
		else
		{
			r_data[0]=0;
     		r_data[1]=1;
		}
	}
}
// error check function
int bercheck(int data, double d_data)
{
	int error=0;
	if ((double)data==d_data)
	{
		error=0;
	}
	else
	{
		error=1;
		//printf("error occured\n");
	}
	return (error);
}

void convolution(int *data, int *codeddata)
{
	int a = data[0];
	int b = 0;
	int c = 0;
	int N =0;
	//int codeddata[12];
	for(N=0;N<=(THENUMBEROFDATA+1);N++)
	{
		codeddata[2*N] = a^b^c;
		codeddata[2*N+1] = a^c;
		a = data[N+1];
		b = data[N];
		if(N == 0)
		{
			c = 0;
		}else
		{
			c = data[N-1];
		}
		if(N>=(THENUMBEROFDATA-1))
		{
			a=0;
		}
		if(N>=THENUMBEROFDATA)
		{
			b=0;
		}
	}
}

void viterbi_decode(int *d_data, int *decode_data)
{
	double s[4][(THENUMBEROFDATA+3)];
	int k1;
	int k2;
	int c[2*((THENUMBEROFDATA+2))];
	int l1;
//	int l2;
	int i;
	int l;
	for (i = 0;i<=3;i++)
	{
		for(l=0;l<=THENUMBEROFDATA+2;l++)
		{
			s[i][l]=10000000000000;
		}
	}
    s[0][0]=0;
	s[0][1]=dis(d_data[0],d_data[1],d_data[2],0,0);
	s[1][1]=dis(d_data[0],d_data[1],d_data[2],1,1);
	s[0][2]=s[0][1]+dis(d_data[3],d_data[4],d_data[5],0,0);
	s[1][2]=s[0][1]+dis(d_data[3],d_data[4],d_data[5],1,1);
	s[2][2]=s[1][1]+dis(d_data[3],d_data[4],d_data[5],1,0);
	s[3][2]=s[1][1]+dis(d_data[3],d_data[4],d_data[5],0,1);
	
	for(k2=3;k2<=(THENUMBEROFDATA+2);k2++)
	{
		for(k1=0;k1<=3;k1++)
		{
			if(k1==0)
			{
				s[k1][k2]=comparison(s[k1][k2-1]+dis(d_data[3*(k2-1)],d_data[3*(k2-1)+1],d_data[3*(k2-1)+2],0,0), s[k1+2][k2-1]+dis(d_data[3*(k2-1)],d_data[3*(k2-1)+1],d_data[3*(k2-1)+2],1,1));
			}
			else if(k1==1)
			{
				s[k1][k2]=comparison(s[k1-1][k2-1]+dis(d_data[3*(k2-1)],d_data[3*(k2-1)+1],d_data[3*(k2-1)+2],1,1), s[k1+1][k2-1]+dis(d_data[3*(k2-1)],d_data[3*(k2-1)+1],d_data[3*(k2-1)+2],0,0));
			}
			else if(k1==2)
			{
				s[k1][k2]=comparison(s[k1-1][k2-1]+dis(d_data[3*(k2-1)],d_data[3*(k2-1)+1],d_data[3*(k2-1)+2],1,0), s[k1+1][k2-1]+dis(d_data[3*(k2-1)],d_data[3*(k2-1)+1],d_data[3*(k2-1)+2],0,1));
			}
			else if(k1==3)
			{
				s[k1][k2]=comparison(s[k1-2][k2-1]+dis(d_data[3*(k2-1)],d_data[3*(k2-1)+1],d_data[3*(k2-1)+2],0,1), s[k1][k2-1]+dis(d_data[3*(k2-1)],d_data[3*(k2-1)+1],d_data[3*(k2-1)+2],1,0));
			}
		}
	}
    
	if((s[0][(THENUMBEROFDATA+1)]+dis(d_data[3*(THENUMBEROFDATA+1)],d_data[3*(THENUMBEROFDATA+1)+1],d_data[3*(THENUMBEROFDATA+1)+2],0,0))==comparison(s[0][(THENUMBEROFDATA+1)]+dis(d_data[3*(THENUMBEROFDATA+1)],d_data[3*(THENUMBEROFDATA+1)+1],d_data[3*(THENUMBEROFDATA+1)+2],0,0) , s[2][(THENUMBEROFDATA+1)]+dis(d_data[3*(THENUMBEROFDATA+1)],d_data[3*(THENUMBEROFDATA+1)+1],d_data[3*(THENUMBEROFDATA+1)+2],1,1)))
	{
		c[(THENUMBEROFDATA+1)] = 0;
		decode_data[(THENUMBEROFDATA-1)]=0;
	}
	else 
	{
		c[(THENUMBEROFDATA+1)] = 2;
		decode_data[(THENUMBEROFDATA-1)]=0;
	}
	for (l1=(THENUMBEROFDATA+1);l1>=1;l1--)
	{
		if (c[l1]==0)
		{
			if((s[0][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],0,0))==comparison(s[0][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],0,0),s[2][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],1,1)))
			{
				c[l1-1] = 0;
				decode_data[l1-1]=0;
			}
			else 
			{
				c[l1-1] = 2;
				decode_data[l1-1]=0;
			}
		}else if(c[l1]==1)
		{
			if((s[0][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],1,1))==comparison(s[0][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],1,1),s[2][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],0,0)))
			{
				c[l1-1] = 0;
				decode_data[l1-1]=1;
			}
			else 
			{
				c[l1-1] = 2;
				decode_data[l1-1]=1;
			}
		}else if(c[l1]==2)
		{
			if((s[1][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],1,0))==comparison(s[1][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],1,0),s[3][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],0,1)))
			{
				c[l1-1] = 1;
				decode_data[l1-1]=0;
			}
			else 
			{
				c[l1-1] = 3;
				decode_data[l1-1]=0;
			}
		}else if(c[l1]==3)
		{
			if((s[1][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],0,1))==comparison(s[1][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],0,1),s[3][l1-1]+dis(d_data[3*(l1-1)],d_data[3*(l1-1)+1],d_data[3*(l1-1)+2],1,0)))
			{
				c[l1-1] = 1;
				decode_data[l1-1]=1;
			} 
			else 
			{
				c[l1-1] = 3;
				decode_data[l1-1]=1;
			}
		}
	}
	/*
		for(l=0;l<=THENUMBEROFDATA+1;l++)
		{
			printf("%d ",c[l]);
		}
	printf("\n");
	
	for(i=0;i<=3;i++)
	{
		for(l=0;l<=THENUMBEROFDATA+2;l++)
		{
			printf("%f ",s[i][l]);	
		}
		printf("\n");
	}*/
}

double comparison(double com1, double com2)
{
	double com3;
	if(com1>=com2)
	{
		com3=com2;
	}else
	{
		com3=com1;
	}
	return (com3);
}

double dis(int data, int data1, int data2, int value1, int value2)
{
	double distance;
	double value;
	if(value1==0)
	{
		if(value2==0)
		{
			value = 0;
		}
		else
		{
			value = 7/3;
		}
	}
	else
	{
		if(value2==0)
		{
			value = 7;
		}else
		{
			value = 14/3;
		}
	}
	distance = sqrt(pow(((data*4+data1*2+data2)-value),2));
	return (distance);
}