tudecsova.c

B3g_phase2_C语言_Matlab程序及说明

/*
* function [dec_out] = TuDecSova(dec_input, puncture, iter_num, alpha, sigma2, channel_type, g1, g2)
*
*	Inputs
*			dec_input:decoder input
*			puncture:puncture flag
*			iter_num:iterative number
*			alpha	:interleaver index
*
*	Output:
*			dec_out	:Turbo decoder output
*
*	Created July 15, 2002
*	Wang Dongming NCRL SEU
*/


#include<stdio.h>
#include <stdlib.h>
#include <string.h>
#include<math.h>
#include "mex.h"

#define N_MAX 17000
#define MAX_INTER_SIZE 32768
#define NU2_MAX 32
/* redefine HUGE */
#undef HUGE
#define HUGE 1e6

int SBG[NU2_MAX][2];
int SFG[NU2_MAX][2];

/* sova */
double delta[2][2][N_MAX]; 
int    yk[2][NU2_MAX];
/* retrieve the number of cells of the code */
double path_metric[N_MAX][NU2_MAX];
/* the metric difference of the two path */
double metric_diff[N_MAX][NU2_MAX];
/* the previous decition bit  */
int prev_bit[N_MAX][NU2_MAX];
int estimate_bit[N_MAX];
int mlstate[N_MAX];




int return_gi (int g, int position)
{
	return ( (g & (1 << (position)) ) != 0 );
}


int memory_length (int g)
{
	int i;
	
	for (i = 16; i >= 0; i--)
  		if (return_gi(g,i) != 0)
    		return i+1;
	return 0;
}


int dot (int x, int y)
{
	int i,l;
	int result;

	l = 32;
	y=y << 1;
	result = 0;
	for (i=0; i<l; i++)
		result = result + (return_gi(x,i) & return_gi(y,i));
	
	result = result % 2;

	return result;
}



int dot2 (int x, int y)
{
	int i,l ;
	int result;
	l = 32;   
	result = 0;

	for (i=0; i<l; i++)
		result = result + (return_gi(x,i) & return_gi(y,i));
	

	result = result % 2;
	return result;
}



int Sbcp(int m,int j,int g1)
{
	int nu;
	nu = memory_length(g1) - 1;

	return ((m >> 1) + (1 << (nu-1)) * ( (j + dot2(g1,m)) % 2 ));
}



int Sfcp(int m, int j, int g1)
{
	int nu;
	nu = memory_length(g1)-1;

	return  (m<<1)%(1<<nu) + ((j + dot(g1,m)) %2);
}



void interleave(double * block, int N, double *alpha)
{
	int j;
	double interleaved_block[MAX_INTER_SIZE];

	/* interleave to buffer */
	for (j=0; j<N; j++)
		interleaved_block[j] = block[(int)alpha[j]];

	
	for (j=0;j<N;j++)
		block[j]=interleaved_block[j];
}



void deinterleave(double * block, int N, double *alpha)
{
	int j;
	double interleaved_block[MAX_INTER_SIZE];

	/* deinterleave in temporary buffer */
	for (j=0; j<N; j++)
		interleaved_block[(int)alpha[j]] = block[j];

	/* copy buffer to final location */
	for (j=0; j<N; j++)
		block[j] = interleaved_block[j];
}



void sova_decoder (double * buffer_in, double * buffer_out, double *ploop,int dec_index,int N,int g1,int g2)
{
	/* loop variables */
	int i,j,k,m;

	/* number of cells of the code */
	int nu;
	/* the state number */
	int state_num;
	int state0,state1;

	/* the next state */
	int temp_state,bit;
	int win_wide = 30;
	
	/* scratch variables */
	double metric0,metric1;
	double temp;
	double llr;

	
	nu = memory_length (g1) - 1;
	state_num = 1 << nu;

	
	for (k=0; k<N; k++)
	for (m=0; m<state_num; m++)
		path_metric[k][m] = -(double)HUGE;

	path_metric[0][0] = 0.0;

	
	for (k=0; k<N; k++)
  	{
		delta[0][0][k] = -buffer_in[2*k] - buffer_in[2*k+1] - ploop[k];
		delta[0][1][k] = -buffer_in[2*k] + buffer_in[2*k+1] - ploop[k];
		delta[1][0][k] =  buffer_in[2*k] - buffer_in[2*k+1] + ploop[k];
		delta[1][1][k] =  buffer_in[2*k] + buffer_in[2*k+1] + ploop[k];
	}

	
	for (k=0; k<N; k++)
	{
		for (m=0; m<state_num; m++)
		{
			state0 = SBG[m][0];
			state1 = SBG[m][1];
			metric0 = (path_metric[k][state0] + delta[0][ yk[0][state0] ][k]);
			metric1 = (path_metric[k][state1] + delta[1][ yk[1][state0] ][k]);
			if (metric0 > metric1)
			{
				path_metric[k+1][m] = metric0;
				metric_diff[k+1][m] = 0.5*(metric0 - metric1);
				prev_bit[k+1][m] = 0;
			}
			else
			{
				path_metric[k+1][m] = metric1;
				metric_diff[k+1][m] = 0.5*(metric1 - metric0);
				prev_bit[k+1][m] = 1;
			}
		}
	}

	
	if (dec_index == 1)
	{
		mlstate[N] = 0;
	}
	else
	{
		temp = -HUGE;
		for (m=0; m<state_num; m++)
		{
			if (path_metric[N][m] > temp)
			{
				temp = path_metric[N][m];
				mlstate[N] = m;
			}
		}
	}

	
	for (k=N-1; k>=0; k--)
	{
		estimate_bit[k] = prev_bit[k+1][mlstate[k+1]];
		mlstate[k] = SBG[mlstate[k+1]][estimate_bit[k]];
	}

	
	for (k=0; k<N; k++)
	{
		llr = (double)HUGE;
		for (i=0; i<=win_wide; i++)
		{
			if ((k+i) < N)
			{
				/* the path that not selected */
				bit = 1 - estimate_bit[k+i];
				temp_state = SBG[ mlstate[k+i+1] ][bit];
				for (j=i-1; j>=0; j--)
				{
					bit = prev_bit[k+j+1][temp_state];
					temp_state = SBG[temp_state][bit];
				}
				if (bit != estimate_bit[k])
				{
					llr = min(llr,metric_diff[k+i+1][ mlstate[k+i+1] ]);
				}
			}
		}
		buffer_out[k] = (double)(2 * estimate_bit[k] -1) * (double)llr;
		
		/* get the extrinsic info. to the next decoder */
		ploop[k] = buffer_out[k] - buffer_in[2*k] - ploop[k];
	}
}


void turbo_decoder_sova(double *dec_input, double *hard_dec, int puncture, int iteration, double *alpha, int N, double sigma2, int channel_type, int poly_g1, int poly_g2)
{
	int i,j,iter;
	int nu, state_num;
	double Lc;
	
	/* The info to decoder one */
	double * pdecoder1_in;
	/* The info to decoder two */
	double * pdecoder2_in;
	/* The ouput info */
	double * pbuffer_out;
	/* The feed back info */
	double * ploop;
	
	/* Rayleigh Fading amplitude A=E(a)=0.8862 */
	double A = 0.8862;
	
	pdecoder1_in = (double *)malloc(N*2*sizeof(double));
	pdecoder2_in = (double *)malloc(N*2*sizeof(double));
	ploop = (double *)malloc(N*sizeof(double));
	pbuffer_out = (double *)malloc(N*sizeof(double));
	if ((pdecoder1_in == NULL)||(pdecoder2_in == NULL)||(ploop == NULL)||(pbuffer_out == NULL))
	{
		printf("not enough memery to alloc!\n");
		exit(0);
	}
	
	if (channel_type ==0)//gausion channel
	{
		Lc=2.0/sigma2;
	}
	else//rayleigh channel
	{
		Lc=A*2.0/sigma2;
	}
	
	/* retrieve the number of blocks of the encoder */
	nu = memory_length(poly_g1) - 1;
	state_num = (1<<nu);
	
	
	for (i=0; i<=1; i++)
  		for (j=0; j<state_num; j++)
			yk[i][j] = (dot(poly_g2,j)+((i+dot(poly_g1,j))%2) ) %2;

	
	for (i=0; i<(state_num); i++)
	{
		SBG[i][0] = Sbcp(i,0,poly_g1);
		SBG[i][1] = Sbcp(i,1,poly_g1);
		
		SFG[i][0] = Sfcp(i,0,poly_g1);
		SFG[i][1] = Sfcp(i,1,poly_g1);
	}
	
	
	if (puncture > 0)	/* unpunctured  rate = 1/3 */
	{
		for (i=0;i<N;i=i+1)
		{
			pdecoder1_in[2*i] = dec_input[i*3];
			pdecoder1_in[2*i+1] = dec_input[i*3+1];
			pdecoder2_in[2*i+1] = dec_input[i*3+2];
			pbuffer_out[i] = pdecoder1_in[2*i];
		}
	}
	else
	{
		for (i=0;i<N;i=i+1)
		{
			pdecoder1_in[2*i] = dec_input[i*2];
			pbuffer_out[i] = pdecoder1_in[2*i];
			if (i%2)
			{
				pdecoder1_in[2*i+1] = dec_input[2*i+1];
				pdecoder2_in[2*i+1] = 0.0;
			}
			else
			{
				pdecoder2_in[2*i+1] = dec_input[2*i+1];
				pdecoder1_in[2*i+1] = 0.0;
			}
		}
	}

	
	interleave(pbuffer_out, N, alpha);
	for (i=0;i<N;i=i+1)
	{
		pdecoder2_in[2*i] = pbuffer_out[i];
	}

	
	for (i=0;i<N;i++)
		ploop[i] = 0.0;
	
	for(iter=1;iter<=iteration;iter++)
	{
		sova_decoder(pdecoder1_in,pbuffer_out,ploop,1,N,poly_g1,poly_g2);
		
		
		interleave(ploop, N, alpha);
		
		
		sova_decoder(pdecoder2_in,pbuffer_out,ploop,2,N,poly_g1,poly_g2);
		
		/* deinterleave the feed back loop */
		deinterleave(ploop, N, alpha);
	}

	
	deinterleave(pbuffer_out, N, alpha);
	for (i=0;i<N;i++)
	{
		if (i<(N-nu))
		{
			if (pbuffer_out[i] < 0.0)
				hard_dec[i]=0;
			else
				hard_dec[i]=1;
		}
	}

	free(pdecoder1_in);
	free(pdecoder2_in);
	free(ploop);
	free(pbuffer_out);
}


/* The gateway routine. */
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
	/* input parameters */
	double *dec_input;
	double puncture;
	double iter_num;
	double *alpha;
	double sigma2;
	double channel_type;
	double poly_g1;
	double poly_g2;
	

	/* output parameters */
	double *hard_dec;
	
	int len, N, info_len;
	
	/* Check for the proper number of arguments */
	if ( nrhs != 8 )
		mexErrMsgTxt("Incorrect number of inputs. eight arguments");
	if ( nlhs != 1)
		mexErrMsgTxt("Incorrect number of outputs.one output");

	/* Get the input parameters */
	dec_input = mxGetPr(prhs[0]);
	puncture = mxGetScalar(prhs[1]);
	iter_num = mxGetScalar(prhs[2]);
	alpha = mxGetPr(prhs[3]);
	sigma2 = mxGetScalar(prhs[4]);
	channel_type = mxGetScalar(prhs[5]);
	poly_g1 = mxGetScalar(prhs[6]);
	poly_g2 = mxGetScalar(prhs[7]);
	
	len = mxGetN(prhs[0]);
	
	/* create a new array and set the output pointer to it */
	if (puncture > 0)
		N = (int) len/3;
	else
		N = (int) len/2;
	
	info_len = N - (memory_length((int)poly_g1) - 1);

	plhs[0]= mxCreateDoubleMatrix(1, info_len, mxREAL);
	
	hard_dec = mxGetPr(plhs[0]);

	turbo_decoder_sova(dec_input, hard_dec, (int)puncture, (int)iter_num, alpha, N, sigma2, (int)channel_type, (int)poly_g1, (int)poly_g2);
	
}