turbodecoder.c

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

* function reference:
* author:
*	Wang Jiaheng, 2004-12-11	created
*	Wang Jiaheng, 2004-02-04	revised
*******************************************************************************
*/
void RunMap(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	numerator;
	double	denominator;
	double	max_alfa;
	double	max_beta;
	double	norm_factor;


	/* compute the value of gamma */
	for (k = 0; k < decode_len; k++)
	{
		temp_0 = exp(0.475 * apriori_llr[k]);
		temp_1 = 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] = 1;
	for (m = 1; m < MAX_TURBO_STATE_NUM; m++)
	{
		turbo_alfa[m][0] = 0;
	}

	for (k = 1; k < decode_len; k++)
	{
		/* compute the alfa at the kth time */
		max_alfa = 1;
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			/* previous states with +1 */
			turbo_alfa[m][k]  = turbo_gamma[ turbo_prev_state[m][0] ][0][k - 1] * turbo_alfa[ turbo_prev_state[m][0] ][k - 1];
			/* previous states with -1 */
			turbo_alfa[m][k] += turbo_gamma[ turbo_prev_state[m][1] ][1][k - 1] * turbo_alfa[ turbo_prev_state[m][1] ][k - 1];

			if ( turbo_alfa[m][k] > max_alfa )
			{
				max_alfa = turbo_alfa[m][k];
			}
		}

		/* normalize the alfa at the current time */
		norm_factor = 1 / max_alfa;
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			turbo_alfa[m][k] = turbo_alfa[m][k] * norm_factor;
		}
	}

	/* 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] = 1;
	for (m = 1; m < MAX_TURBO_STATE_NUM; m++)
	{
		turbo_beta[m][decode_len - 1] = 0;
	}

	for (k = decode_len - 1; k > 0; k--)
	{
		/* compute the beta at the kth time */
		max_beta = 1;
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			/* posterior states with +1 */
			turbo_beta[m][k - 1]  = turbo_gamma[m][0][k] * turbo_beta[ turbo_post_state[m][0] ][k];
			/* posterior states with -1 */
			turbo_beta[m][k - 1] += turbo_gamma[m][1][k] * turbo_beta[ turbo_post_state[m][1] ][k];

			if ( turbo_beta[m][k - 1] > max_beta )
			{
				max_beta = turbo_beta[m][k - 1];
			}
		}

		/* normalize the beta */
		norm_factor = 1 / max_beta;
		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			turbo_beta[m][k - 1] = turbo_beta[m][k - 1] * norm_factor;
		}
	}

	/* 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++)
	{
		numerator   = 0;
		denominator = 0;

		for (m = 0; m < MAX_TURBO_STATE_NUM; m++)
		{
			/* states with +1 */
			numerator   += turbo_alfa[m][k] * turbo_gamma[m][0][k] * turbo_beta[ turbo_post_state[m][0] ][k];
			/* states with -1 */
			denominator += turbo_alfa[m][k] * turbo_gamma[m][1][k] * turbo_beta[ turbo_post_state[m][1] ][k];
		}

		aposteriori_llr[k] = log(numerator / denominator);
	}

	/* 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-11	created
*	Wang Jiaheng, 2005-02-04	revised
*******************************************************************************
*/
int RunTurboDecoder(double sym_in[], int len_in, int encode_rate, double ch_rely, 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;
	double	half_lc;


	/* compute and check the packet size */
	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 exponential values for gamma in the MAP algorithm */
	half_lc = ch_rely / 2;
	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] = exp(half_lc * sum_0);
			turbo_gamma_mid_1[m][1][k] = exp(half_lc * 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] = exp(half_lc * sum_0);
			turbo_gamma_mid_2[m][1][k] = exp(half_lc * 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 */
		RunMap(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.95 * turbo_apriori_llr_1[k] - ch_rely * 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 */
		RunMap(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.95 * turbo_apriori_llr_2[k] - ch_rely * 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;
}