turbo_code_log_map.c

myturbo 简单使用的turbo系统方针程序 对于理解turbo码的编解码很有帮助

	free(tempindex);
}



void randominterleaver_int(int *data_unintlvr, int *interleaverddata,
						   int length)
{
	int i;
	int *index_random =  index_randomintlvr ;

	for (i=0; i<length; i++)
	{
		*(interleaverddata+i) = *(data_unintlvr+ (*(index_random+i)));
	}
}

void randominterleaver_float(float *data_unintlvr, float *interleaverddata, int length)
{
	int i;
	int *index_random =  index_randomintlvr ;

	for (i=0; i<length; i++)
	{
		*(interleaverddata+i) = *(data_unintlvr+ (*(index_random+i)));
	}
}


void random_deinterlvr_int(int *data_unintlvr, int *interleaverddata, int length)
{
	int i;
	int *index_random =  index_randomintlvr;

	for (i=0; i<length; i++)
	{
		*(data_unintlvr+(*(index_random+i))) = *(interleaverddata+i);
	}
}


void random_deinterlvr_float(float *data_unintlvr, float *interleaverddata, int length)
{
	int i;
	int *index_random = index_randomintlvr ;

	for (i=0; i<length; i++)
	{
		*(data_unintlvr+(*(index_random+i))) = *(interleaverddata+i);
	}
}


void encoderm_turbo(int *source, int *send_turbo, int len_info)
{
	int i;
	int len_total = len_info + M_num_reg;

	int *rsc1, *rsc2;
	
	int *input2;
	int *source2;

	if ((rsc1=(int *)malloc(2*len_total*sizeof(int)))==NULL)
	{
	  printf("\n fail to allocate memory of rsc1 \n");
	  exit(1);  
	}
	if ((rsc2=(int *)malloc(2*len_total*sizeof(int)))==NULL)
	{
	  printf("\n fail to allocate memory of rsc2 \n");
	  exit(1);  
	}

	if ((input2=(int *)malloc(len_total*sizeof(int)))==NULL)
	{
	  printf("\n fail to allocate memory of input2 \n");
	  exit(1);  
	}

	/* the first RSC encoder */
	rsc_encode(source, rsc1, 1, len_info);

      if((source2=(int *)malloc(len_total*sizeof(int)))==NULL)
      	{
	  printf("\n fail to allocate memory of input2 \n");
	  exit(1);  
	}
      		   
	  for(i=0;i<len_total;i++)
	   *(source2+i)=*(rsc1+2*i);
	   
	randominterleaver_int(source2, input2,  len_total);

	/* the second RSC encoder */
	rsc_encode(input2, rsc2, 0, len_total);

	for (i=0; i<len_total; i++)
	{
		*(send_turbo+2*i) = *(rsc1+2*i) *2 - 1;
            if(i%2)
		*(send_turbo+2*i+1) = *(rsc1+2*i+1) *2 - 1;
	     else
		*(send_turbo+2*i+1) = *(rsc2+2*i+1) *2 - 1;
	}
	
		
	free(rsc1);
	free(rsc2);
	free(input2);
}



void rsc_encode(int *source, int *rsc, int terminated, int len_source)
{
	int i, j;

	int *state;
	int dk_turbo, ak_turbo, outbit;

	int len_total;
	int len_info;

	if ((state=(int *)malloc(M_num_reg*sizeof(int)))==NULL)
	{
	  printf("\n fail to allocate memory of state \n");
	  exit(1);  
	}
	for (i=0; i<M_num_reg; i++)
	{
		*(state+i) = 0;
	}

	   if(terminated)
	   {
	     len_info=len_source;
	     len_total = len_source+M_num_reg;

	   }
	   else
	   	{len_total=len_source;
	   	len_info=len_total-M_num_reg;
		}
			
	
       
	for (i=0; i<len_total; i++)			/* encoding bit by bit */
	{
		if (!terminated || (terminated && i<len_info))	/* for information bits */
		{
			dk_turbo = *(source+i);
		}
		else											/* terminate the trellis */
		{
			if (terminated && i>=len_info)
			{
				dk_turbo = 0;
				for (j=1; j<turbo_g.K_num_col; j++)
				{
					dk_turbo = dk_turbo + (*(turbo_g.g_matrix+j)) * (*(state+j-1));
				}
				dk_turbo = dk_turbo%2;
			}
		}

		ak_turbo = *(turbo_g.g_matrix+0) * dk_turbo;
		for (j=1; j<turbo_g.K_num_col; j++)
		{
			ak_turbo = ak_turbo + (*(turbo_g.g_matrix+j))*(*(state+j-1));
		}

		ak_turbo = ak_turbo%2;

		outbit = encode_bit(ak_turbo, state);

		*(rsc+2*i) = dk_turbo;
		*(rsc+2*i+1) = outbit;
	}				/* end of encoding bit by bit */

	free(state);
}




int encode_bit(int inbit, int *stat)
{
	int j;
	int output;

	output = (*(turbo_g.g_matrix+turbo_g.K_num_col+0)) * inbit;

	for (j=1; j<turbo_g.K_num_col; j++)
	{
		output = (output + (*(turbo_g.g_matrix+turbo_g.K_num_col+j)) * (*(stat+j-1)))%2;
	}

	for (j=turbo_g.K_num_col-2; j>0; j--)
	{
		*(stat+j)=*(stat+j-1);
	}

	*(stat+0) = inbit;

	return output;
}




void demultiplex(float *rec_turbo, int len_info, float *yk_turbo)
{
	int i;

	int len_total = len_info+M_num_reg;

	float *info2, *inted_info2;

	if ((info2=(float *)malloc(len_total*sizeof(float)))==NULL)
	{
	  printf("\n fail to allocate memory of info2 \n");
	  exit(1);  
	}
	if ((inted_info2=(float *)malloc(len_total*sizeof(float)))==NULL)
	{
	  printf("\n fail to allocate memory of inted_info2 \n");
	  exit(1);  
	}

       for(i=0;i<2*len_total;i++)
        {*(yk_turbo+2*i+1)=0;}

	
	/* for information bits */
	for(i=0; i<len_total; i++)
	{    
		*(info2+i) = *(yk_turbo+2*i) = *(rec_turbo+2*i);
                      if(i%2)
		*(yk_turbo+2*i+1) = *(rec_turbo+2*i+1);
			  else
		*(yk_turbo+2*len_total+2*i+1) = *(rec_turbo+2*i+1);
	}

	randominterleaver_float(info2, inted_info2,  len_total);

	for (i=0; i<len_total; i++)
	{
		*(yk_turbo+2*len_total+2*i) = *(inted_info2+i);
	}
         
	
	free(info2);
	free(inted_info2);
}


void Log_MAP_decoder(float *recs_turbo, float *La_turbo, int terminated, float *LLR_all_turbo, int len_total)
{
	int i, j;

	float *alpha_Log, *beta_Log, *gama_Log;

	float *tempmax;
	float *temp0, *temp1;
	float tempx, tempy;

	if ((alpha_Log=(float *)malloc(n_states*(len_total+1)*sizeof(float)))==NULL)
	{
	  printf("\n fail to allocate memory of alpha_Log \n");
	  exit(1);  
	}
	if ((beta_Log=(float *)malloc(n_states*(len_total+1)*sizeof(float)))==NULL)
	{
	  printf("\n fail to allocate memory of beta_Log \n");
	  exit(1);  
	}
	if ((gama_Log=(float *)malloc(n_states*len_total*2*sizeof(float)))==NULL)
	{
	  printf("\n fail to allocate memory of gama_Log \n");
	  exit(1);  
	}

	if ((tempmax=(float *)malloc((len_total+1)*sizeof(float)))==NULL)
	{
	  printf("\n fail to allocate memory of tempmax \n");
	  exit(1);  
	}

	if ((temp0=(float *)malloc(n_states*sizeof(float)))==NULL)
	{
	  printf("\n fail to allocate memory of temp0 \n");
	  exit(1);  
	}
	if ((temp1=(float *)malloc(n_states*sizeof(float)))==NULL)
	{
	  printf("\n fail to allocate memory of temp1 \n");
	  exit(1);  
	}
	
	/*===== Initialization of alpha_Log and beta_Log =====*/
	*(alpha_Log+0) = 0;
	*(beta_Log+len_total) = 0;

	for (i=1; i<n_states; i++)
	{
		*(alpha_Log+i*(len_total+1)+0) = (float)-INFTY;

		if (terminated)
		{
			*(beta_Log+i*(len_total+1)+len_total) = (float)-INFTY;
		}
		else
		{
			*(beta_Log+i*(len_total+1)+len_total) = 0;
		}
	}
	/*========compute Gama_Log========*/
	for (i=0; i<len_total; i++)			/* 0--len_total-1 代表 1--len_total */
	{
		for (j=0; j<n_states; j++)		/* j->k */
		{
			*(gama_Log+j*len_total*2+i*2+0) 
				= -*(recs_turbo+2*i) + *(recs_turbo+2*i+1)*(*(turbo_trellis.mx_nextout+j*4+1)) - *(La_turbo+i)/2;
			*(gama_Log+j*len_total*2+i*2+1)
				= *(recs_turbo+2*i) + *(recs_turbo+2*i+1)*(*(turbo_trellis.mx_nextout+j*4+3)) + *(La_turbo+i)/2;
		}
	}

	/*========Trace forward, compute Alpha_Log========*/
	for (i=1; i<len_total+1; i++)
	{
		for (j=0; j<n_states; j++)		/* 以j为中心*/
		{
			tempx = *(gama_Log+(*(turbo_trellis.mx_laststat+j*2+0))*len_total*2+(i-1)*2+0)
					+ *(alpha_Log+(*(turbo_trellis.mx_laststat+j*2+0))*(len_total+1)+i-1);
			tempy = *(gama_Log+(*(turbo_trellis.mx_laststat+j*2+1))*len_total*2+(i-1)*2+1)
					+ *(alpha_Log+(*(turbo_trellis.mx_laststat+j*2+1))*(len_total+1)+i-1);
			*(alpha_Log+j*(len_total+1)+i) = E_algorithm(tempx, tempy);
		}

		for (j=0; j<n_states; j++)
		{
			if (*(tempmax+i) < *(alpha_Log+j*(len_total+1)+i))
			{
				*(tempmax+i) = *(alpha_Log+j*(len_total+1)+i);
			}
		}

		for(j=0; j<n_states; j++)
		{
			*(alpha_Log+j*(len_total+1)+i) = *(alpha_Log+j*(len_total+1)+i) - *(tempmax+i);
		}

	}

	/*========Trace backward, compute Beta_Log========*/
	for (i=len_total-1; i>=0; i--)
	{
		for (j=0; j<n_states; j++)		/*j为中心*/
		{
			tempx = *(gama_Log+j*len_total*2+i*2+0) 
					+ *(beta_Log+(*(turbo_trellis.mx_nextstat+j*2+0))*(len_total+1)+i+1);
			tempy = *(gama_Log+j*len_total*2+i*2+1) 
					+ *(beta_Log+(*(turbo_trellis.mx_nextstat+j*2+1))*(len_total+1)+i+1);

			*(beta_Log+j*(len_total+1)+i) = E_algorithm(tempx, tempy);
		}

		for (j=0; j<n_states; j++)
		{
			*(beta_Log+j*(len_total+1)+i) = *(beta_Log+j*(len_total+1)+i) - *(tempmax+i+1);
		}
	}

	/*===Compute the soft output,log-likelihood ratio of symbols in the frame===*/
	for (i=0; i<len_total; i++)	
	{
		for (j=0; j<n_states; j++)
		{
			*(temp0+j) = *(gama_Log+(*(turbo_trellis.mx_laststat+j*2+0))*len_total*2+i*2+0) 
				+ *(alpha_Log+*(turbo_trellis.mx_laststat+j*2+0)*(len_total+1)+i)
				+ *(beta_Log+ j*(len_total+1)+i+1);

			*(temp1+j) = *(gama_Log+(*(turbo_trellis.mx_laststat+j*2+1))*len_total*2+i*2+1) 
				+ *(alpha_Log+*(turbo_trellis.mx_laststat+j*2+1)*(len_total+1)+i)
				+ *(beta_Log+j*(len_total+1)+i+1);
		}

		*(LLR_all_turbo+i) = E_algorithm_seq(temp1, n_states) - E_algorithm_seq(temp0, n_states);
	}

	free(alpha_Log);
	free(beta_Log);
	free(gama_Log);
	free(tempmax);
	free(temp0);
	free(temp1);
}


float E_algorithm(float x, float y) // log(exp(x) + exp(y)) = max(x,y)+f(-|y-x|) 
{
	float temp = (y-x)>0? (y-x):(x-y);
	int i;

	if (temp>=4.3758)
	{
		temp = 0;
	}
	else
	{
		for (i=0; i<16 && temp>=lookup_index_Log_MAP[i]; i++)
		{
			;
		}
		temp = (float)lookup_table_Log_MAP[i-1];
	}
	
	return ( (x>y?x:y) + temp );
}


float E_algorithm_seq(float *data_seq, int length)
{
	int i;
	float temp;
	temp = E_algorithm(*(data_seq+0), *(data_seq+1));
	for (i=2; i<length; i++)
	{
		temp = E_algorithm(temp, *(data_seq+i));
	}
	return temp;
}


void gen_source(int *data, int length)
{
	double temp;
	int i;

	for (i=0; i<length; i++)
	{
		temp = (double)rand()/RAND_MAX;
		//temp = random_turbo();
		if (temp <= 0.5)
		{
			*(data+i) = 0;
		}
		else 
		{
			*(data+i) = 1;
		}
	}

}



void AWGN(float *send, float *r, float sigma, int totallength)
{
	int i;
	double *noise = (double *)malloc(sizeof(double)*totallength);
	double seed =  3.0 - (double)((rand() & RAND_MAX)/(double)RAND_MAX)/10e6;
	mgrns(0,sigma,seed,totallength,noise);
	for(i=0; i<totallength; i++)
	{
		*(r+i) = (float)( *(send+i) + *(noise+i) );
	}
	free(noise);
}


void mgrns(double mean,double sigma,double seed,int n,double *a)
{ int i,k,m;
    double s,w,v,t;
    s=65536.0; w=2053.0; v=13849.0;
    for (k=0; k<=n-1; k++)
	{
		t=0.0;
		for (i=1; i<=12; i++)
        { 
			seed=seed*w+v; m=(int)(seed/s);
            seed=seed-m*s; t=t+(seed)/s;
        }/*按照中心极限定理产生服从高斯分布的随机数*/
        *(a+k)=mean+sigma*(t-6.0);
    }
    return;
}