codec_function.cpp

DPCM编码的联合信源信道译码

#include "codec_function.h"
#include "UVLC.h"
#include "init_UVLC.h"
#include "sim_config.h"
/*
**********************************************************************************
  	FILENAME: ldpc_codec_function.c
	version 0.1
	(C) Pan Yu, 2004
	Declaration:
		All right is deserved by author
		Anyone use this source codes and remove author's name is prohibited   
		Without agreement,noboby can copy this codes to others. 
	Fucntion Description:
		All function need by simulate are included in this file!!
*********************************************************************************
        ChangeLog:
	2004-07-10 : start
	2004-07-11 : all debug complete
	2004-07-19 : add fixed decoding function and modified ADconverter()
	2004-08-01 : add BER statistic function for every iteration in decoding program
	2004-10-22 : add min-sum iteration function by ZhuJia
**********************************************************************************
*/


void  SourceGen(int i)
{
	int k;
	for (k=sys_UVLC.bpl[i];k<sys_spt.message_bits;k++)
	{
		sys_spt.source_bits[k] = 0;	  //add 0 to the end of UVLC edcoded sequence
	
	}
	
}


void LdpcEncoder()
{
	int i,k;
	
	int *codetemp;
	codetemp = (int *)malloc(sys_spt.col_num*sizeof(int));
	
	for(i=0;i<sys_spt.message_bits;i++)
	{
		sys_spt.source_bits_in[i] = sys_spt.source_bits[i];
	}

	for (i=0,k=0;i<sys_spt.col_num;i++)
	{
		//get the row number in i-th col
		int sum=0;
		while (sys_spt.G_sp_col[3*k+1]==i+1)
		{
			sum += sys_spt.source_bits[sys_spt.G_sp_col[3*k]-1];
			k++;		
		}
		codetemp[i] = (sum%2);
	}

	int col=0;
	int	flag_col=sys_spt.col_num-sys_spt.message_bits;
	for (i=0;i<sys_spt.col_num;i++)
	{
		if (sys_spt.message_locations[i]==1)
		{	
			sys_spt.source_code[i] = codetemp[flag_col];
			flag_col++;
		}
		else 
		{
			sys_spt.source_code[i] = codetemp[col];
			col++;
		}
	}

	free(codetemp);
}


void LdpcEncoder2()
{
	int i,k;
	
	int *codetemp;
	codetemp = (int *)malloc(sys_spt.col_num*sizeof(int));

	for (i=0,k=0;i<sys_spt.col_num;i++)
	{
		//get the row number in i-th col
		int sum=0;
		while (sys_spt.G_sp_col[3*k+1]==i+1)
		{
			sum += sys_spt.source_bits_in[sys_spt.G_sp_col[3*k]-1];
			k++;		
		}
		codetemp[i] = (sum%2);
	}

	int col=0;
	int	flag_col=sys_spt.col_num-sys_spt.message_bits;
	for (i=0;i<sys_spt.col_num;i++)
	{
		if (sys_spt.message_locations[i]==1)
		{	
			sys_spt.bit_nodes.llr_in[i] = codetemp[flag_col];
			flag_col++;
		}
		else 
		{
			sys_spt.bit_nodes.llr_in[i] = codetemp[col];
			col++;
		}

	}

	free(codetemp);
}




void Modulation(int blocknum)
{
	int i;
	for (i=0;i<sys_UVLC.bpl[blocknum];i++)
	{
//		sys_spt.send_code[i] = 1-2*sys_spt.source_bits[i];		//modulate BPSK
		sys_spt.send_code[i] = sys_spt.source_bits[i];	
	}
}



double randn(double sigma)
{
	double x,y;

	do
		x=(double)rand()/32767.;
	while(x<1.0e-10);		//note we need to eliminate number instability

 	y=(double)rand()/32767.;
 	
	return (sqrt(-2.*log(x))*cos(2.*3.1415926*y)*sigma);
}


double gallager_f(double x){
	double interm;                
	if (x < 0.00001)
		interm = 12.21; // 12.21=gallager_f(0.00001)
	else if (x > 30.)
		interm = 0.0;
	else 
		interm = log((exp(x)+1.0)/(exp(x)-1.0));
	return(interm);
}

void Init_LdpcIteration()
{
	int i,j;
	
	//----------------------------initialization for bit_nodes_llr
	for (i=0;i<sys_spt.col_num;i++)		
	{
		for (j=0;j<sys_spt.bit_nodes.size[i];j++)	
		{
			sys_spt.bit_nodes.llr[i][j] = sys_spt.llr_init[i];
		}
	}
	


}

void LdpcIteration(int block_num){
	int i;
	int j;
	int iter;
	int stop_flag;
	int bit_error_flag=0;
	int row_index;
	double llr_q;

	int correct_flag=0;

	int *S;
	S = (int *)calloc(sys_spt.row_num,sizeof(int));

	for (iter=1; iter<=sys_spt.iter_times; iter++) 
	{
		//------------------------clear check_nodes.llr	
		//------------------------set initial sign for check sum
		for (i=0;i<sys_spt.row_num;i++)
		{
			sys_spt.check_nodes.llr[i] = 0.	;	
			sys_spt.check_nodes.llr_sign[i] = 1;
		}
		//------------------------clear bit nodes llr sum
		for (i=0;i<sys_spt.col_num;i++)	
		{
			sys_spt.bit_nodes.llr_sum[i] = 0.;		 
		}
		//------------------------calculate  parity check  sum
		for (i=0;i<sys_spt.col_num;i++)		
		{
			for (j=0;j<sys_spt.bit_nodes.size[i];j++)
			{
				row_index = sys_spt.bit_nodes.index[i][j];
				llr_q = sys_spt.bit_nodes.llr[i][j];
				sys_spt.check_nodes.llr[row_index] += gallager_f(fabs(llr_q)); 
				sys_spt.check_nodes.llr_sign[row_index] *= ((llr_q>=0)?1:-1);
			}		
		}		
		//------------------------calculate   check nodes llr and sum bit_nodes.llr
		for (i=0;i<sys_spt.col_num;i++)		
		{
			for (j=0;j<sys_spt.bit_nodes.size[i];j++)
			{
				row_index = sys_spt.bit_nodes.index[i][j];
				llr_q = sys_spt.bit_nodes.llr[i][j];
				sys_spt.bit_nodes.llr_sign[i][j] = sys_spt.check_nodes.llr_sign[row_index]*((llr_q>=0)?1:-1);

				sys_spt.bit_nodes.llr[i][j] = gallager_f(sys_spt.check_nodes.llr[row_index] - gallager_f(fabs(llr_q)))\
						*sys_spt.bit_nodes.llr_sign[i][j];

				sys_spt.bit_nodes.llr_sum[i] += sys_spt.bit_nodes.llr[i][j]; 		
			}		
		}
		//------------------------update bit_nodes.llr
		for (i=0;i<sys_spt.col_num;i++)
		{
			for (j=0;j<sys_spt.bit_nodes.size[i];j++)
			{
				row_index = sys_spt.bit_nodes.index[i][j];

				sys_spt.bit_nodes.llr[i][j] = sys_spt.bit_nodes.llr_sum[i] - \
					sys_spt.bit_nodes.llr[i][j] + sys_spt.llr_init[i];				
				sys_spt.bit_nodes.llr_out[i] = ((sys_spt.bit_nodes.llr_sum[i] + sys_spt.llr_init[i])>=0)?0:1;
//				sys_spt.bit_nodes.llr_out[i] = sys_spt.bit_nodes.llr_sum[i] + sys_spt.llr_init[i]; 
			}		
		}

		//-------------------------re-generate the decoded source bits

		for (j=0,i=0;i<sys_spt.col_num;i++)
		{
			if (sys_spt.message_locations[i]==1)
			{	
				sys_spt.source_bits_out[j]=sys_spt.bit_nodes.llr_out[i];
				sys_UVLC.received[j] = sys_spt.source_bits_out[j];  //soft information
				j++;

				
			}
			
			
		}

		
		//------------------------use source_bits to compute the error

		for (i=0;i<sys_spt.message_bits;i++)
		{
			if (sys_spt.source_bits_in[i]!=sys_spt.source_bits_out[i])
			{
					sys_spt.error_iter[iter-1]++;
			}
		}

				
		//------------------------stop iteration judgement
		for (i=0;i<sys_spt.row_num;i++)
		{
			S[i] = 0;				//clear check sum hard value
		}
		for (i=0;i<sys_spt.col_num;i++)		
		{			
			for (j=0;j<sys_spt.bit_nodes.size[i];j++)
			{
				S[sys_spt.bit_nodes.index[i][j]] += sys_spt.bit_nodes.llr_out[i]; 
			}
		}
		stop_flag = 0;
		for (i=0;i<sys_spt.row_num;i++)
		{
			if (S[i]%2!=0)
			{
				stop_flag = 1;
				break;
			}
		}
		if (!stop_flag)		//means all check sums are satisfied
		{
			correct_flag=1;
			break;
		}

		//Ldpc_Ber_Cal();

////////////////////////////////////joint uvlc decoding and update the ldpc		
/*        if((iter%2)==0)
		{
          UVLCDecoder(block_num);
        //UVLC_Ber_Cal(block_num);
		//printf("\n");
          UpdateBit();
		}
*/
/////////////////////////////////////////////////////
	}

	free(S);
//	return(correct_flag);
}



void PassChannel(int blocknum){
	int i,k;
	double receive_signal;
	int r;
	int e,d=0;

	e=1/BSC_P;
	for (i=0;i<sys_UVLC.bpl[blocknum];i++)
	{
		for(k=0;k<BSC_PASS;k++)
		{
			r=rand();
			r=r%e;
		}
	   if(r==1)	  sys_spt.received[i]=(sys_spt.send_code[i]+1)%2;
	   else		  sys_spt.received[i]=sys_spt.send_code[i];

//	   if(sys_spt.received[i]==0)  sys_spt.llr_init[i]=log10((1-BSC_P)/BSC_P);
//	   else	                       sys_spt.llr_init[i]=log10(BSC_P/(1-BSC_P));

	}

/*	sys_spt.noisevar = 1/(2*sys_spt.code_rate*pow(10,sys_spt.EbNo/10.));

/*    for(i=0;i<sys_spt.message_bits;i++)
	{
     	receive_signal= sys_spt.send_code[i] + randn(sqrt(sys_spt.noisevar));
        sys_UVLC.received[i]=receive_signal;
	}

	for (i=0;i<sys_spt.col_num;i++){
		receive_signal		= sys_spt.send_code[i] + randn(sqrt(sys_spt.noisevar));
		sys_spt.llr_init[i] = 2.0*receive_signal/sys_spt.noisevar;		//calculate the initial inputLLR
	}
*/
}

void BerCalc(long blocknum){
	int i;
	int block_error_flag = 0;

	sys_UVLC.sym_error = 0; 

//	FILE *fq;

//	fq=fopen("temp2.txt","wb");

//   sys_spt.bit_error=0;
//	sys_spt.hard_error=0;

	//----------------------error statistic
	for (i=0;i<sys_spt.col_num;i++)
	{
		if (sys_spt.source_code[i]!=((sys_spt.llr_init[i]>=0)?0:1))
		{
			sys_spt.hard_error++;
		}
	}



//	if(flag)
	{
		for(i=0;i<sys_UVLC.bpl[blocknum];i++)
		{
			if (sys_spt.source_bits[i]!=sys_spt.source_bits_in[i])  //zhynxn SSCD-2
		//	if (sys_spt.source_bits[i]!=sys_spt.source_bits_out[i])   //SSCD-1
			{
				sys_spt.bit_error++;
				block_error_flag = 1;
			}
		}

		for(i=0;i<256;i++)
		{
			if(sys_UVLC.decoded[i] != sys_UVLC.encoded[i])
				sys_UVLC.sym_error++;
		}	
	}


    sys_UVLC.sum_symbol_error+=sys_UVLC.sym_error; 
//	sys_UVLC.sym_ber = (double)sys_UVLC.sym_error/128;

	for (i=0;i<sys_spt.iter_times;i++)
	{
		sys_spt.ber_iter[i] = (double)sys_spt.error_iter[i]/(sys_spt.col_num);	
	}
//	sys_spt.ber			= (double)sys_spt.bit_error/(sys_spt.col_num);
//	sys_spt.bler		= (double)sys_spt.block_error;
//	sys_spt.hard_ber	= (double)sys_spt.hard_error/(sys_spt.col_num);
//	sys_spt.undet_ber	= (double)sys_spt.undetected_bit_error/(sys_spt.col_num);
//	sys_spt.undet_bler	= (double)sys_spt.undetected_block_error;

//	sys_spt.fp_result= fopen(RESULT_FILE,"w");

	printf("BSC_P:%1.4f blockCnt=%3d harderror =%5d  biterror =%3d  symerror = %3d \n",\
		BSC_P,blocknum,sys_spt.hard_error,sys_spt.bit_error,sys_UVLC.sym_error);

	fprintf(sys_spt.fp_result,"BSC_P:%1.4f blockCnt=%3d harderror =%5d  biterror =%3d  symerror = %3d \n",\
		BSC_P,blocknum,sys_spt.hard_error,sys_spt.bit_error,sys_UVLC.sym_error);


/*	fprintf(sys_spt.fp_result, "SNR:%1.2f blockCnt =%3d  biterror =%3d  bitBER =%1.4e symerror = %3d  symBER =%1.4e\n",\
		sys_spt.EbNo,blocknum,sys_spt.bit_error,sys_spt.ber,sys_UVLC.sym_error,sys_UVLC.sym_ber);
*/
//	fprintf(sys_spt.fp_result, "BSC_P:%1.4f blockCnt =%3d  biterror =%3d bitBER =%1.4e symerror = %3d  symBER =%1.4e\n",\
		BSC_P,blocknum,sys_spt.bit_error,sys_spt.ber,sys_UVLC.sym_error,sys_UVLC.sym_ber);
	
/*	printf("SNR:%1.2f  blockCnt =%6d  error =%5d  Undet = %3d  BER =%1.4e\n",\
		sys_spt.EbNo,blocknum,sys_spt.bit_error,sys_spt.undetected_bit_error,sys_spt.ber);
	sys_spt.fp_result= fopen(RESULT_FILE,"a+");
	fprintf(sys_spt.fp_result, "SNR:%1.2f  blockCnt =%6d  error =%5d  Undet = %3d  BER =%1.4e\n",\
		sys_spt.EbNo,blocknum,sys_spt.bit_error,sys_spt.undetected_bit_error,sys_spt.ber);
*/
}

void Ldpc_Ber_Cal(int blocknum)
{
	int i,j=0;
	int block_error_flag = 0;

//	sys_spt.bit_error=0;
//	sys_spt.hard_error=0;
	
		//----------------------error statistic
	for (i=0;i<sys_spt.col_num;i++)
	{
		if (sys_spt.source_code[i]!=((sys_spt.llr_init[i]>=0)?0:1))
		{
			sys_spt.hard_error++;
		}
        if (sys_spt.source_code[i]!=sys_UVLC.received[i])
		{
			j++;
		}
		if (sys_spt.source_code[i]!=sys_spt.bit_nodes.llr_out[i])
		{
//			printf("%d",sys_spt.source_bits_out[i]);

			sys_spt.bit_error++;
			block_error_flag = 1;
		}
	}
	if (block_error_flag==1)
	{
		sys_spt.block_error++;
	}

//	for (i=0;i<sys_spt.iter_times;i++)
//	{
//		sys_spt.ber_iter[i] = (double)sys_spt.error_iter[i]/(blocknum*sys_spt.col_num);	
//	}
//	sys_spt.ber			= (double)sys_spt.bit_error/(blocknum*sys_spt.col_num);
//	sys_spt.bler		= (double)sys_spt.block_error/blocknum;
//	sys_spt.hard_ber	= (double)sys_spt.hard_error/(blocknum*sys_spt.col_num);
//	sys_spt.undet_ber	= (double)sys_spt.undetected_bit_error/(blocknum*sys_spt.col_num);
//	sys_spt.undet_bler	= (double)sys_spt.undetected_block_error/blocknum;

	printf("BSC_P:%1.4f blockCnt=%3d harderror =%5d  biterror =%3d  \n",\
		BSC_P,blocknum,sys_spt.hard_error, sys_spt.bit_error);


}
void UVLC_Ber_Cal(int blocknum)
{
	int i,j,k;
	int block_error_flag = 0;
	

	sys_spt.hard_error=0;
	sys_spt.bit_error=0;

	for(i=0;i<sys_spt.message_bits;i++)
	{
	if (sys_spt.source_bits[i]!=(sys_UVLC.received[i]) )
		{
			sys_spt.hard_error++;
		}
	}
	for(i=0;i<sys_spt.message_bits;i++)
	{
		sys_UVLC.decoded_out[i]=0;
	}
	for(k=0,i=0;i<128;i++)
	{
		for(j=0;j<sys_UVLC.code[sys_UVLC.decoded[i]][0];j++)
		{
			sys_UVLC.decoded_out[k++]=sys_UVLC.code[sys_UVLC.decoded[i]][j+1];
		}

	}
	for(i=0;i<sys_spt.message_bits;i++)
	{
		if (sys_spt.source_bits[i]!=sys_UVLC.decoded_out[i])  
		{
			sys_spt.bit_error++;
			
		}
	}
	    
   printf("BSC_P:%1.4f blockCnt=%3d  harderr=%3d biterror=%3d \n", BSC_P,blocknum,sys_spt.hard_error,sys_spt.bit_error);

}
void PrintResult(long blocknum){
	
	sys_spt.run_time = (sys_spt.stop_time - sys_spt.start_time)/CLOCKS_PER_SEC;
	sys_spt.fp_result= fopen(RESULT_FILE,"a+");
	fprintf(sys_spt.fp_result,"\n------------   SNR   =   %2.2f ----------\n",sys_spt.EbNo);
	fprintf(sys_spt.fp_result,"hard     error       =   %8d\n",sys_spt.hard_error);
	fprintf(sys_spt.fp_result,"decode   error       =   %8d\n",sys_spt.bit_error);
	fprintf(sys_spt.fp_result,"undet    error       =   %8d\n",sys_spt.undetected_bit_error);
	fprintf(sys_spt.fp_result,"block    error       =   %8d\n",sys_spt.block_error);	 
	fprintf(sys_spt.fp_result,"undetbl  error       =   %8d\n",sys_spt.undetected_block_error);	 
	fprintf(sys_spt.fp_result,"max      iters       =   %8d\n",sys_spt.max_iterations);
	fprintf(sys_spt.fp_result,"min      iters       =   %8d\n",sys_spt.min_iterations);
	fprintf(sys_spt.fp_result,"average  iters       =   %8.2f\n",sys_spt.aver_iterations/blocknum);
	fprintf(sys_spt.fp_result,"hard     BER         =   %1.4e\n",sys_spt.hard_ber);
	fprintf(sys_spt.fp_result,"decoding BER         =   %1.4e\n",sys_spt.ber);
	fprintf(sys_spt.fp_result,"block    BER         =   %1.4e\n",sys_spt.bler);	 	
	fprintf(sys_spt.fp_result,"undet    BER         =   %1.4e\n",sys_spt.undet_ber);
	fprintf(sys_spt.fp_result,"undetbl  BER         =   %1.4e\n",sys_spt.undet_bler);
//add for test iteration times needed
	for (int i=0;i<sys_spt.iter_times;i++)
	{
		fprintf(sys_spt.fp_result,"%d\t",sys_spt.error_iter[i]);	
	}
	fprintf(sys_spt.fp_result,"\n");
	for (i=0;i<sys_spt.iter_times;i++)
	{
		fprintf(sys_spt.fp_result,"%1.4e\t",sys_spt.ber_iter[i]);	
	}
	fprintf(sys_spt.fp_result,"\n");

	fprintf(sys_spt.fp_result,"Time ellapsed for %d Block processing is %d seconds!!\n", \
			blocknum, sys_spt.run_time );
	fclose(sys_spt.fp_result);
}