turbo.c

zigzag纠错码是一种非常简单但却具有非常优秀性能的码。非常适合用于中码率或者高码率编码的系统。此为一位教授提供的源码。

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h>
#include "input.h"
#define TINY		1.0e-50
#define TOOSMALL	1.0e-7
#define TOOLARGE	1.0e7
#define STEP(x)		((x)>1.0? 0:1)

int	DataLength;
int	CodeLength;

int	IValue		= IVALUE;
int	JValue		= JVALUE; 
int	KValue		= KVALUE;
int	DataLength	= DATALENGTH;

int	ShaffleRule[KVALUE][DATALENGTH];
float	ExInf[KVALUE][DATALENGTH]; 
float	Step1Buffer[IVALUE][JVALUE]; 
float	Step1Out[IVALUE];  

extern int	IterationNum;
extern float	CodingRate;

void Prepare( );

void Channel( float *sb, double sigma );
void Coder( int *r, float *sb );
void Coding( int *v, int *d, int *mem );

void Decoder( int *g, float *sb, double sigma );
void Turbo( float *z, float *sb );
void TurboDecoding( float *sb, float *zin );
void Step1( float *zin, float *mapout, int blength );
void Step2( float *w, float *sb );
void Step3( float *w, float *zin );
void Clip( float *x, int length );
float Plr( float x, float y ); 

void SetShafflerule( );
void Shaffle( int *vin, int *vout, int *rule, int n );
void Scramble( int *rule, int length );

double GNoise( double s );

extern void Error( char *message );


// Initialisation of parameters and interleavers.
void Prepare( )
{
	if( DataLength == 0 ) Error( "Data length = 0" );
	if( JValue == 0 ) Error( "J value = 0" );
	if( DataLength % JValue != 0 ) Error( "Data length is not divisible by J value" );

	CodeLength = DataLength + IValue * KValue;
	CodingRate = (float) DataLength / (float) CodeLength;
	SetShafflerule( );
}

// Global multi-dimensional turbo-type decoder.
void Turbo( float *z, float *sb )
{
	int i, j;
	float zin[DATALENGTH], buffer[DATALENGTH];

	for( i=0; i<KValue; i++ )
	{
		for( j=0; j<DataLength; j++ )
			buffer[j] = zin[j] = z[ShaffleRule[i][j]] / ExInf[i][j];

		TurboDecoding( &sb[IValue*i], zin );

		for( j=0; j<DataLength; j++)
		{
			ExInf[i][j] = zin[j] / buffer[j];
			z[ShaffleRule[i][j]] = zin[j];
		}
	}
}

// AWGN channel.
void Channel( float *sb, double sigma )
{
	int i;
	for( i=0; i<CodeLength; i++ ) sb[i] += GNoise( sigma );
}


// Global encoder.
void Coder( int *r, float *sb )
{
	int i, j, n=0, mem;
	int OneDimv[IVALUE]; 
	int ShaffledData[DATALENGTH]; 

	for( i=0; i<KValue; i++ )
	{
		mem = 0;
		Shaffle( r, ShaffledData, &ShaffleRule[i][0], DataLength );
		Coding( OneDimv, ShaffledData, &mem );
		for( j=0; j<IValue; j++ )
			sb[n++] = (float) ( 1 - ( OneDimv[j] << 1 ) );
	}

	for( j=0; j<DataLength; j++ )
		sb[n++] = (float) ( 1 - ( r[j] << 1 ) );
}


//Local encoder.
void Coding( int *v, int *d, int *mem )
{
	int i, j, n=0;

	for( i=0; i<IValue; i++ )
	{
		for( j=0; j<JValue; j++ ) *mem = *mem ^ d[n++];
		v[i] = *mem;
	}
}

//Local decoder.
void TurboDecoding( float *sb, float *zin )
{
	float w[IVALUE];
	
	Step1( zin, Step1Out, JValue );
	Step2( w, sb );
	Step3( w, zin );
}


//This is the decoder for individul SPC code.
//The trellis part is not considered here.
//The outputs of step 1 are stored for step 3.
//Also the outputs for the parity bits will be used
//in Step 2 to decode the trellis part. 
void Step1( float *zin, float *mapout, int blength )
{
        int i, j, m = -1;
        float a;
        if( blength > 1 )
        {
                for( j=0; j<IValue; j++ )
                {
                        m += blength;
                        Step1Buffer[j][blength-2] = zin[m--];
                        for( i=blength-3; i>=0; i-- )
                                 Step1Buffer[j][i] = Plr( Step1Buffer[j][i+1], zin[m--] );

                        mapout[j] = Plr( Step1Buffer[j][0], zin[m] );

                        a = zin[m++];
                        for( i=1; i<blength-1; i++ )
                        {
                                Step1Buffer[j][i] = Plr( Step1Buffer[j][i], a );
                                a = Plr(zin[m++], a);
                        }
                        Step1Buffer[j][i] = a;
                }
        }
        else // If blength = 1, the processing is simpler.
        {
                for( j=0; j<IValue; j++ )
                {
                        Step1Buffer[j][0] = TOOLARGE;
                        mapout[j] = zin[j];
                }
        }
}


//Forward and backward search
void Step2( float *w, float *sb )
{
	int i;
	float a;

	//Backward search
        w[IValue-1] = sb[IValue-1];
        for( i=IValue-1; i>0; i-- )
                w[i-1] = sb[i-1] * Plr( w[i], Step1Out[i] );
	
	//Forward search
	a = sb[0] * Step1Out[0];
	for( i=1; i<IValue; i++ )
	{
		w[i] = Plr( w[i], a );
		a = sb[i] * Plr( a, Step1Out[i] );
	}
}


//This is the decoder for individul SPC code again. 
//The APP output from the trellis part is used here. 
void Step3( float *w, float *zin )
{
        int i, j, m;

        for( i=0; i<JValue; i++ )
                zin[i] *= Plr( w[0], Step1Buffer[0][i] );

        m = JValue;
        for( i=1; i<IValue; i++ )
                for( j=0; j<JValue; j++ )
                        zin[m++] *= Plr( w[i], Step1Buffer[i][j] );
	
	Clip( zin, DataLength );
}


void Clip( float *x, int length )
{
	int i;
	for( i=0; i<length; i++ )
	{
      		if( x[i] > TOOLARGE ) x[i] = TOOLARGE;
		else if( x[i] < TOOSMALL ) x[i] = TOOSMALL;
	}
}

void Decoder( int *g, float *sb, double sigma )
{
	int i, j, k = KValue * IValue;
	float s2 = sigma * sigma / 2.0, z[DATALENGTH];

	for( i=0; i<k; i++ )
		sb[i] = exp( sb[i] / s2 ); 
	Clip( sb, k ); 


	for( j=0; j<DataLength; j++ )
	{
		z[j] = exp( sb[k++] / s2 );
		for( i=0; i<KValue; i++ ) 
			ExInf[i][j] = 1.0;
	}
	Clip( z, DataLength ); 

	for( i=0; i<IterationNum; i++ ) Turbo( z, sb );

	for( i=0; i<DataLength; i++ ) g[i] = STEP( z[i] );
}

//AWGN noise generator.
double GNoise( double s )
{
	double x;
	x = drand48();
	if( x < TINY ) x = TINY;
	return (s * sqrt( -2.0 * log( x ) ) * cos( 6.283185 * drand48( ) ) );
}

//Interleaving
void Shaffle( int *vin, int *vout, int *rule, int n )
{
	int i;
	for( i=0; i<n; i++ ) vout[i] = vin[rule[i]];
}

//The basic parity APP decoding function. 
float Plr( float x, float y ) 
{
	return ( x * y + 1.0 ) / ( x + y ); 
}

//Generate K random interleavers.
void SetShafflerule( )
{
	int i, j;
	for( i=0; i<KValue; i++ )
	{
		for( j=0; j<DataLength; j++ ) ShaffleRule[i][j] = j;
		if( i > 0 ) Scramble( ShaffleRule[i], DataLength );
	}
}

//Random interleaver generator.
void Scramble( int *rule, int length )
{
	int i, j, k;
	for( i=0; i<length; i++ ) rule[i] = i;
	for( i=0; i<length-1; i++ )
	{
		k = i + lrand48( ) % ( length - i );
		j = rule[i];
		rule[i] = rule[k];
		rule[k] = j;
	}
}