rs.c

应用于无线通信的RS编码和OFDM调制源程序。64点FFT

#include <stdio.h>
#include <stdlib.h>
#include "Setting.H"
#include "RS.H"

/* Primitive polynomials */
#if(RSInforBit == 2)
WordType Pp[RSInforBit+1] = { 1, 1, 1 };

#elif(RSInforBit == 3)
/* 1 + x + x^3 */
WordType Pp[RSInforBit+1] = { 1, 1, 0, 1 };

#elif(RSInforBit == 4)
/* 1 + x + x^4 */
WordType Pp[RSInforBit+1] = { 1, 1, 0, 0, 1 };

#elif(RSInforBit == 5)
/* 1 + x^2 + x^5 */
WordType Pp[RSInforBit+1] = { 1, 0, 1, 0, 0, 1 };

#elif(RSInforBit == 6)
/* 1 + x + x^6 */
WordType Pp[RSInforBit+1] = { 1, 1, 0, 0, 0, 0, 1 };

#elif(RSInforBit == 7)
/* 1 + x^3 + x^7 */
WordType Pp[RSInforBit+1] = { 1, 0, 0, 1, 0, 0, 0, 1 };

#elif(RSInforBit == 8)
/* 1+x^2+x^3+x^4+x^8 */
WordType Pp[RSInforBit+1] = { 1, 0, 1, 1, 1, 0, 0, 0, 1 };

#elif(RSInforBit == 9)
/* 1+x^4+x^9 */
WordType Pp[RSInforBit+1] = { 1, 0, 0, 0, 1, 0, 0, 0, 0, 1 };

#elif(RSInforBit == 10)
/* 1+x^3+x^10 */
WordType Pp[RSInforBit+1] = { 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1 };

#elif(RSInforBit == 11)
/* 1+x^2+x^11 */
WordType Pp[RSInforBit+1] = { 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 };

#elif(RSInforBit == 12)
/* 1+x+x^4+x^6+x^12 */
WordType Pp[RSInforBit+1] = { 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1 };

#elif(RSInforBit == 13)
/* 1+x+x^3+x^4+x^13 */
WordType Pp[RSInforBit+1] = { 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 };

#elif(RSInforBit == 14)
/* 1+x+x^6+x^10+x^14 */
WordType Pp[RSInforBit+1] = { 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 };

#elif(RSInforBit == 15)
/* 1+x+x^15 */
WordType Pp[RSInforBit+1] = { 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 };

#elif(RSInforBit == 16)
/* 1+x+x^3+x^12+x^16 */
WordType Pp[RSInforBit+1] = { 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1 };

#else
#error "RSInforBit must be in range 2 to 16"
#endif

/* Index->polynomial form conversion table */
DType Alpha_to[RSTotalSym + 1];

/* Polynomial->index form conversion table */
DType Index_of[RSTotalSym + 1];

/* No legal value in index form represents zero, so
 * we need a special value for this purpose
 */
#define A0 (RSTotalSym)

/* Alpha exponent for the first root of the generator polynomial */
#define B0	1

/* Generator polynomial g(x)
 * Degree of g(x) = RSTotalSym - RSInforSym
 * has roots @^B0, @^(B0+1), ... ,@^(B0+RSTotalSym-RSInforSym-1)
 */
DType Gg[RSTotalSym - RSInforSym + 1];

/* Compute x (mod) RSTotalSym, where RSTotalSym is 2^RSInforBit - 1,
 * without a slow division
 */
static WordType ModRSTotalSym(WordType x)
{
	while (x >= RSTotalSym)
	{
		x -= RSTotalSym;
		x = (WordType) ((x >> RSInforBit) + (x & RSTotalSym));
	}
	return x;
}

#define	CLEAR(a, n) {\
	WordType ci;\
	for(ci = (n) - 1; ci >= 0; ci--)\
		(a)[ci] = 0;\
	}

#define	COPY(a, b, n) {\
	WordType ci;\
	for(ci = (n) - 1;ci >= 0; ci--)\
		(a)[ci] = (b)[ci];\
	}
#define	COPYDOWN(a, b, n) {\
	WordType ci;\
	for(ci = (n) - 1; ci >= 0; ci--)\
		(a)[ci] = (b)[ci];\
	}

void Init_RS(void)                          /* Initialization function       */
{
	Generate_GF();
	Generate_Poly();
}

void Generate_GF(void)                      /* Generate Galois Field         */
{
	register DType i, mask;

	mask = 1;
	Alpha_to[RSInforBit] = 0;
	for (i = 0; i < RSInforBit; i++)
	{
		Alpha_to[i] = mask;
		Index_of[Alpha_to[i]] = i;
		/* If Pp[i] == 1 then, term @^i occurs in poly-repr of @^RSInforBit  */
		if (Pp[i] != 0)
		{
			Alpha_to[RSInforBit] ^= mask;	      /* Bit-wise EXOR operation */
		}
		mask <<= 1;                               /* Single left-shift       */
	}
	Index_of[Alpha_to[RSInforBit]] = RSInforBit;
	/* Have obtained poly-repr of @^RSInforBit. Poly-repr of @^(i+1) is given by
	 * poly-repr of @^i shifted left one-bit and accounting for any @^RSInforBit
	 * term that may occur when poly-repr of @^i is shifted.
	 */
	mask >>= 1;
	for (i = RSInforBit + 1; i < RSTotalSym; i++)
	{
		if (Alpha_to[i - 1] >= mask)
		{
			Alpha_to[i] = Alpha_to[RSInforBit] ^ ((Alpha_to[i - 1] ^ mask) << 1);
		}
		else
		{
			Alpha_to[i] = Alpha_to[i - 1] << 1;
		}
		Index_of[Alpha_to[i]] = i;
	}
	Index_of[0] = A0;
	Alpha_to[RSTotalSym] = 0;
}

void Generate_Poly(void)                    /* Generate generator polynomial */
{
	register WordType i, j;

	Gg[0] = Alpha_to[B0];
	Gg[1] = 1;                              /* g(x) = (X+@^B0) initially    */
	for (i = 2; i < RSTotalSym - RSInforSym +1; i++)
	{
		Gg[i] = 1;
		/* Below multiply (Gg[0] + Gg[1]*x + ... + Gg[i]*x^i) by
		 * (@**(B0+i-1) + x)
		 */
		for (j = i - 1; j > 0; j--)
		{
			if (Gg[j] != 0)
			{
				Gg[j] = Gg[j - 1] ^ 
					Alpha_to[ModRSTotalSym((WordType) (Index_of[Gg[j]] + B0 + i - 1))];
			}
			else
			{
				Gg[j] = Gg[j - 1];
			}
		}
		/* Gg[0] can never be zero */
		Gg[0] = Alpha_to[ModRSTotalSym((WordType) (Index_of[Gg[0]] + B0 + i - 1))];
	}
	/* convert Gg[] to index form for quicker encoding */
	for (i = 0; i < RSTotalSym - RSInforSym + 1; i++)
	{
		Gg[i] = Index_of[Gg[i]];
	}
}

void Encode_RS(DType SourceData[RSInforSym], DType PariData[RSTotalSym - RSInforSym])
{
	register WordType i, j;
	DType FeedBack;
    
	CLEAR(PariData, RSTotalSym - RSInforSym);

	for (i = RSInforSym - 1; i >= 0; i--)
	{
		if(SourceData[i] > RSTotalSym)      /* Illegal symbol                */
		{
			printf("Illegal symbol occurs when encoding !\n");
			exit(0);
		}
		/* Change to Index_of form */
		FeedBack = Index_of[SourceData[i] ^ PariData[RSTotalSym - RSInforSym - 1]];
		if (FeedBack != A0)
		{
			/* FeedBack term is non-zero */
			for (j = RSTotalSym - RSInforSym - 1; j > 0; j--)
			{
				if (Gg[j] != A0)
				{
					PariData[j] = PariData[j - 1] ^ 
						Alpha_to[ModRSTotalSym((WordType) (Gg[j] + FeedBack))];
				}
				else
				{
					PariData[j] = PariData[j - 1];
				}
			}
			PariData[0] = Alpha_to[ModRSTotalSym((WordType) (Gg[0] + FeedBack))];
		}
		else
		{
		    /* FeedBack term is zero. Encoder becomes a single-byte shifter  */
			for (j = RSTotalSym - RSInforSym - 1; j > 0; j--)
			{
				PariData[j] = PariData[j - 1];
			}
			PariData[0] = 0;
		}
	}
}

WordType Decode_RS(DType ComData[RSTotalSym])
{
	WordType deg_lambda, el, deg_omega;
	WordType i, j, r;
	DType q, tmp, num1, num2, den, discr_r;
	DType recd[RSTotalSym];
	/* Error Locator poly and syndrome poly */
	DType lambda[RSTotalSym - RSInforSym + 1], s[RSTotalSym - RSInforSym + 1];
	DType b[RSTotalSym - RSInforSym + 1], t[RSTotalSym - RSInforSym + 1];
	DType omega[RSTotalSym - RSInforSym + 1], reg[RSTotalSym - RSInforSym + 1];
	WordType root[RSTotalSym - RSInforSym], loc[RSTotalSym - RSInforSym];
	DType syn_error;
	WordType count;

	/* ComData[] is in polynomial form, copy and convert to index form */
	for (i = RSTotalSym - 1; i >= 0; i--)
	{
		if(ComData[i] > RSTotalSym)         /* Illegal symbol                */
		{
			printf("Illegal symbol occurs when decoding !\n");
			exit(0);
		}
		recd[i] = Index_of[ComData[i]];
	}
	/* First form the syndromes; i.e., evaluate recd(x) at roots of g(x)
	 * namely @^(B0+i), i = 0, ... ,(RSTotalSym - RSInforSym - 1)
	 */
	syn_error = 0;
	for (i = 1; i <= RSTotalSym - RSInforSym; i++)
	{
		tmp = 0;
		for (j = 0; j < RSTotalSym; j++)
		{
			if (recd[j] != A0)              /* recd[j] in index form         */
				tmp ^= Alpha_to[ModRSTotalSym((WordType) (recd[j] + (B0 + i - 1) * j))];
		}
		syn_error |= tmp;                   /* Set flag if non-zero syndrome =>
		                                     * error                         */
		/* Store syndrome in index form */
		s[i] = Index_of[tmp];
	}
	if (!syn_error)
	{
		/* If syndrome is zero, ComData[] is a codeword and there is no
		 * error to correct. So return ComData[] unmodified
		 */
		return 0;
	}

	CLEAR(&lambda[1], RSTotalSym - RSInforSym);
	lambda[0] = 1;

	for(i=0; i < RSTotalSym - RSInforSym + 1; i++)
	{
		b[i] = Index_of[lambda[i]];
	}

	/* Begin Berlekamp-Massey algorithm to determine error+erasure
	 * locator polynomial
	 */
	r = 0;
	el = 0;
	while (++r <= RSTotalSym - RSInforSym)
	{
		/* r is the step number */
		/* Compute discrepancy at the r-th step in poly-form */
		discr_r = 0;
		for (i = 0; i < r; i++)
		{
			if ((lambda[i] != 0) && (s[r - i] != A0))
			{
				discr_r ^= Alpha_to[ModRSTotalSym((WordType) (Index_of[lambda[i]] + s[r - i]))];
			}
		}
		discr_r = Index_of[discr_r];        /* Stored in index form          */
		if (discr_r == A0)
		{
			/* 2 lines below: B(x) <-- x*B(x) */
			COPYDOWN(&b[1], b, RSTotalSym - RSInforSym);
			b[0] = A0;
		}
		else
		{
			/* 7 lines below: T(x) <-- lambda(x) - discr_r * x*b(x) */
			t[0] = lambda[0];
			for (i = 0 ; i < RSTotalSym - RSInforSym; i++)
			{
				if(b[i] != A0)
					t[i+1] = lambda[i+1] ^ 
					Alpha_to[ModRSTotalSym((WordType) (discr_r + b[i]))];
				else
					t[i+1] = lambda[i+1];
			}
			if (2 * el <= r - 1)
			{
				el = r - el;
				/* 2 lines below: B(x) <-- inv(discr_r) * lambda(x) */
				for (i = 0; i <= RSTotalSym - RSInforSym; i++)
					b[i] = (lambda[i] == 0) ? A0 : 
				ModRSTotalSym((WordType) (Index_of[lambda[i]] - discr_r + RSTotalSym));
			}
			else
			{
				/* 2 lines below: B(x) <-- x*B(x) */
				COPYDOWN(&b[1], b, RSTotalSym - RSInforSym);
				b[0] = A0;
			}
			COPY(lambda, t, RSTotalSym - RSInforSym + 1);
		}
	}

	/* Convert lambda to index form and compute deg(lambda(x)) */
	deg_lambda = 0;
	for(i = 0; i < RSTotalSym - RSInforSym + 1; i++)
	{
		lambda[i] = Index_of[lambda[i]];
		if(lambda[i] != A0)
			deg_lambda = i;
	}

	/* Find roots of the error locator polynomial by Chien search */
	COPY(&reg[1], &lambda[1], RSTotalSym - RSInforSym);
	count = 0;                              /* Number of roots of lambda(x)  */
	for (i = 1; i <= RSTotalSym; i++)
	{
		q = 1;
		for (j = deg_lambda; j > 0; j--)
			if (reg[j] != A0)
			{
				reg[j] = ModRSTotalSym((WordType) (reg[j] + j));
				q ^= Alpha_to[reg[j]];
			}
		if (!q)
		{
			/* Store root (index-form) and error location number */
			root[count] = i;
			loc[count] = RSTotalSym - i;
			count++;
		}
	}

	if (deg_lambda != count)
	{
		/*
		 * deg(lambda) unequal to number of roots => uncorrectable
		 * error detected
		 */
		return -1;
	}
	if (deg_lambda == 0)
	{
		/*
		 * deg(lambda) = 0, there's no root => uncorrectable
		 * error detected
		 */
		return -1;
	}

	/*
	 * Compute error evaluator poly omega(x) = s(x)*lambda(x) (modulo
	 * x^(RSTotalSym-RSInforSym)). in index form. Also find deg(omega).
	 */
	deg_omega = 0;
	for (i = 0; i < RSTotalSym - RSInforSym; i++)
	{
		tmp = 0;
		j = (deg_lambda < i) ? deg_lambda : i;
		for(; j >= 0; j--)
		{
			if ((s[i + 1 - j] != A0) && (lambda[j] != A0))
				tmp ^= Alpha_to[ModRSTotalSym((WordType) (s[i + 1 - j] + lambda[j]))];
		}
		if(tmp != 0)
			deg_omega = i;
		omega[i] = Index_of[tmp];
	}
	omega[RSTotalSym-RSInforSym] = A0;

	/* Compute error values in poly-form.
	 * num1 = omega(inv(X(l))), num2 = inv(X(l))^(B0-1)
	 * and den = lambda_pr(inv(X(l))) all in poly-form
	 */
	for (j = count - 1; j >= 0; j--)
	{
		num1 = 0;
		for (i = deg_omega; i >= 0; i--)
		{
			if (omega[i] != A0)
				num1  ^= Alpha_to[ModRSTotalSym((WordType) (omega[i] + i * root[j]))];
		}
		num2 = Alpha_to[ModRSTotalSym((WordType) (root[j] * (B0 - 1) + RSTotalSym))];
		den = 0;

		/* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
		for (i = min(deg_lambda, RSTotalSym - RSInforSym - 1) & ~1; i >= 0; i -= 2)
		{
			if(lambda[i+1] != A0)
				den ^= Alpha_to[ModRSTotalSym((WordType) (lambda[i+1] + i * root[j]))];
		}
		if (den == 0)
		{
			return -1;
		}
		/* Apply error to ComData */
		if (num1 != 0)
		{
			ComData[loc[j]] ^= Alpha_to[ModRSTotalSym(
				(WordType) (Index_of[num1] + Index_of[num2] + RSTotalSym - Index_of[den]))];
		}
	}

	return count;
}