rs.cpp

QAM module to use in Java with an easy interface and powerful performance

//---------------------------------------------------------------------------

#pragma hdrstop

#include "RS.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
using namespace std;

__fastcall RS::RS()
{
        NPAR=4;
        MAXDEG=2*NPAR;
        DEBUG = FALSE;
        pBytes=new int[MAXDEG];
        synBytes=new int[MAXDEG];
        genPoly=new int[2*MAXDEG];
        Lambda=new int[MAXDEG];
        Omega=new int[MAXDEG];
        initialize_ecc(NPAR);
}

void __fastcall RS::initialize_ecc(int NumberOfParityBytes)
{

        if(NumberOfParityBytes!=NPAR)
        {
                delete [] pBytes;
                delete [] synBytes;
                delete [] genPoly;
                delete [] Lambda;
                delete [] Omega;
                NPAR=NumberOfParityBytes;
                MAXDEG=2*NPAR;
                pBytes=new int[MAXDEG];
                synBytes=new int[MAXDEG];
                genPoly=new int[2*MAXDEG];
                Lambda=new int[MAXDEG];
                Omega=new int[MAXDEG];
        }

        /* Initialize the galois field arithmetic tables */
        init_galois_tables();

        /* Compute the encoder generator polynomial */
        compute_genpoly(NPAR, genPoly);
}

void __fastcall RS::init_galois_tables()
{
        /* initialize the table of powers of alpha */
        init_exp_table();
}

void __fastcall RS::init_exp_table()
{
        int i, z;
        int pinit,p1,p2,p3,p4,p5,p6,p7,p8;

        pinit = p2 = p3 = p4 = p5 = p6 = p7 = p8 = 0;
        p1 = 1;

        gexp[0] = 1;
        gexp[255] = gexp[0];
        glog[0] = 0;			/* shouldn't log[0] be an error? */

        for (i = 1; i < 256; i++)
        {
                pinit = p8;
                p8 = p7;
                p7 = p6;
                p6 = p5;
                p5 = p4 ^ pinit;
                p4 = p3 ^ pinit;
                p3 = p2 ^ pinit;
                p2 = p1;
                p1 = pinit;
                gexp[i] = p1 + p2*2 + p3*4 + p4*8 + p5*16 + p6*32 + p7*64 + p8*128;
                gexp[i+255] = gexp[i];
        }
	
        for (i = 1; i < 256; i++)
        {
                for (z = 0; z < 256; z++)
                {
                        if (gexp[z] == i)
                        {
	                        glog[i] = z;
	                        break;
                        }
                }
        }
}

/* multiplication using logarithms */
int __fastcall RS::gmult(int a, int b)
{
        int i,j;
        if (a==0 || b == 0) return (0);
        i = glog[a];
        j = glog[b];
        return (gexp[i+j]);
}

int __fastcall RS::ginv(int elt)
{
        return (gexp[255-glog[elt]]);
}

/* Create a generator polynomial for an n byte RS code.
 * The coefficients are returned in the genPoly arg.
 * Make sure that the genPoly array which is passed in is
 * at least n+1 bytes long.
 */
void __fastcall RS::compute_genpoly(int nbytes, int genpoly[])
{
        int i, tp[256], tp1[256];
        /* multiply (x + a^n) for n = 1 to nbytes */
        zero_poly(tp1);
        tp1[0] = 1;
        for (i = 1; i <= nbytes; i++)
        {
                zero_poly(tp);
                tp[0] = gexp[i];		/* set up x+a^n */
                tp[1] = 1;

                mult_polys(genpoly, tp, tp1);
                copy_poly(tp1, genpoly);
        }
}

/* polynomial multiplication */
void __fastcall RS::mult_polys(int dst[], int p1[], int p2[])
{
        int i, j;
        int *tmp1=new int[MAXDEG*2];
        for (i=0; i < (MAXDEG*2); i++) dst[i] = 0;
        for (i = 0; i < MAXDEG; i++)
        {
                for(j=MAXDEG; j<(MAXDEG*2); j++) tmp1[j]=0;
                // scale tmp1 by p1[i]
                for(j=0; j<MAXDEG; j++) tmp1[j]=gmult(p2[j], p1[i]);
                // and mult (shift) tmp1 right by i
                for (j = (MAXDEG*2)-1; j >= i; j--) tmp1[j] = tmp1[j-i];
                for (j = 0; j < i; j++) tmp1[j] = 0;
                // add into partial product
                for(j=0; j < (MAXDEG*2); j++) dst[j] ^= tmp1[j];
        }
        delete [] tmp1;
}

void __fastcall RS::zero_poly(int poly[])
{
        int i;
        for (i = 0; i < MAXDEG; i++) poly[i] = 0;
}

void __fastcall RS::copy_poly(int dst[], int src[])
{
        int i;
        for (i = 0; i < MAXDEG; i++) dst[i] = src[i];
}


/* Simulate a LFSR with generator polynomial for n byte RS code. 
 * Pass in a pointer to the data array, and amount of data. 
 *
 * The parity bytes are deposited into pBytes[], and the whole message
 * and parity are copied to dest to make a codeword.
 * 
 */
void __fastcall RS::encode_data(unsigned char msg[], int nbytes, unsigned char dst[])
{
        int i,dbyte, j;
        int *LFSR=new int[NPAR+1];
        for(i=0; i < NPAR+1; i++) LFSR[i]=0;
        for (i = 0; i < nbytes; i++)
        {
                dbyte = msg[i] ^ LFSR[NPAR-1];
                for (j = NPAR-1; j > 0; j--)
                {
                        LFSR[j] = LFSR[j-1] ^ gmult(genPoly[j], dbyte);
                }
                LFSR[0] = gmult(genPoly[0], dbyte);
        }

        for (i = 0; i < NPAR; i++) pBytes[i] = LFSR[i];
        build_codeword(msg, nbytes, dst);
        delete [] LFSR;
}

/* Append the parity bytes onto the end of the message */
void __fastcall RS::build_codeword(unsigned char msg[], int nbytes, unsigned char dst[])
{
        int i;
        for (i = 0; i < nbytes; i++) dst[i] = msg[i];
        for (i = 0; i < NPAR; i++) dst[i+nbytes] = pBytes[NPAR-1-i];
}

void __fastcall RS::encode_data(int msg[], int nbytes)
{
        int i,dbyte, j;
        int *LFSR=new int[NPAR+1];
        for(i=0; i < NPAR+1; i++) LFSR[i]=0;
        for (i = 0; i < nbytes; i++)
        {
                dbyte = msg[i] ^ LFSR[NPAR-1];
                for (j = NPAR-1; j > 0; j--)
                {
                        LFSR[j] = LFSR[j-1] ^ gmult(genPoly[j], dbyte);
                }
                LFSR[0] = gmult(genPoly[0], dbyte);
        }

        for (i = 0; i < NPAR; i++)
        {
                pBytes[i] = LFSR[i];
                msg[i+nbytes] = LFSR[NPAR-1-i];
        }
        delete [] LFSR;
}

void __fastcall RS::build_codeword(int msg[], int nbytes, int dst[])
{
        int i;
        for (i = 0; i < nbytes; i++) dst[i] = msg[i];
        for (i = 0; i < NPAR; i++) dst[i+nbytes] = pBytes[NPAR-1-i];
}


/**********************************************************
 * Reed Solomon Decoder 
 *
 * Computes the syndrome of a codeword. Puts the results
 * into the synBytes[] array.
 */
void __fastcall RS::decode_data(unsigned char data[], int nbytes)
{
        int i, j, sum;
        for (j = 0; j < NPAR;  j++)
        {
                sum = 0;
                for (i = 0; i < nbytes; i++)
                {
                        sum = data[i] ^ gmult(gexp[j+1], sum);
                }
                synBytes[j]  = sum;
        }
        for(j=NPAR;j<MAXDEG;j++)synBytes[j]=0;
}

void __fastcall RS::decode_data(int data[], int nbytes)
{
        int i, j, sum, tmp;
        for (j = 0; j < NPAR;  j++)
        {
                sum = 0;
                for (i = 0; i < nbytes; i++)
                {
                        tmp = data[i];
                        if(tmp<0)tmp=0;
                        sum = (unsigned char)(tmp ^ gmult(gexp[j+1], sum));
                        //sum = (unsigned char)tmp ^ gmult(gexp[j+1], sum);
                }
                synBytes[j]  = sum;
        }
        for(j=NPAR;j<MAXDEG;j++)synBytes[j]=0;
}

int __fastcall RS::check_syndrome()
{
        int i, nz = 0;
        for (i =0 ; i < NPAR; i++)
        {
                if (synBytes[i] != 0) nz = 1;
        }
        return nz;
}

/* Combined Erasure And Error Magnitude Computation
 *
 * Pass in the codeword, its size in bytes, as well as
 * an array of any known erasure locations, along the number
 * of these erasures.
 *
 * Evaluate Omega(actually Psi)/Lambda' at the roots
 * alpha^(-i) for error locs i.
 *
 * Returns 1 if everything ok, or 0 if an out-of-bounds error is found
 *
 */
int __fastcall RS::correct_errors_erasures(unsigned char codeword[], int csize, int nerasures, int erasures[])
{
        int r, i, j, err;
        /* If you want to take advantage of erasure correction, be sure to
           set NErasures and ErasureLocs[] with the locations of erasures.
        */
        NErasures = nerasures;
        for (i = 0; i < NErasures; i++) ErasureLocs[i] = erasures[i];

        Modified_Berlekamp_Massey();
        Find_Roots();


        if ((NErrors <= NPAR) && NErrors > 0)
        {
                UncorrectableCodeword=false;
                /* first check for illegal error locs */
                for (r = 0; r < NErrors; r++)
                {
                        if (ErrorLocs[r] >= csize)
                        {
	                        //if (DEBUG) fprintf(stderr, "Error loc i=%d outside of codeword length %d\n", i, csize);
	                        return(0);
                        }
                }

                for (r = 0; r < NErrors; r++)
                {
                        int num, denom;
                        i = ErrorLocs[r];
                        /* evaluate Omega at alpha^(-i) */

                        num = 0;
                        for (j = 0; j < MAXDEG; j++) num ^= gmult(Omega[j], gexp[((255-i)*j)%255]);

                        /* evaluate Lambda' (derivative) at alpha^(-i) ; all odd powers disappear */
                        denom = 0;
                        for (j = 1; j < MAXDEG; j += 2)
                        {
	                        denom ^= gmult(Lambda[j], gexp[((255-i)*(j-1)) % 255]);
                        }

                        err = gmult(num, ginv(denom));
                        //if (DEBUG) fprintf(stderr, "Error magnitude %#x at loc %d\n", err, csize-i);

                        codeword[csize-i-1] ^= err;
                }
                return(1);
        }
         else
         {