Reed Solomon Codec C code implement!
//this code is from Internet, not mine. thank for the author
//C code begin
/
/*
* Reed-Solomon coding and decoding
* This code supports a symbol size from 2 bits to 16 bits,
* implying a block size of 32-bit symbols (6 bits) up to 65535
* 16-bit symbols (1,048,560 bits).
* Note that if symbols larger than 8 bits are used, the type of each
* data array element switches from unsigned char to unsigned int. The
* caller must ensure that elements larger than the symbol range are
* not passed to the encoder or decoder.
*/
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <time.h>
#define MM 8 /* RS code over GF(2**MM) - change to suit */
#define KK 239 /* KK = number of information symbols */
//#define MM 4
//#define KK 9
#define NN ((1 << MM) - 1)
#define TT 8
//#define TT 3
#define DEBUG 1
#define PS 1
#if (MM < 8)
typedef unsigned char dtype;
#else
typedef unsigned int dtype;
#endif
#if (KK >= NN)
#error "KK must be less than 2**MM - 1"
#endif
/*primitive polynomials */
#if(MM == 2)
int pp[MM+1] = { 1, 1, 1 };
#elif(MM == 3)
/* 1 + x + x^3 */
int pp[MM+1] = { 1, 1, 0, 1 };
#elif(MM == 4)
/* 1 + x + x^4 */
int pp[MM+1] = { 1, 1, 0, 0, 1 };
#elif(MM == 5)
/* 1 + x^2 + x^5 */
int pp[MM+1] = { 1, 0, 1, 0, 0, 1 };
#elif(MM == 6)
/* 1 + x + x^6 */
int pp[MM+1] = { 1, 1, 0, 0, 0, 0, 1 };
#elif(MM == 7)
/* 1 + x^3 + x^7 */
int pp[MM+1] = { 1, 0, 0, 1, 0, 0, 0, 1 };
#elif(MM == 8)
/* 1+x^2+x^3+x^4+x^8 */
int pp[MM+1] = { 1, 0, 1, 1, 1, 0, 0, 0, 1 };
#elif(MM == 9)
/* 1+x^4+x^9 */
int pp[MM+1] = { 1, 0, 0, 0, 1, 0, 0, 0, 0, 1 };
#elif(MM == 10)
/* 1+x^3+x^10 */
int pp[MM+1] = { 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1 };
#elif(MM == 11)
/* 1+x^2+x^11 */
int pp[MM+1] = { 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 };
#elif(MM == 12)
/* 1+x+x^4+x^6+x^12 */
int pp[MM+1] = { 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1 };
#elif(MM == 13)
/* 1+x+x^3+x^4+x^13 */
int pp[MM+1] = { 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 };
#elif(MM == 14)
/* 1+x+x^6+x^10+x^14 */
int pp[MM+1] = { 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 };
#elif(MM == 15)
/* 1+x+x^15 */
int pp[MM+1] = { 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 };
#elif(MM == 16)
/* 1+x+x^3+x^12+x^16 */
int pp[MM+1] = { 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1 };
#else
#error "MM must be in range 2-16"
#endif
int alpha_to[NN+1];
int index_of[NN+1];
int gg[NN-KK+1];
int recd[NN];
int err_flag;
dtype data[KK],bb[NN-KK];
// int elp[NN-KK+2][NN-KK], d[NN-KK+2], l[NN-KK+2], u_lu[NN-KK+2], s[NN-KK+1] ;
// int count=0, syn_error=0, root[TT], loc[TT], z[TT+1], err[NN], reg[TT+1] ;
// int i,j,u,q ;
void generate_gf()
/* generate GF(2**MM) from the irreducible polynomial p(X) in pp[0]..pp[MM]
lookup tables: index->polynomial form alpha_to[] contains j=alpha**i;
polynomial form -> index form index_of[j=alpha**i] = i
alpha=2 is the primitive element of GF(2**MM)
*/
{
register int i, mask ;
mask = 1 ;
alpha_to[MM] = 0 ;
for (i=0; i<MM; i++)
{ alpha_to[i] = mask ;
index_of[alpha_to[i]] = i ;
if (pp[i]!=0)
alpha_to[MM] ^= mask ;
mask <<= 1 ;
}
index_of[alpha_to[MM]] = MM ;
mask >>= 1 ;
for (i=MM+1; i<NN; i++)
{ if (alpha_to[i-1] >= mask)
alpha_to[i] = alpha_to[MM] ^ ((alpha_to[i-1]^mask)<<1) ;
else alpha_to[i] = alpha_to[i-1]<<1 ;
index_of[alpha_to[i]] = i ;
}
index_of[0] = -1 ;
}
void gen_poly()
/* Obtain the generator polynomial of the TT-error correcting, length
NN=(2**MM -1) Reed Solomon code from the product of (X+alpha**i), i=1..2*TT
*/
{
register int i,j ;
gg[0] = 2 ; /* primitive element alpha = 2 for GF(2**MM) */
gg[1] = 1 ; /* g(x) = (X+alpha) initially */
for (i=2; i<=NN-KK; i++)
{ gg[i] = 1 ;
for (j=i-1; j>0; j--)
if (gg[j] != 0) gg[j] = gg[j-1]^ alpha_to[(index_of[gg[j]]+i)%NN] ;
else gg[j] = gg[j-1] ;
gg[0] = alpha_to[(index_of[gg[0]]+i)%NN] ; /* gg[0] can never be zero */
}
/* convert gg[] to index form for quicker encoding */
for (i=0; i<=NN-KK; i++) gg[i] = index_of[gg[i]] ;
}
void encode_rs()
/* take the string of symbols in data[i], i=0..(k-1) and encode systematically
to produce 2*TT parity symbols in bb[0]..bb[2*TT-1]
data[] is input and bb[] is output in polynomial form.
Encoding is done by using a feedback shift register with appropriate
coNNections specified by the elements of gg[], which was generated above.
Codeword is c(X) = data(X)*X**(NN-KK)+ b(X) */
{
register int i,j ;
int feedback ;
for (i=0; i<NN-KK; i++) bb[i] = 0 ;
for (i=KK-1; i>=0; i--)
{ feedback = index_of[data[i]^bb[NN-KK-1]] ;
if (feedback != -1)
{ for (j=NN-KK-1; j>0; j--)
if (gg[j] != -1)
bb[j] = bb[j-1]^alpha_to[(gg[j]+feedback)%NN] ;
else
bb[j] = bb[j-1] ;
bb[0] = alpha_to[(gg[0]+feedback)%NN] ;
}
else
{ for (j=NN-KK-1; j>0; j--)
bb[j] = bb[j-1] ;
bb[0] = 0 ;
} ;
} ;
} ;
void decode_rs()
/* assume we have received bits grouped into MM-bit symbols in recd[i],
i=0..(NN-1), and recd[i] is index form (ie as powers of alpha).
We first compute the 2*TT syndromes by substituting alpha**i into rec(X) and
evaluating, storing the syndromes in s[i], i=1..2TT (leave s[0] zero) .
Then we use the Berlekamp iteration to find the error location polynomial
elp[i]. If the degree of the elp is >TT, we cannot correct all the errors
and hence just put out the information symbols uncorrected. If the degree of
elp is <=TT, we substitute alpha**i , i=1..n into the elp to get the roots,
hence the inverse roots, the error location numbers. If the number of errors
located does not equal the degree of the elp, we have more than TT errors
and cannot correct them. Otherwise, we then solve for the error value at
the error location and correct the error. The procedure is that found in
Lin and Costello. For the cases where the number of errors is known to be too
large to correct, the information symbols as received are output (the
advantage of systematic encoding is that hopefully some of the information
symbols will be okay and that if we are in luck, the errors are in the
parity part of the transmiTTed codeword). Of course, these insoluble cases
can be returned as error flags to the calling routine if desired. */
{
register int i,j,u,q ;
int elp[NN-KK+2][NN-KK], d[NN-KK+2], l[NN-KK+2], u_lu[NN-KK+2], s[NN-KK+1] ;
int count=0, syn_error=0, root[TT], loc[TT], z[TT+1], err[NN], reg[TT+1] ;
#ifdef PS
int syndrome[NN-KK+1];
int omega[NN-KK+2][NN-KK];
int zx[TT+1];
#endif
err_flag=0;
/* first form the syndromes */
for (i=1; i<=NN-KK; i++)
{ s[i] = 0 ;
for (j=0; j<NN; j++)
if (recd[j]!=-1)
s[i] ^= alpha_to[(recd[j]+i*j)%NN] ; /* recd[j] in index form */
/* convert syndrome from polynomial form to index form */
if (s[i]!=0) syn_error=1 ; /* set flag if non-zero syndrome => error */
#ifdef PS
syndrome[i]=s[i];
#endif
s[i] = index_of[s[i]] ;
} ;
if (syn_error) /* if errors, try and correct */
{
/* compute the error location polynomial via the Berlekamp iterative algorithm,
following the terminology of Lin and Costello : d[u] is the 'mu'th
discrepancy, where u='mu'+1 and 'mu' (the Greek leTTer!) is the step number
ranging from -1 to 2*TT (see L&C), l[u] is the
degree of the elp at that step, and u_l[u] is the difference between the
step number and the degree of the elp.
*/
/* initialise table entries */
err_flag=1;
d[0] = 0 ; /* index form */
d[1] = s[1] ; /* index form */
elp[0][0] = 0 ; /* index form */
elp[1][0] = 1 ; /* polynomial form */
for (i=1; i<NN-KK; i++)
{ elp[0][i] = -1 ; /* index form */
elp[1][i] = 0 ; /* polynomial form */
}
l[0] = 0 ;
l[1] = 0 ;
u_lu[0] = -1 ;
u_lu[1] = 0 ;
u = 0 ;
do
{
u++ ;
if (d[u]==-1)
{ l[u+1] = l[u] ;
for (i=0; i<=l[u]; i++)
{ elp[u+1][i] = elp[u][i] ;
#ifdef PS
omega[u][i]=elp[u][i];
#endif
elp[u][i] = index_of[elp[u][i]] ;
}
}
else
/* search for words with greatest u_lu[q] for which d[q]!=0 */
{ q = u-1 ;
while ((d[q]==-1) && (q>0)) q-- ;
/* have found first non-zero d[q] */
if (q>0)
{ j=q ;
do
{ j-- ;
if ((d[j]!=-1) && (u_lu[q]<u_lu[j]))
q = j ;
}while (j>0) ;
} ;
/* have now found q such that d[u]!=0 and u_lu[q] is maximum */
/* store degree of new elp polynomial */
if (l[u]>l[q]+u-q) l[u+1] = l[u] ;
else l[u+1] = l[q]+u-q ;
/* form new elp(x) */
// for (i=0; i<NN-KK; i++) elp[u+1][i] = 0 ;
for (i=0; i<9; i++) elp[u+1][i] = 0 ;
for (i=0; i<=l[q]; i++)
if (elp[q][i]!=-1)
elp[u+1][i+u-q] = alpha_to[(d[u]+NN-d[q]+elp[q][i])%NN] ;
for (i=0; i<=l[u]; i++)
{ elp[u+1][i] ^= elp[u][i] ;
#ifdef PS
omega[u][i]=elp[u][i];
#endif
elp[u][i] = index_of[elp[u][i]] ; /*convert old elp value to index*/
}
}
u_lu[u+1] = u-l[u+1] ;
/* form (u+1)th discrepancy */
if (u<NN-KK) /* no discrepancy computed on last iteration */
{
if (s[u+1]!=-1)
d[u+1] = alpha_to[s[u+1]] ;
else
d[u+1] = 0 ;
for (i=1; i<=l[u+1]; i++)
if ((s[u+1-i]!=-1) && (elp[u+1][i]!=0))
d[u+1] ^= alpha_to[(s[u+1-i]+index_of[elp[u+1][i]])%NN] ;
d[u+1] = index_of[d[u+1]] ; /* put d[u+1] into index form */
}
} while ((u<NN-KK) && (l[u+1]<=TT)) ;
//------------------------------------------------------------------------------------------------------------------
u++ ;
if (l[u]<=TT) /* can correct error *///*****************************************
{
/* put elp into index form */
for (i=0; i<=l[u]; i++) elp[u][i] = index_of[elp[u][i]] ;
/* find roots of the error location polynomial */
for (i=1; i<=l[u]; i++)
reg[i] = elp[u][i] ;
count = 0 ;
for (i=1; i<=NN; i++)
{ q = 1 ;
for (j=1; j<=l[u]; j++)
if (reg[j]!=-1)
{ reg[j] = (reg[j]+j)%NN ;
q ^= alpha_to[reg[j]] ;
} ;
if (!q) /* store root and error location number indices */
{ root[count] = i;
loc[count] = NN-i ;
count++ ;
};
} ;
if (count==l[u]) /* no. roots = degree of elp hence <= TT errors */
{
/* form polynomial z(x) */
for (i=1; i<=l[u]; i++) /* Z[0] = 1 always - do not need */
{ if ((s[i]!=-1) && (elp[u][i]!=-1))
z[i] = alpha_to[s[i]] ^ alpha_to[elp[u][i]] ;
else if ((s[i]!=-1) && (elp[u][i]==-1))
z[i] = alpha_to[s[i]] ;
else if ((s[i]==-1) && (elp[u][i]!=-1))
z[i] = alpha_to[elp[u][i]] ;
else
z[i] = 0 ;
for (j=1; j<i; j++)
if ((s[j]!=-1) && (elp[u][i-j]!=-1))
z[i] ^= alpha_to[(elp[u][i-j] + s[j])%NN] ;
#ifdef PS
zx[i]=z[i];
#endif
z[i] = index_of[z[i]] ; /* put into index form */
} ;
/* evaluate errors at locations given by error location numbers loc[i] */
for (i=0; i<NN; i++)
{ err[i] = 0 ;
if (recd[i]!=-1) /* convert recd[] to polynomial form */
recd[i] = alpha_to[recd[i]] ;
else recd[i] = 0 ;
}
for (i=0; i<l[u]; i++) /* compute numerator of error term first */
{ err[loc[i]] = 1; /* accounts for z[0] */
for (j=1; j<=l[u]; j++)
if (z[j]!=-1)
err[loc[i]] ^= alpha_to[(z[j]+j*root[i])%NN] ;
if (err[loc[i]]!=0)
{ err[loc[i]] = index_of[err[loc[i]]] ;
q = 0 ; /* form denominator of error term */
for (j=0; j<l[u]; j++)
if (j!=i)
q += index_of[1^alpha_to[(loc[j]+root[i])%255]] ;
q = q % NN ;
err[loc[i]] = alpha_to[(err[loc[i]]-q+NN)%NN] ;
recd[loc[i]] ^= err[loc[i]] ; /*recd[i] must be in polynomial form */
}
}
}
else /* no. roots != degree of elp => >TT errors and cannot solve */
{ err_flag=2;
for (i=0; i<NN; i++) /* could return error flag if desired */
{if (recd[i]!=-1) /* convert recd[] to polynomial form */
recd[i] = alpha_to[recd[i]] ;
else recd[i] = 0 ; /* just output received codeword as is */
}
}
}
else /* elp has degree has degree >TT hence cannot solve *///********************
{
err_flag=2;
for (i=0; i<NN; i++) /* could return error flag if desired */
{
if (recd[i]!=-1) /* convert recd[] to polynomial form */
recd[i] = alpha_to[recd[i]] ;
else recd[i] = 0 ; /* just output received codeword as is */
}
}
}
else /* no non-zero syndromes => no errors: output received codeword */
for (i=0; i<NN; i++)
if (recd[i]!=-1) /* convert recd[] to polynomial form */
recd[i] = alpha_to[recd[i]] ;
else recd[i] = 0 ;
}
void main(void)
{
register int i;
int acterr_num,err_pos[NN];
// int goon=1;
/* generate the Galois Field GF(2**MM) */
generate_gf() ;
#ifdef DEBUG
printf("Look-up tables for GF(2**%2d)\n",MM);
printf("i alpha_to[i] index_of[i]\n");
for (i=0; i<=NN; i++)
printf("%d\t%d\t\t%d\n",i,alpha_to[i],index_of[i]);
printf("\n\n") ;
// for(i=0; i<9;i++)
// printf("%x",pp[i]);
#endif
/*compute the generator polynomial*/
gen_poly();
#ifdef DEBUG
printf("the generator polynomial g(x)(index form)=\n");
printf(" i gg[i]\n");
for(i=NN-KK;i>=0;i--)
printf(" %d %d\n",i,gg[i]);
#endif
do{
printf("Reed-Solomon code is (%d,%d) over F(%d).\n",NN,KK,(1<<MM));
printf("This pragramm can correct %d errors, and delect %d errors.\n",TT,2*TT);
printf("\n");
printf("Please determine the number of the actual errors: ");
scanf("%d",&acterr_num);
printf("Please determine the position of the actual errors\n");
printf("(total %d ,range from 0 to %d): ",acterr_num,NN-1);
for(i=0;i<acterr_num;i++)
scanf("%d",&err_pos[i]);
printf("\n");
/*generate the random data*/
/* srand(unsigned(time(NULL)));
for(i=0;i<KK;i++)
data[i]=rand()&NN;
*/
//input the data manully
#ifdef PS
for (i=0;i<KK;i++)
data[i]=KK-i;
#else
for (i=0;i<KK;i++)
data[i]=KK-i;
#endif
printf("The to be transmitted data is:\n");
for(i=0;i<KK;i++)
printf("%d ",data[i]);
printf("\n\n");
encode_rs();
/*form the transmitted code,with parity plus data.*/
for(i=0;i<NN-KK;i++)
recd[i]=bb[i];
for(i=NN-KK;i<NN;i++)
recd[i]=data[i-NN+KK];
printf("Encoded data is:\n");
for(i=0;i<NN;i++)
printf("%d ",recd[i]);
printf("\n\n");
printf("The recieved code with errors is:\n");
for(i=0;i<acterr_num;i++)
//--------------------------------------------------------
#ifdef PS
recd[err_pos[i]]=(recd[err_pos[i]]^16)%(NN+1);
//for (i=0;i<8;i++)
//only one error
for (i=0;i<7;i++)
recd[254-i]=16;
// for (i=0;i<247;i++)
// recd[i] = 0;
/*
recd[254]=55;
recd[253]=55;
recd[252]=55;
recd[251]=55;
recd[246]=111;
recd[245]=111;
recd[244]=111;
recd[243]=111;
*/
//recd[240]=66;
#else
recd[err_pos[i]]=(recd[err_pos[i]]^rand())%(NN+1);
#endif
//----------------------------------------------------------
for(i=0;i<NN;i++)
printf("%d ",recd[i]);
printf("\n\n");
//****************************************************
for(i=0;i<NN;i++)
recd[i]=index_of[recd[i]];
decode_rs();
switch(err_flag)
{
case(0):
printf("This code has been correctly transmitted.\n\n");
break;
case(1):
printf("This code can be corrected.\n\n");
break;
case(2):
printf("This code cannot be corrected,\n");
printf("while the number of errors is greater than TT\n\n");
break;
}
printf("The result of the decoding is :\n");
for(i=0;i<NN;i++)
printf("%d ",recd[i]);
printf("\n\n");
printf("The correct information is maybe:\n");
for(i=NN-KK;i<NN;i++)
printf("%d ",recd[i]);
printf("\n\n");
printf("Press ESC to EXIT, press any other key to CONTINUE...\n");
}while(_getch()!=27);
}
///
//C code end
///
if you want to get the detailed file, contact with me via EastDragonInc@hotmail.com
RS16
最新推荐文章于 2023-02-07 18:00:28 发布