codec.cpp

Turbo码性能仿真源程序信道编码方面的

#include <math.h>
#include <string.h>
#include "Global.h"
//#include "System.h"
#include <iostream.h>
#include "CoDec.h"

extern IO			IO;

/******************************************************************************************
function : Initialization and parameters check of the turbo-SPC CODEC
input    :	Nuser			-- number of user
			Datalen			-- length of data
			Codewordlen		-- length of codeword
			Enum			-- number of recursive bit
			Fnum			-- number of non-recursive bit
			Nestdim			-- number of component code
			Spctrellislen	-- length of the trellis
			Nodeorder		-- node order of the trellis
			Spcnpoly		-- numerator polynome -- less significate bit represents 2^0
			Spcdpoly 		-- denominator polynome 
******************************************************************************************/
void SpcCoDec::Init( int Nuser, int Datalen, int Pun, int Enum, int Fnum, int Nestdim, int Spctrellislen, int Nodeorder, int Spcnpoly, int Spcdpoly, int Spcinternum )
{
	NUser					= Nuser;
	DataLen					= Datalen;	
	Puncture				= Pun;
	NestDim					= Nestdim;
	EfoldNum				= Enum;
	FfoldNum				= Fnum;
	FoldNum					= EfoldNum + FfoldNum;
	TrellisLen				= Spctrellislen;
	NodeOrder				= Nodeorder;
	NodeNum					= 1 << Nodeorder;
	Npoly					= Spcnpoly;
	Dpoly					= Spcdpoly;
	CodewordLen				= DataLen + NestDim * ( TrellisLen + 2 * NodeOrder );
//	PuncturedCodewordLen	= DataLen + NestDim * ( TrellisLen + 2 * NodeOrder ) * ( NestDim - Puncture ) / NestDim;
	PuncturedCodewordLen	= DataLen + ( TrellisLen + 2 * NodeOrder ) * ( NestDim - Puncture );
	SpcInterNum				= Spcinternum;
	Rate					= (float) DataLen / (float) PuncturedCodewordLen;

	SpcCoDec::Check( );
	ConvSpc.Init( Nodeorder, TrellisLen+NodeOrder, Npoly, Dpoly );// Initialization of the convolutional SPC

	//set up the interleaver
	if( FfoldNum == 0 ) EInterleaver( );
	else FInterleaver( );
}

/******************************************************************************************
function : interleaver of the code with only recrsive component
******************************************************************************************/
void SpcCoDec::EInterleaver( )
{
	int i, j;
	
	for( j=0; j<DataLen; j++ ) ShaffleRule[0][j] = j;
	for( i=0; i<NestDim; i++ ) 
	{
		for( j=0; j<DataLen; j++ ) ShaffleRule[i][j] = j;
		Scramble( &ShaffleRule[i][0], DataLen, SpcInterNum );
	}
}
/******************************************************************************************
function : interleaver of the code with both recrsive and non-recursive component
******************************************************************************************/
void SpcCoDec::FInterleaver( )
{
	int i, j, k, m, n, kk;
	int Length		= TrellisLen;
	int Elength		= ( DataLen * EfoldNum ) / FoldNum;
	int Flength		= DataLen - Elength;
	int SegLength	= FoldNum / NestDim;
	int a			= Length * SegLength;
	int temp[DATALEN];

	for( i=0; i<NestDim; i++)
	{
		for( n=0, k=0; k<FoldNum; k+=SegLength ) for( j=0; j<Length; j++ )
		{
			m = j * FoldNum + k;
			for( kk=k; kk<k+SegLength; kk++ ) temp[n++] = m++;
		}
		CyclicShift( temp, DataLen, i*a );
		Scramble( temp, Elength, SpcInterNum );
		Scramble( &temp[Elength], Flength, SpcInterNum );
		for( n=0, k=0; k<FoldNum; k+=SegLength ) for( j=0; j<Length; j++ )
		{
			m = j * FoldNum + k;
			for( kk=k; kk<k+SegLength; kk++ ) ShaffleRule[i][m++] = temp[n++];
		}
	}
}
/******************************************************************************************
function : parameters check of the turbo-SPC
******************************************************************************************/
void SpcCoDec::Check( )
{
	if( DataLen % FoldNum != 0 ) IO.ErrorMessage( "turbo-SPC: DataLen % FoldNum != 0" );
	if( EfoldNum == 0 ) IO.ErrorMessage( "EfoldNum=0 ==> SIMPLE parity check -- use another prog." );
	if( FfoldNum != 0 )
	{
		if( FoldNum % NestDim != 0 ) IO.ErrorMessage( "turbo-SPC: FoldNum % NestDim != 0" ); // NOTE: there is something wrong here
		if( DataLen % NestDim != 0 ) IO.ErrorMessage( "turbo-SPC: DataLen % NestDim != 0" );
	}
}

/******************************************************************************************
function : encoder of turbo-SPC
input    : NUser, _RawData[NUser][]
output   : *Codeword
******************************************************************************************/
void SpcCoDec::Encoder( int Nuser, int *Codeword, int *RawData )
{
	int i, j, k, m, p;
	int state;
	int t = DataLen + NestDim * TrellisLen;

	for( i=0; i<DataLen; i++ ) Codeword[i] = RawData[i];
	for( p=DataLen, i=0; i<NestDim; i++ )
	{
		state = 0;
		for( m=0, j=0; j<TrellisLen; j++ )
		{
			DataIn[j] = RawData[ShaffleRule[i][m++]];
			for( k=1; k<EfoldNum; k++ ) DataIn[j] ^= RawData[ShaffleRule[i][m++]];
			m += FfoldNum;
		}
		ConvSpc.Encoder( &DataIn[0], &ParityOut[0], &state );
		for( m=EfoldNum, j=0; j<TrellisLen; j++ )
		{
			Codeword[p] = ParityOut[j];
			for( k=EfoldNum; k<FoldNum; k++ ) Codeword[p] ^= RawData[ShaffleRule[i][m++]];
			m += EfoldNum;
			p++;
		}
		for( j=0; j<2*NodeOrder; j++ ) Codeword[t++] = ConvSpc.TermTable[state][j];
	}
}
/******************************************************************************************
function : decoder of turbo-SPC
input    : NUser, *Codeword -- LR value
output   : *Codeword
******************************************************************************************/
void SpcCoDec::Decoder( int Nuser, float *Codeword, char *Flag )
{
	int i, j, p=DataLen+NestDim*TrellisLen;

	for( i=0; i<NestDim; i++ )
	{
		BorderReSet( );
		for( j=0; j<DataLen; j++ ) Zin[j] = RClip3( Codeword[ShaffleRule[i][j]] / ExInf[Nuser][i][j] );
		MapDecoding( &Codeword[DataLen+(i*TrellisLen)], &Codeword[p+(2*i*NodeOrder)], Nuser, i, Flag );// the major decoding step 
		for( j=0; j<DataLen; j++ ) Codeword[ShaffleRule[i][j]] = Zin[j];
	}
}

void SpcCoDec::BorderReSet( )
{
	int i;

	for( i=1; i<NodeNum; i++ )
		Border[0][i] = Border[1][i] = 0;
	Border[0][0] = Border[1][0] = 1;
}

/******************************************************************************************
function : decoder of convolutional-SPC
******************************************************************************************/
void SpcCoDec::MapDecoding( float *Parity, float *Tail, int Nuser, int dim, char *Flag )
{
	Apriori( Parity, Tail );
	ConvSpc.BCJR( &Border[0][0], &Border[1][0] );
	Aposteriori( Nuser, dim, Parity, Tail, Flag );
}
/******************************************************************************************
function : handle the information output from the BCJR
******************************************************************************************/
void SpcCoDec::Aposteriori( int Nuser, int dim, float *Parity, float *Tail, char *Flag )
{
	int i, m;

	for( i=0; i<DataLen; i++ ) tmpbf[i] = Zin[i];
	for( i=0; i<TrellisLen; i++ )
	{
		Eparity[i] = RClip3( RClip3((ConvSpc.aPP[i][0]+ConvSpc.aPP[i][1])/(ConvSpc.aPP[i][2]+ConvSpc.aPP[i][3]))/Eparity[i] );
		Fparity[i] = RClip3( RClip3((ConvSpc.aPP[i][0]+ConvSpc.aPP[i][2])/(ConvSpc.aPP[i][1]+ConvSpc.aPP[i][3]))/Fparity[i] );
	}

	m = TrellisLen;
	for( i=0; i<NodeOrder; i++ )
	{
		Tail[2*i] = RClip3( RClip3((ConvSpc.aPP[m][0]+ConvSpc.aPP[m][1])/(ConvSpc.aPP[m][2]+ConvSpc.aPP[m][3]))/Tail[2*i] );
		Tail[2*i+1] = RClip3( RClip3((ConvSpc.aPP[m][0]+ConvSpc.aPP[m][2])/(ConvSpc.aPP[m][1]+ConvSpc.aPP[m][3]))/Tail[2*i+1] );
		m++;
	}

	for( m=0, i=0; i<TrellisLen; i++, m += FoldNum )
		PostPreDec( &Zin[m], Eparity[i], &EparityBf[i][0], EfoldNum );

	for( i=0; i<TrellisLen; i++ ) Eparity[i] = Parity[i];

	// "ParityRequire" -- extrinstic information of parity are provided
	// "ParityNotRequire" -- extrinstic information of parity are NOT provided
	if( strcmp( Flag, "ParityRequire" ) == 0 )
	{
		if( FfoldNum == 0 ) for( i=0; i<TrellisLen; i++ ) Parity[i] = RClip3( Parity[i] * Fparity[i] );
		else if( FfoldNum == 1 ) for( m=EfoldNum, i=0; i<TrellisLen; i++, m+=FoldNum ) 
			Parity[i] = RClip3( Parity[i] * Plr( Fparity[i], Zin[m] ) );
		else for( m=EfoldNum, i=0; i<TrellisLen; i++, m+=FoldNum )
			Parity[i] = RClip3( Parity[i] * Plr( Fparity[i], Plr( Zin[m], FparityBf[i][0] ) ) );
	}

	if( FfoldNum != 0 ) for( m=EfoldNum, i=0; i<TrellisLen; i++, m+=FoldNum )
		PostPreDec( &Zin[m], Plr( Eparity[i], Fparity[i] ), &FparityBf[i][0], FfoldNum );

	for( i=0; i<DataLen; i++ )	ExInf[Nuser][dim][i] = RClip3( Zin[i] / tmpbf[i] );
}
/******************************************************************************************
function : sub-function of SpcCoDec::Aposteriori( int Nuser, int dim, float *Parity )
******************************************************************************************/
void SpcCoDec::PostPreDec( float *Zin, float Parity, float *bf, int blength )
{
	int i;
	for( i=0; i<blength; i++ ) Zin[i] = RClip3( Zin[i] * Plr( Parity, bf[i] ) );
}
/******************************************************************************************
function : prepare the information for the BCJR
******************************************************************************************/
void SpcCoDec::Apriori( float *Parity, float *Tail )
{
	int i, j, m;

	for( m=0, i=0; i<TrellisLen; i++ )
	{	
		PrePrtDec( &Zin[m], &Eparity[i], &EparityBf[i][0], EfoldNum );
		m += EfoldNum;
		PrePrtDec( &Zin[m], &Fparity[i], &FparityBf[i][0], FfoldNum );
		m += FfoldNum;
	}

	if( FfoldNum == 0 ) for( i=0; i<TrellisLen; i++ ) Fparity[i] = Parity[i];
	else for( i=0; i<TrellisLen; i++ ) Fparity[i] = Plr( Fparity[i], Parity[i] );

	for( i=0; i<TrellisLen; i++ )
	{
		for( j=0; j<NodeNum; j++ )
		{
			if( ConvSpc.Bcode[j][0] == 0 ) ConvSpc.CodeInf[i][j][0] = RClip3( Eparity[i] * Fparity[i] );
			else ConvSpc.CodeInf[i][j][0] = Eparity[i];
			if( ConvSpc.Bcode[j][1] == 0 ) ConvSpc.CodeInf[i][j][1] = Fparity[i];
			else ConvSpc.CodeInf[i][j][1] = 1;