/**********************************************/
/* De-modulator
/* double * DemodulatingF (struct BasicParaS * ctrl, struct Complex * input, double Eb_No_dB)
/*      Written by: Ouyang Ziyue,
/*            Date: Dec 21st, 2007,
/*        Function: It de-modulates the input complex symbols, outputs the LLR of bits
/* Input parameter:
/*        The modulation type is included in the ctrl,
/*        The input includes all the input complex symbols with AWGN,
/*        The Eb_No_dB is used to define the SNR or Eb/No.
/* Output parameter:
/*        An DOUBLE set which consists all the LLR of bits is outputed
/* Note:
/*        ctrl should be built before this function is called,
/**********************************************/

#include "parameter_sets.h"

const double sqrt2 = 1.414213562373095;

double * DemodulatingF (struct BasicParaS * ctrl, struct Complex * input, double Eb_No_dB)
{
	//////////////////////////////////////////////////////////////////////////
	//Declaration
	double * output;
	double inv_sigma2;
	double currentI, currentQ;
	int i,j;
	double temp[6];
	
	//////////////////////////////////////////////////////////////////////////
	//Allocating
	i = ctrl->numModuOut*ctrl->bitsPerSym;
	output = new double[i];

	//////////////////////////////////////////////////////////////////////////
	// Set the parameter about the noise
	if (ctrl->noiseMode == 0)
	{
			if(ctrl->typeModu == 0){
			inv_sigma2 = 2*pow(10.0,(Eb_No_dB/10))/ctrl->r;
		}	else if(ctrl->typeModu == 1){
			inv_sigma2 = 2*pow(10.0,(Eb_No_dB/10))*2/ctrl->r;
		}	else if(ctrl->typeModu == 2){
			inv_sigma2 = 2*pow(10.0,(Eb_No_dB/10))*4/ctrl->r;
		}	else if(ctrl->typeModu == 3){
			inv_sigma2 = 2*pow(10.0,(Eb_No_dB/10))*6/ctrl->r;
		}
	}
	else
	{
			inv_sigma2 = 2*pow(10.0,(Eb_No_dB/10));
	}

	
	//////////////////////////////////////////////////////////////////////////
	//Soft Demodulating
	switch(ctrl->typeDecode)
	{
	case 0: // SPA
		switch (ctrl->typeModu) {
		case 0: // BPSK
			for (i=0; i<ctrl->numModuOut; i++) {
				currentI = (*(input+i)).real;
				currentQ = (*(input+i)).imag;
				// calculate the soft output
				*(output+i) = 1/(1+exp(2*inv_sigma2*currentI)); // the probability of 1
			}
			break;
		case 1: // QPSK
			for (i=0; i<ctrl->numModuOut; i++) {
				currentI = (*(input+i)).real;
				currentQ = (*(input+i)).imag;
				// calculate the soft output
				*(output+2*i) = 1/(1+exp(sqrt2*inv_sigma2*currentI));
				*(output+2*i+1) = 1/(1+exp(sqrt2*inv_sigma2*currentQ));
			}
			break;
		case 2: // 16QAM
			printf("SPA is not supported in 16QAM!\n");
			exit(EXIT_FAILURE);
			break;
		case 3: // 64QAM
			printf("SPA is not supported in 64QAM!\n");
			exit(EXIT_FAILURE);
			break;
		default: // unsupported modulation types
			exit(EXIT_FAILURE);
		}
		break;
	case 1: // LSPA
		switch (ctrl->typeModu) {
		case 0: // BPSK
			for (i=0; i<ctrl->numModuOut; i++) {
				currentI = (*(input+i)).real;
				currentQ = (*(input+i)).imag;
				// calculate the soft output
				*(output+i) = 2*currentI*inv_sigma2; // the logarithm of probability of 0/1
			}
			break;
		case 1: // QPSK
			for (i=0; i<ctrl->numModuOut; i++) {
				currentI = (*(input+i)).real;
				currentQ = (*(input+i)).imag;
				// calculate the soft output
				*(output+2*i) = sqrt2*currentI*inv_sigma2;
				*(output+2*i+1) = sqrt2*currentQ*inv_sigma2;
			}
			break;
		case 2: // 16QAM
			for (i=0; i<ctrl->numModuOut; i++) {
				currentI = (*(input+i)).real;
				currentQ = (*(input+i)).imag;
				if (MaxLogDeMap(0, 2, currentI, temp) != 0)
					exit(EXIT_FAILURE);
				if (MaxLogDeMap(2, 2, currentQ, temp) != 0)
					exit(EXIT_FAILURE);
				for (j=0; j<ctrl->bitsPerSym; j++) {
					*(output+ctrl->bitsPerSym*i+j) = *(temp+j);
				}
			}
			break;
		case 3: // 64QAM
			for (i=0; i<ctrl->numModuOut; i++) {
				currentI = (*(input+i)).real;
				currentQ = (*(input+i)).imag;
				if (MaxLogDeMap(0, 3, currentI, temp) != 0)
					exit(EXIT_FAILURE);
				if (MaxLogDeMap(2, 3, currentQ, temp) != 0)
					exit(EXIT_FAILURE);
				for (j=0; j<ctrl->bitsPerSym; j++) {
					*(output+ctrl->bitsPerSym*i+j) = *(temp+j);
				}
			}
			break;
		default: // unsupported modulation types
			exit(EXIT_FAILURE);
		}
		break;
	case 2: // MSA
		switch (ctrl->typeModu) {
		case 0 : // BPSK
			for (i=0; i<ctrl->numModuOut; i++) {
				currentI = (*(input+i)).real;
				currentQ = (*(input+i)).imag;
				// calculate the soft output
				*(output+i) = 2*currentI*inv_sigma2;
			}
			break;
		case 1: // QPSK
			for (i=0; i<ctrl->numModuOut; i++) {
				currentI = (*(input+i)).real;
				currentQ = (*(input+i)).imag;
				// calculate the soft output
				*(output+2*i) = sqrt2*currentI*inv_sigma2;
				*(output+2*i+1) = sqrt2*currentQ*inv_sigma2;
			}
			break;
		case 2: // 16QAM
			for (i=0; i<ctrl->numModuOut; i++) {
				currentI = (*(input+i)).real;
				currentQ = (*(input+i)).imag;
				if (MaxLogDeMap(0, 2, currentI, temp) != 0)
					exit(EXIT_FAILURE);
				if (MaxLogDeMap(2, 2, currentQ, temp) != 0)
					exit(EXIT_FAILURE);
				for (j=0; j<ctrl->bitsPerSym; j++) {
					*(output+ctrl->bitsPerSym*i+j) = *(temp+j);
				}
			}
			break;
		case 3: // 64QAM
			for (i=0; i<ctrl->numModuOut; i++) {
				currentI = (*(input+i)).real;
				currentQ = (*(input+i)).imag;
				if (MaxLogDeMap(0, 3, currentI, temp) != 0)
					exit(EXIT_FAILURE);
				if (MaxLogDeMap(2, 3, currentQ, temp) != 0)
					exit(EXIT_FAILURE);
				for (j=0; j<ctrl->bitsPerSym; j++) {
					*(output+ctrl->bitsPerSym*i+j) = *(temp+j);
				}
			}
			break;
		default: // unsupported modulation types
			exit(EXIT_FAILURE);
		}
		break;
	default: // unsupported decoding algorithms
		exit(EXIT_FAILURE);		
	}


	//////////////////////////////////////////////////////////////////////////
	//DEBUG
#ifdef DEBUG
	printf("The output of soft demodulation is...\n");
	int h;
	h=0;
	for (i=0; i<ctrl->codeN; i++)
	{
		printf(" %2f", *(output+i));
		h++;
		if (h==8)
		{
			printf("\n");
			h=0;

		}
	}
#endif

	return output;
}