turbodecoder.c

信道编解码中的turbo编解码程序

* function reference:
* author:
*	Wang Jiaheng, 2004-12-23	created
*	Wang Jiaheng, 2004-02-04	revised
*******************************************************************************
*/
void RunMaxLogMap(int decode_len, double apriori_llr[], double gamma_mid[MAX_TURBO_STATE_NUM][2][MAX_TURBO_PACKET_SIZE], double aposteriori_llr[])
{
	int		k, m;
	double	temp_0;
	double	temp_1;
	double	metric_0;
	double	metric_1;


	/* compute the value of gamma */
	for (k = 0; k < decode_len; k++)
	{
		temp_0 = 0.375 * apriori_llr[k];
		temp_1 = -temp_0;

		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			/* previous state with +1 */
			turbo_gamma[m][0][k] = temp_0 + gamma_mid[m][0][k];
			/* previous state with -1 */
			turbo_gamma[m][1][k] = temp_1 + gamma_mid[m][1][k];
		}
	}

	/* for test 
	pf = fopen("gamma.txt", "w");
	for (k = 0; k < decode_len; k++)
	{
		fprintf(pf, "%-4d  |  ", k);
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			fprintf(pf, "%12f  |  ", turbo_gamma[m][0][k]);
			fprintf(pf, "%12f  |  ", turbo_gamma[m][1][k]);
		}
		fprintf(pf, "\n\n");
	}
	fclose(pf); */


	/* initialize the first alfa */
	turbo_alfa[0][0] = 0;
	for (m = 1; m < MAX_TURBO_STATE_NUM; m++)
	{
		turbo_alfa[m][0] = MAX_NEGATIVE;
	}

	for (k = 1; k < decode_len; k++)
	{
		/* compute the alfa at the kth time */
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			/* previous states with +1 */
			temp_0 = turbo_gamma[ turbo_prev_state[m][0] ][0][k - 1] + turbo_alfa[ turbo_prev_state[m][0] ][k - 1];
			/* previous states with -1 */
			temp_1 = turbo_gamma[ turbo_prev_state[m][1] ][1][k - 1] + turbo_alfa[ turbo_prev_state[m][1] ][k - 1];
            turbo_alfa[m][k]=log(exp(temp_0)+exp(temp_1));
			/* select the bigger one */
		/*	if ( temp_0 > temp_1 )
			{
				turbo_alfa[m][k] = temp_0;
			}
			else
			{
				turbo_alfa[m][k] = temp_1;
			}*/
		}
	}

	/* for test 
	pf = fopen("alfa.txt", "w");
	for (k = 0; k < decode_len; k++)
	{
		fprintf(pf, "%-4d  |  ", k);
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			fprintf(pf, "%12f  |  ", turbo_alfa[m][k]);
		}
		fprintf(pf, "\n\n");
	}
	fclose(pf);*/


	/* initialize the last beta */
	turbo_beta[0][decode_len - 1] = 0;
	for (m = 1; m < MAX_TURBO_STATE_NUM; m++)
	{
		turbo_beta[m][decode_len - 1] = MAX_NEGATIVE;
	}

	for (k = decode_len - 1; k > 0; k--)
	{
		/* compute the beta at the kth time */
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			/* posterior states with +1 */
			temp_0 = turbo_gamma[m][0][k] + turbo_beta[ turbo_post_state[m][0] ][k];
			/* posterior states with -1 */
			temp_1 = turbo_gamma[m][1][k] + turbo_beta[ turbo_post_state[m][1] ][k];
           turbo_beta[m][k - 1] =log(exp(temp_0)+exp(temp_1));
			/* select the bigger one */
		/*	if ( temp_0 > temp_1 )
			{
				turbo_beta[m][k - 1] = temp_0;
			}
			else
			{
				turbo_beta[m][k - 1] = temp_1;
			}*/
		}
	}

	/* for test 
	pf = fopen("beta.txt", "w");
	for (k = 0; k < decode_len; k++)
	{
		fprintf(pf, "%-4d  |  ", k);
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			fprintf(pf, "%12f  |  ", turbo_beta[m][k]);
		}
		fprintf(pf, "\n\n");
	}
	fclose(pf);*/


	/* compute the aposteriori LLR */
	for (k = 0; k < decode_len; k++)
	{
		/* states with +1 */
		temp_0 = exp(turbo_alfa[0][k] + turbo_gamma[0][0][k] + turbo_beta[ turbo_post_state[0][0] ][k]);
		/* states with -1 */
		temp_1 = exp(turbo_alfa[0][k] + turbo_gamma[0][1][k] + turbo_beta[ turbo_post_state[0][1] ][k]);
		
		for (m = 1; m < MAX_TURBO_STATE_NUM; m++)
		{
			/* states with +1 */
			temp_0 += exp(turbo_alfa[m][k] + turbo_gamma[m][0][k] + turbo_beta[ turbo_post_state[m][0] ][k]);
			/* states with -1 */
			temp_1 += exp(turbo_alfa[m][k] + turbo_gamma[m][1][k] + turbo_beta[ turbo_post_state[m][1] ][k]);
            
			/* select the biggest one */
		/*	if ( temp_0 > metric_0 )
			{
				metric_0 = temp_0;
			}

			/* select the biggest one */
		//	if ( temp_1 > metric_1 )
		//	{
		//		metric_1 = temp_1;
		//	}
		}

		aposteriori_llr[k] =log(temp_0)-log(temp_1);//metric_0 - metric_1;
	}

	/* for test 
	pf = fopen("llr.txt", "w");
	for (k = 0; k < decode_len; k++)
	{
		fprintf(pf, "%-4d  |  ", k);
		fprintf(pf, "%12f\n\n", turbo_aposteriori_llr[k]);
	}
	fclose(pf);*/
}


/*
******************************************************************************* 
* description:
*	run the turbo decoder with a encoded bit sequence input and a packet output
* input: 
*	sym_in: input symbols sequence
*	len_in: length of input symbols sequence
*	encode_rate: the encoding rate (2, 3, 4 or 5) 
*	ch_rely: channel reliability (2a/noise_var)
*	iterate_num: number of iterating decoding
* output:
*	bit_out: encoded output bit sequence
*	*p_len_out: length of the output bit sequence
*	if ok, return 1; otherwise retrun 0.
* function reference:
*	EV-DO specification
* author:
*	Wang Jiaheng, 2004-12-23	created
*	Wang Jiaheng, 2004-02-04	revised
*******************************************************************************
*/
int RunTurboDecoder(double sym_in[], int len_in, int encode_rate, int iterate_num, int bit_out[], int *p_len_out)
{
	int		k, m;
	int		packet_size;
	int		decode_len;
	int		valid_len;
	double	sum_0, sum_1;


	/* compute the packet size and check it */
	packet_size = len_in / encode_rate;
	switch ( packet_size )
	{
	case 256:
		break;
	case 512:
		break;
	case 1024:
		break;
	case 2048:
		break;
	case 3072:
		break;
	case 4096:
		break;
	default:
		return 0;
	}

	/* compute other lengths */
	decode_len = packet_size - 3;
	valid_len  = packet_size - 6;
	if ( valid_len != turbo_interleaved_len )
	{
		return 0;
	}


	/* rebuild the received symbols for decoder 1 and 2 */
	if ( !RebuildSymbol(sym_in, valid_len, decode_len, encode_rate) )
	{
		return 0;
	}


	/* compute the middle values for gamma in the Max-log-Map algorithm */
	for (k = 0; k < decode_len; k++)
	{
		/* for the constituent encoder 1 */
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			sum_0  = turbo_recv_sym_1[0][k] * turbo_trellis_output[m][0][0];
			sum_0 += turbo_recv_sym_1[1][k] * turbo_trellis_output[m][0][1];
			sum_0 += turbo_recv_sym_1[2][k] * turbo_trellis_output[m][0][2];

			sum_1  = turbo_recv_sym_1[0][k] * turbo_trellis_output[m][1][0];
			sum_1 += turbo_recv_sym_1[1][k] * turbo_trellis_output[m][1][1];
			sum_1 += turbo_recv_sym_1[2][k] * turbo_trellis_output[m][1][2];

			turbo_gamma_mid_1[m][0][k] = 0.5 * sum_0;
			turbo_gamma_mid_1[m][1][k] = 0.5 * sum_1;
		}

		/* for the constituent encoder 2 */
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			sum_0  = turbo_recv_sym_2[0][k] * turbo_trellis_output[m][0][0];
			sum_0 += turbo_recv_sym_2[1][k] * turbo_trellis_output[m][0][1];
			sum_0 += turbo_recv_sym_2[2][k] * turbo_trellis_output[m][0][2];

			sum_1  = turbo_recv_sym_2[0][k] * turbo_trellis_output[m][1][0];
			sum_1 += turbo_recv_sym_2[1][k] * turbo_trellis_output[m][1][1];
			sum_1 += turbo_recv_sym_2[2][k] * turbo_trellis_output[m][1][2];

			turbo_gamma_mid_2[m][0][k] = 0.5 * sum_0;
			turbo_gamma_mid_2[m][1][k] = 0.5 * sum_1;
		}
	}


	/* initialize the a-priori LLR sequence to 0 */
	memset(turbo_apriori_llr_1, 0, MAX_TURBO_PACKET_SIZE * sizeof(double));
	memset(turbo_apriori_llr_2, 0, MAX_TURBO_PACKET_SIZE * sizeof(double));


	/* start iterative decoding */
	for (m = 0; m < iterate_num; m++)
	{
		/* run the decoder 1 */
		RunMaxLogMap(decode_len, turbo_apriori_llr_1,  turbo_gamma_mid_1, turbo_aposteriori_llr_1);

		/* compute the extrinsic information as the a-priori LLR for the edcoder 2 */
		for (k = 0; k < decode_len; k++)
		{
			turbo_apriori_llr_1[k] = turbo_aposteriori_llr_1[k] - 0.75 * turbo_apriori_llr_1[k] - turbo_recv_sym_1[0][k];
		}

		/* interleave the extrinsic LLR */
		for (k = 0; k < valid_len; k++)
		{
			turbo_apriori_llr_2[k] = turbo_apriori_llr_1[ turbo_interleaved_addr[k] ];
		}
		for (k = valid_len; k < decode_len; k++)
		{
			turbo_apriori_llr_2[k] = 0;
		}

		/* run the decoder 2 */
		RunMaxLogMap(decode_len, turbo_apriori_llr_2, turbo_gamma_mid_2, turbo_aposteriori_llr_2);

		/* compute the extrinsic LLR as the a-priori LLR for the edcoder 1 */
		for (k = 0; k < decode_len; k++)
		{
			turbo_apriori_llr_2[k] = turbo_aposteriori_llr_2[k] - 0.75 * turbo_apriori_llr_2[k] - turbo_recv_sym_2[0][k];
		}

		/* deinterleave the extrinsic LLR */
		for (k = 0; k < valid_len; k++)
		{
			turbo_apriori_llr_1[ turbo_interleaved_addr[k] ] = turbo_apriori_llr_2[k];
		}
		for (k = valid_len; k < decode_len; k++)
		{
			turbo_apriori_llr_1[k] = 0;
		}
	}

	/* deinterleave the a-posteriori LLRs generated by the decoder 2 */
	for (k = 0; k < valid_len; k++)
	{
		turbo_aposteriori_llr_1[ turbo_interleaved_addr[k] ] = turbo_aposteriori_llr_2[k];
	}


	/* make a hard decision according to the a-posteriori LLR sequence */
	for (k = 0; k < valid_len; k++)
	{
		if ( turbo_aposteriori_llr_1[k] >= 0 )
		{
			bit_out[k] = 0;
		}
		else
		{
			bit_out[k] = 1;
		}
	}

	for (k = valid_len; k < packet_size; k++)
	{
		bit_out[k] = 0;
	}

	*p_len_out = packet_size;

	return 1;
}