www.pudn.com > ReedSolomon_codec.rar > ReedSolomon.h
// ReedSolomon.h: interface for the CReedSolomon class.
//
//////////////////////////////////////////////////////////////////////
#if !defined(AFX_REEDSOLOMON_H__7A4D7FA4_C47A_4880_9E03_833194109EA0__INCLUDED_)
#define AFX_REEDSOLOMON_H__7A4D7FA4_C47A_4880_9E03_833194109EA0__INCLUDED_
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
#include
#include
#include
#include
#include
using namespace std;
#define DEBUG 1
/************************************************************************
类CReedSolomon用于RS码的解码及编码
有两种方式初始化: 使用下面的_INIT_DATA或者使用预定义的值(只指定symsize).
symsize表示2的多少次方, 码长为2**symsize.
genpoly生成多项式, 应该为不可约多项式, 可参照预定义值
fcs 应该小于(1 << symsize).
prim 不能为0, 且小于(1 << symsize)
nroots 应该小于(1 << symsize), 表示校验码的长度.
纠错能力为0~nroots/2.
************************************************************************/
struct _INIT_DATA
{
int symsize; // RS code over GF(2**symsize) = mm
int genpoly; // Generator polynomial
int fcs; // First consecutive root, index form
int prim;// Primitive element, index form
int nroots;// Number of generator roots = number of parity symbols
int ntrials; // no use
};
// 预定义的初始值数组
extern _INIT_DATA Tab[];
template
class CReedSolomon
{
public:
CReedSolomon(unsigned int symsize, unsigned int gfpoly, unsigned int m_fcr,
unsigned int m_prim, unsigned int m_nroots);
CReedSolomon(_INIT_DATA data);
CReedSolomon(unsigned int m);
void Encode(T* data, T* parity);
int Decode(T* data, T* eras_pos, int no_eras);
virtual ~CReedSolomon();
private:
inline int modnn(int x)
{
while (x >= m_nn)
{
x -= m_nn;
x = (x >> m_mm) + (x & m_nn);
}
return x;
}
public:
unsigned int m_mm; /* Bits per symbol */
unsigned int m_nn; /* Symbols per block (= (1<
CReedSolomon::CReedSolomon(unsigned int m) {
assert(m >= 2 && m <= 16);
_INIT_DATA& data = Tab[m - 2];
this->CReedSolomon < T >::CReedSolomon(data.symsize, data.genpoly,
data.fcs, data.prim, data.nroots);
}
template
CReedSolomon::CReedSolomon(_INIT_DATA data) {
this->CReedSolomon < T >::CReedSolomon(data.symsize, data.genpoly,
data.fcs, data.prim, data.nroots);
}
#define A0 m_nn
//#define T unsigned char
template
CReedSolomon::CReedSolomon(
unsigned int symsize,
unsigned int gfpoly,
unsigned int fcr,
unsigned int prim,
unsigned int nroots)
{
int i;
int j;
int sr;
int root;
int iprim;
/* Need version with ints rather than chars */
assert(symsize <= 8 * sizeof(T));
assert(fcr < (1 << symsize));
assert(prim != 0 && prim < (1 << symsize));
/* Can't have more roots than symbol values! */
assert(nroots < (1 << symsize));
m_mm = symsize;
m_nn = (1 << symsize) - 1;
m_alpha_to = new T[m_nn + 1];
m_index_of = new T[m_nn + 1];
/* Generate Galois field lookup tables */
m_index_of[0] = A0; /* log(zero) = -inf */
m_alpha_to[A0] = 0; /* alpha**-inf = 0 */
sr = 1;
for (i = 0; i < m_nn; i++)
{
m_index_of[sr] = i;
m_alpha_to[i] = sr;
sr <<= 1;
if (sr & (1 << symsize)) sr ^= gfpoly;
sr &= m_nn;
}
/* field generator polynomial is not primitive! */
assert(sr == 1);
/* Form RS code generator polynomial from its roots */
m_genpoly = new T[nroots + 1];
m_fcr = fcr;
m_prim = prim;
m_nroots = nroots;
/* Find m_prim-th root of 1, used in decoding */
for (iprim = 1; (iprim % prim) != 0; iprim += m_nn);
m_iprim = iprim / prim;
m_genpoly[0] = 1;
for (i = 0, root = fcr * prim; i < nroots; i++, root += prim)
{
m_genpoly[i + 1] = 1;
/* Multiply m_genpoly[] by @**(root + x) */
for (j = i; j > 0; j--)
{
if (m_genpoly[j] != 0)
{
m_genpoly[j] = m_genpoly[j - 1] ^ m_alpha_to[modnn(m_index_of[m_genpoly[j]] + root)];
}
else
m_genpoly[j] = m_genpoly[j - 1];
}
/* m_genpoly[0] can never be zero */
m_genpoly[0] = m_alpha_to[modnn(m_index_of[m_genpoly[0]] + root)];
}
/* convert m_genpoly[] to index form for quicker encoding */
for (i = 0; i <= nroots; i++)
m_genpoly[i] = m_index_of[m_genpoly[i]];
}
template
CReedSolomon::~CReedSolomon() {
delete m_alpha_to;
delete m_index_of;
delete m_genpoly;
}
template
void CReedSolomon::Encode(T* data, T* bb) {
unsigned int i;
unsigned int j;
T feedback;
memset(bb, 0, m_nroots * sizeof(T));
for (i = 0; i < m_nn - m_nroots; i++)
{
feedback = m_index_of[data[i] ^ bb[0]];
if (feedback != A0)
{ /* feedback term is non-zero */
#ifdef UNNORMALIZED
/* This line is unnecessary when m_genpoly[m_nroots] is unity, as it must
* always be for the polynomials constructed by init_rs()
*/
feedback = modnn(m_nn - m_genpoly[m_nroots] + feedback);
#endif
for (j = 1; j < m_nroots; j++)
{
bb[j] ^= m_alpha_to[modnn(feedback + m_genpoly[m_nroots - j])];
}
}
/* Shift */
memmove(&bb[0], &bb[1], sizeof(T) * (m_nroots - 1));
if (feedback != A0)
{
bb[m_nroots - 1] = m_alpha_to[modnn(feedback + m_genpoly[0])];
}
else
bb[m_nroots - 1] = 0;
}
}
template
int CReedSolomon::Decode(T* data, T* eras_pos, int no_eras) {
int deg_lambda;
int el;
int deg_omega;
int i;
int j;
int r;
int k;
T u;
T q;
T tmp;
T num1;
T num2;
T den;
T discr_r;
/* Err+Eras Locator poly * and syndrome poly */
vector < T > lambda(m_nroots + 1), s(m_nroots);
vector < T > b(m_nroots + 1), t(m_nroots + 1), omega(m_nroots + 1);
vector root(m_nroots);
vector reg(m_nroots + 1);
vector loc(m_nroots);
int syn_error;
int count;
/* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
for (i = 0; i < m_nroots; i++)
s[i] = data[0];
for (j = 1; j < m_nn; j++)
{
for (i = 0; i < m_nroots; i++)
{
if (s[i] == 0)
{
s[i] = data[j];
}
else
{
s[i] = data[j] ^ m_alpha_to[modnn(m_index_of[s[i]] + (m_fcr + i) * m_prim)];
}
}
}
/* Convert syndromes to index form, checking for nonzero condition */
syn_error = 0;
for (i = 0; i < m_nroots; i++)
{
syn_error |= s[i];
s[i] = m_index_of[s[i]];
}
if (!syn_error)
{
/* if syndrome is zero, data[] is a codeword and there are no
* errors to correct. So return data[] unmodified
*/
count = 0;
goto finish;
}
memset(&lambda[1], 0, m_nroots * sizeof(lambda[0]));
lambda[0] = 1;
if (no_eras > 0)
{
/* Init lambda to be the erasure locator polynomial */
lambda[1] = m_alpha_to[modnn(m_prim * (m_nn - 1 - eras_pos[0]))];
for (i = 1; i < no_eras; i++)
{
u = modnn(m_prim * (m_nn - 1 - eras_pos[i]));
for (j = i + 1; j > 0; j--)
{
tmp = m_index_of[lambda[j - 1]];
if (tmp != A0)
lambda[j] ^= m_alpha_to[modnn(u + tmp)];
}
}
}
for (i = 0; i < m_nroots + 1; i++)
b[i] = m_index_of[lambda[i]];
/*
* Begin Berlekamp-Massey algorithm to determine error+erasure
* locator polynomial
*/
r = no_eras;
el = no_eras;
while (++r <= m_nroots)
{ /* 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 - 1] != A0))
{
discr_r ^= m_alpha_to[modnn(m_index_of[lambda[i]] + s[r - i - 1])];
}
}
discr_r = m_index_of[discr_r]; /* Index form */
if (discr_r == A0)
{
/* 2 lines below: B(x) <-- x*B(x) */
//memmove(&b[1],b,m_nroots*sizeof(b[0]));
char_traits < T >::move(b.begin() + 1, b.begin(), b.size() - 1);
b[0] = A0;
}
else
{
/* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
t[0] = lambda[0];
for (i = 0; i < m_nroots; i++)
{
if (b[i] != A0)
{
t[i + 1] = lambda[i + 1] ^ m_alpha_to[modnn(discr_r + b[i])];
}
else
t[i + 1] = lambda[i + 1];
}
if (2 * el <= r + no_eras - 1)
{
el = r + no_eras - el;
/*
* 2 lines below: B(x) <-- inv(discr_r) *
* lambda(x)
*/
for (i = 0; i <= m_nroots; i++)
{
b[i] = (lambda[i] == 0) ? A0 : modnn(m_index_of[lambda[i]] - discr_r + m_nn);
}
}
else
{
/* 2 lines below: B(x) <-- x*B(x) */
//memmove(&b[1],b,m_nroots*sizeof(b[0]));
char_traits < T >::move(b.begin() + 1, b.begin(), b.size() - 1);
b[0] = A0;
}
// memcpy(lambda,t,(m_nroots+1)*sizeof(t[0]));
char_traits < T >::copy(lambda.begin(), t.begin(), t.size());
}
}
/* Convert lambda to index form and compute deg(lambda(x)) */
deg_lambda = 0;
for (i = 0; i < m_nroots + 1; i++)
{
lambda[i] = m_index_of[lambda[i]];
if (lambda[i] != A0)
deg_lambda = i;
}
/* Find roots of the error+erasure locator polynomial by Chien search */
// memcpy(®[1],&lambda[1],m_nroots*sizeof(reg[0]));
char_traits < T >::copy(reg.begin() + 1, lambda.begin() + 1, reg.size() - 1);
count = 0; /* Number of roots of lambda(x) */
for (i = 1, k = m_iprim - 1; i <= m_nn; i++, k = modnn(k + m_iprim))
{
q = 1; /* lambda[0] is always 0 */
for (j = deg_lambda; j > 0; j--)
{
if (reg[j] != A0)
{
reg[j] = modnn(reg[j] + j);
q ^= m_alpha_to[reg[j]];
}
}
if (q != 0) continue; /* Not a root */
/* store root (index-form) and error location number */
root[count] = i;
loc[count] = k;
/* If we've already found max possible roots,
* abort the search to save time
*/
if (++count == deg_lambda) break;
}
if (deg_lambda != count)
{
/*
* deg(lambda) unequal to number of roots => uncorrectable
* error detected
*/
count = -1;
goto finish;
}
/*
* Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
* x**m_nroots). in index form. Also find deg(omega).
*/
deg_omega = 0;
for (i = 0; i < m_nroots; i++)
{
tmp = 0;
j = (deg_lambda < i) ? deg_lambda : i;
for (; j >= 0; j--)
{
if ((s[i - j] != A0) && (lambda[j] != A0))
{
tmp ^= m_alpha_to[modnn(s[i - j] + lambda[j])];
}
}
if (tmp != 0) deg_omega = i;
omega[i] = m_index_of[tmp];
}
omega[m_nroots] = A0;
/*
* Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
* inv(X(l))**(m_fcr-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 ^= m_alpha_to[modnn(omega[i] + i * root[j])];
}
}
num2 = m_alpha_to[modnn(root[j] * (m_fcr - 1) + m_nn)];
den = 0;
/* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
for (i = _cpp_min((unsigned int) deg_lambda, m_nroots - 1) &~1; i >= 0;
i -= 2)
{
if (lambda[i + 1] != A0)
{
den ^= m_alpha_to[modnn(lambda[i + 1] + i * root[j])];
}
}
if (den == 0)
{
count = -1;
goto finish;
}
/* Apply error to data */
if (num1 != 0)
{
data[loc[j]] ^= m_alpha_to[modnn(m_index_of[num1] + m_index_of[num2] + m_nn - m_index_of[den])];
}
}
finish:
if (eras_pos != NULL)
{
for (i = 0; i < count; i++) eras_pos[i] = loc[i];
}
return count;
}
#endif // !defined(AFX_REEDSOLOMON_H__7A4D7FA4_C47A_4880_9E03_833194109EA0__INCLUDED_)