Cleanups to the Modem and Reed-Solomon code.

This commit is contained in:
Jonathan Naylor 2016-09-27 07:56:00 +01:00
parent 7447eee328
commit 018b6e4dec
3 changed files with 112 additions and 79 deletions

View file

@ -562,7 +562,7 @@ void CModem::clock(unsigned int ms)
m_txP25Data.getData(m_buffer, len); m_txP25Data.getData(m_buffer, len);
if (m_debug) { if (m_debug) {
if (m_buffer[3U] == MMDVM_P25_HDR) if (m_buffer[2U] == MMDVM_P25_HDR)
CUtils::dump(1U, "TX P25 HDR", m_buffer, len); CUtils::dump(1U, "TX P25 HDR", m_buffer, len);
else else
CUtils::dump(1U, "TX P25 LDU", m_buffer, len); CUtils::dump(1U, "TX P25 LDU", m_buffer, len);

13
RS.cpp
View file

@ -17,6 +17,7 @@
*/ */
#include "RS.h" #include "RS.h"
#include "Log.h"
#include <cstdio> #include <cstdio>
#include <cassert> #include <cassert>
@ -91,8 +92,10 @@ bool CRS362017::decode(unsigned char* data)
// Now we can call decode on the base class // Now we can call decode on the base class
int irrecoverable_errors = CReedSolomon63<8>::decode(input, output); int irrecoverable_errors = CReedSolomon63<8>::decode(input, output);
if (irrecoverable_errors != 0) if (irrecoverable_errors != 0) {
LogWarning("Unrecoverable errors in the RS(36,20,17) code");
return false; return false;
}
// Convert it back to binary and put it into hex_data. // Convert it back to binary and put it into hex_data.
offset = 0U; offset = 0U;
@ -168,8 +171,10 @@ bool CRS241213::decode(unsigned char* data)
// Now we can call decode on the base class // Now we can call decode on the base class
int irrecoverable_errors = CReedSolomon63<6>::decode(input, output); int irrecoverable_errors = CReedSolomon63<6>::decode(input, output);
if (irrecoverable_errors != 0) if (irrecoverable_errors != 0) {
LogWarning("Unrecoverable errors in the RS(24,12,13) code");
return false; return false;
}
// Convert it back to binary and put it into hex_data. // Convert it back to binary and put it into hex_data.
offset = 0U; offset = 0U;
@ -245,8 +250,10 @@ bool CRS24169::decode(unsigned char* data)
// Now we can call decode on the base class // Now we can call decode on the base class
int irrecoverable_errors = CReedSolomon63<4>::decode(input, output); int irrecoverable_errors = CReedSolomon63<4>::decode(input, output);
if (irrecoverable_errors != 0) if (irrecoverable_errors != 0) {
LogWarning("Unrecoverable errors in the RS(24,16,9) code");
return false; return false;
}
// Convert it back to binary and put it into hex_data. // Convert it back to binary and put it into hex_data.
offset = 0U; offset = 0U;

176
RS.h
View file

@ -61,79 +61,98 @@ Simon Rockliff, 26th June 1991
template<int TT> class CReedSolomon63 template<int TT> class CReedSolomon63
{ {
private: private:
static const int MM = 6; /* RS code over GF(2**mm) */ static const int MM = 6; // RS code over GF(2**mm)
static const int NN = 63; /* nn=2**mm -1 length of codeword */ static const int NN = 63; // nn=2**mm -1 length of codeword
//int tt; /* number of errors that can be corrected */ //int tt; number of errors that can be corrected
//int kk; /* kk = nn-2*tt */ //int kk; kk = nn-2*tt
static const int KK = NN - 2 * TT; static const int KK = NN - 2 * TT;
// distance = nn-kk+1 = 2*tt+1 // distance = nn-kk+1 = 2*tt+1
int* alpha_to;
int* index_of;
int* gg; int* gg;
void generate_gf(int* generator_polinomial) // generate GF(2**mm) from the irreducible polynomial p(X) in pp[0]..pp[mm]
/* 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;
lookup tables: index->polynomial form alpha_to[] contains j=alpha**i; // polynomial form -> index form index_of[j=alpha**i] = i
polynomial form -> index form index_of[j=alpha**i] = i // alpha=2 is the primitive element of GF(2**mm)
alpha=2 is the primitive element of GF(2**mm) void generate_gf(const int* generator_polinomial)
*/
{ {
register int i, mask; register int i, mask;
mask = 1; mask = 1;
alpha_to[MM] = 0; alpha_to[MM] = 0;
for (i = 0; i < MM; i++) { for (i = 0; i < MM; i++) {
alpha_to[i] = mask; alpha_to[i] = mask;
index_of[alpha_to[i]] = i; index_of[alpha_to[i]] = i;
if (generator_polinomial[i] != 0) if (generator_polinomial[i] != 0)
alpha_to[MM] ^= mask; alpha_to[MM] ^= mask;
mask <<= 1; mask <<= 1;
} }
index_of[alpha_to[MM]] = MM; index_of[alpha_to[MM]] = MM;
mask >>= 1; mask >>= 1;
for (i = MM + 1; i < NN; i++) { for (i = MM + 1; i < NN; i++) {
if (alpha_to[i - 1] >= mask) if (alpha_to[i - 1] >= mask)
alpha_to[i] = alpha_to[MM] ^ ((alpha_to[i - 1] ^ mask) << 1); alpha_to[i] = alpha_to[MM] ^ ((alpha_to[i - 1] ^ mask) << 1);
else else
alpha_to[i] = alpha_to[i - 1] << 1; alpha_to[i] = alpha_to[i - 1] << 1;
index_of[alpha_to[i]] = i; index_of[alpha_to[i]] = i;
} }
index_of[0] = -1; index_of[0] = -1;
} }
// 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
void gen_poly() 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; register int i, j;
gg[0] = 2; /* primitive element alpha = 2 for GF(2**mm) */ gg[0] = 2; // primitive element alpha = 2 for GF(2**mm)
gg[1] = 1; /* g(x) = (X+alpha) initially */ gg[1] = 1; // g(x) = (X+alpha) initially
for (i = 2; i <= NN - KK; i++) { for (i = 2; i <= NN - KK; i++) {
gg[i] = 1; gg[i] = 1;
for (j = i - 1; j > 0; j--) for (j = i - 1; j > 0; j--)
if (gg[j] != 0) if (gg[j] != 0)
gg[j] = gg[j - 1] ^ alpha_to[(index_of[gg[j]] + i) % NN]; gg[j] = gg[j - 1] ^ alpha_to[(index_of[gg[j]] + i) % NN];
else else
gg[j] = gg[j - 1]; gg[j] = gg[j - 1];
gg[0] = alpha_to[(index_of[gg[0]] + i) % NN]; /* gg[0] can never be zero */
gg[0] = alpha_to[(index_of[gg[0]] + i) % NN]; // gg[0] can never be zero
} }
/* convert gg[] to index form for quicker encoding */
// convert gg[] to index form for quicker encoding
for (i = 0; i <= NN - KK; i++) for (i = 0; i <= NN - KK; i++)
gg[i] = index_of[gg[i]]; gg[i] = index_of[gg[i]];
} }
public: protected:
int* alpha_to;
int* index_of;
CReedSolomon63() CReedSolomon63()
{ {
alpha_to = new int[NN + 1]; alpha_to = new int[NN + 1];
index_of = new int[NN + 1]; index_of = new int[NN + 1];
gg = new int[NN - KK + 1]; gg = new int[NN - KK + 1];
for (unsigned int i = 0U; i < (NN + 1); i++) {
alpha_to[i] = 0;
index_of[i] = 0;
}
for (unsigned int i = 0U; i < (NN - KK + 1); i++)
gg[i] = 0;
// Polynom used in P25 is alpha**6+alpha+1 // Polynom used in P25 is alpha**6+alpha+1
int generator_polinomial[] = { 1, 1, 0, 0, 0, 0, 1 }; /* specify irreducible polynomial coeffts */ const int generator_polinomial[] = {1, 1, 0, 0, 0, 0, 1}; // specify irreducible polynomial coeffts
generate_gf(generator_polinomial); generate_gf(generator_polinomial);
@ -147,19 +166,20 @@ public:
delete[] alpha_to; delete[] alpha_to;
} }
// 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)
void encode(const int* data, int* bb) void encode(const int* data, int* bb)
/* 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; register int i, j;
int feedback; int feedback;
for (i = 0; i < NN - KK; i++) for (i = 0; i < NN - KK; i++)
bb[i] = 0; bb[i] = 0;
for (i = KK - 1; i >= 0; i--) { for (i = KK - 1; i >= 0; i--) {
feedback = index_of[data[i] ^ bb[NN - KK - 1]]; feedback = index_of[data[i] ^ bb[NN - KK - 1]];
if (feedback != -1) { if (feedback != -1) {
@ -168,41 +188,40 @@ public:
bb[j] = bb[j - 1] ^ alpha_to[(gg[j] + feedback) % NN]; bb[j] = bb[j - 1] ^ alpha_to[(gg[j] + feedback) % NN];
else else
bb[j] = bb[j - 1]; bb[j] = bb[j - 1];
bb[0] = alpha_to[(gg[0] + feedback) % NN]; bb[0] = alpha_to[(gg[0] + feedback) % NN];
} } else {
else {
for (j = NN - KK - 1; j > 0; j--) for (j = NN - KK - 1; j > 0; j--)
bb[j] = bb[j - 1]; bb[j] = bb[j - 1];
bb[0] = 0; bb[0] = 0;
} }
} }
} }
/* assume we have received bits grouped into mm-bit symbols in recd[i],
i=0..(nn-1), and recd[i] is polynomial form.
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. */
int decode(const int* input, int* recd) int decode(const int* input, int* recd)
/* assume we have received bits grouped into mm-bit symbols in recd[i],
i=0..(nn-1), and recd[i] is polynomial form.
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; 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 int elp[NN - KK + 2][NN - KK], d[NN - KK + 2], l[NN - KK + 2], u_lu[NN - KK + 2], s[NN - KK + 1];
+ 2], s[NN - KK + 1]; int count = 0, syn_error = 0, root[TT], loc[TT], z[TT + 1], err[NN], reg[TT + 1];
int count = 0, syn_error = 0, root[TT], loc[TT], z[TT + 1], err[NN], reg[TT
+ 1];
int irrecoverable_error = 0; int irrecoverable_error = 0;
@ -212,12 +231,15 @@ public:
/* first form the syndromes */ /* first form the syndromes */
for (i = 1; i <= NN - KK; i++) { for (i = 1; i <= NN - KK; i++) {
s[i] = 0; s[i] = 0;
for (j = 0; j < NN; j++) for (j = 0; j < NN; j++)
if (recd[j] != -1) if (recd[j] != -1)
s[i] ^= alpha_to[(recd[j] + i * j) % NN]; /* recd[j] in index form */ s[i] ^= alpha_to[(recd[j] + i * j) % NN]; /* recd[j] in index form */
/* convert syndrome from polynomial form to index form */ /* convert syndrome from polynomial form to index form */
if (s[i] != 0) if (s[i] != 0)
syn_error = 1; /* set flag if non-zero syndrome => error */ syn_error = 1; /* set flag if non-zero syndrome => error */
s[i] = index_of[s[i]]; s[i] = index_of[s[i]];
} }
@ -230,15 +252,18 @@ public:
degree of the elp at that step, and u_l[u] is the difference between 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. step number and the degree of the elp.
*/ */
/* initialise table entries */ /* initialise table entries */
d[0] = 0; /* index form */ d[0] = 0; /* index form */
d[1] = s[1]; /* index form */ d[1] = s[1]; /* index form */
elp[0][0] = 0; /* index form */ elp[0][0] = 0; /* index form */
elp[1][0] = 1; /* polynomial form */ elp[1][0] = 1; /* polynomial form */
for (i = 1; i < NN - KK; i++) { for (i = 1; i < NN - KK; i++) {
elp[0][i] = -1; /* index form */ elp[0][i] = -1; /* index form */
elp[1][i] = 0; /* polynomial form */ elp[1][i] = 0; /* polynomial form */
} }
l[0] = 0; l[0] = 0;
l[1] = 0; l[1] = 0;
u_lu[0] = -1; u_lu[0] = -1;
@ -247,19 +272,21 @@ public:
do { do {
u++; u++;
if (d[u] == -1) { if (d[u] == -1) {
l[u + 1] = l[u]; l[u + 1] = l[u];
for (i = 0; i <= l[u]; i++) { for (i = 0; i <= l[u]; i++) {
elp[u + 1][i] = elp[u][i]; elp[u + 1][i] = elp[u][i];
elp[u][i] = index_of[elp[u][i]]; elp[u][i] = index_of[elp[u][i]];
} }
} } else
else
/* search for words with greatest u_lu[q] for which d[q]!=0 */ /* search for words with greatest u_lu[q] for which d[q]!=0 */
{ {
q = u - 1; q = u - 1;
while ((d[q] == -1) && (q > 0)) while ((d[q] == -1) && (q > 0))
q--; q--;
/* have found first non-zero d[q] */ /* have found first non-zero d[q] */
if (q > 0) { if (q > 0) {
j = q; j = q;
@ -280,15 +307,17 @@ public:
/* form new elp(x) */ /* form new elp(x) */
for (i = 0; i < NN - KK; i++) for (i = 0; i < NN - KK; i++)
elp[u + 1][i] = 0; elp[u + 1][i] = 0;
for (i = 0; i <= l[q]; i++) for (i = 0; i <= l[q]; i++)
if (elp[q][i] != -1) if (elp[q][i] != -1)
elp[u + 1][i + u - q] = alpha_to[(d[u] + NN - d[q] elp[u + 1][i + u - q] = alpha_to[(d[u] + NN - d[q] + elp[q][i]) % NN];
+ elp[q][i]) % NN];
for (i = 0; i <= l[u]; i++) { for (i = 0; i <= l[u]; i++) {
elp[u + 1][i] ^= elp[u][i]; elp[u + 1][i] ^= elp[u][i];
elp[u][i] = index_of[elp[u][i]]; /*convert old elp value to index*/ elp[u][i] = index_of[elp[u][i]]; /*convert old elp value to index*/
} }
} }
u_lu[u + 1] = u - l[u + 1]; u_lu[u + 1] = u - l[u + 1];
/* form (u+1)th discrepancy */ /* form (u+1)th discrepancy */
@ -298,10 +327,12 @@ public:
d[u + 1] = alpha_to[s[u + 1]]; d[u + 1] = alpha_to[s[u + 1]];
else else
d[u + 1] = 0; d[u + 1] = 0;
for (i = 1; i <= l[u + 1]; i++) for (i = 1; i <= l[u + 1]; i++)
if ((s[u + 1 - i] != -1) && (elp[u + 1][i] != 0)) if ((s[u + 1 - i] != -1) && (elp[u + 1][i] != 0))
d[u + 1] ^= alpha_to[(s[u + 1 - i] d[u + 1] ^= alpha_to[(s[u + 1 - i]
+ index_of[elp[u + 1][i]]) % NN]; + index_of[elp[u + 1][i]]) % NN];
d[u + 1] = index_of[d[u + 1]]; /* put d[u+1] into index form */ d[u + 1] = index_of[d[u + 1]]; /* put d[u+1] into index form */
} }
} while ((u < NN - KK) && (l[u + 1] <= TT)); } while ((u < NN - KK) && (l[u + 1] <= TT));
@ -316,14 +347,18 @@ public:
/* find roots of the error location polynomial */ /* find roots of the error location polynomial */
for (i = 1; i <= l[u]; i++) for (i = 1; i <= l[u]; i++)
reg[i] = elp[u][i]; reg[i] = elp[u][i];
count = 0; count = 0;
for (i = 1; i <= NN; i++) { for (i = 1; i <= NN; i++) {
q = 1; q = 1;
for (j = 1; j <= l[u]; j++) for (j = 1; j <= l[u]; j++)
if (reg[j] != -1) { if (reg[j] != -1) {
reg[j] = (reg[j] + j) % NN; reg[j] = (reg[j] + j) % NN;
q ^= alpha_to[reg[j]]; q ^= alpha_to[reg[j]];
}; };
if (!q) /* store root and error location number indices */ if (!q) /* store root and error location number indices */
{ {
root[count] = i; root[count] = i;
@ -345,9 +380,11 @@ public:
z[i] = alpha_to[elp[u][i]]; z[i] = alpha_to[elp[u][i]];
else else
z[i] = 0; z[i] = 0;
for (j = 1; j < i; j++) for (j = 1; j < i; j++)
if ((s[j] != -1) && (elp[u][i - j] != -1)) if ((s[j] != -1) && (elp[u][i - j] != -1))
z[i] ^= alpha_to[(elp[u][i - j] + s[j]) % NN]; z[i] ^= alpha_to[(elp[u][i - j] + s[j]) % NN];
z[i] = index_of[z[i]]; /* put into index form */ z[i] = index_of[z[i]]; /* put into index form */
}; };
@ -359,38 +396,39 @@ public:
else else
recd[i] = 0; recd[i] = 0;
} }
for (i = 0; i < l[u]; i++) /* compute numerator of error term first */ for (i = 0; i < l[u]; i++) /* compute numerator of error term first */
{ {
err[loc[i]] = 1; /* accounts for z[0] */ err[loc[i]] = 1; /* accounts for z[0] */
for (j = 1; j <= l[u]; j++) for (j = 1; j <= l[u]; j++)
if (z[j] != -1) if (z[j] != -1)
err[loc[i]] ^= alpha_to[(z[j] + j * root[i]) % NN]; err[loc[i]] ^= alpha_to[(z[j] + j * root[i]) % NN];
if (err[loc[i]] != 0) { if (err[loc[i]] != 0) {
err[loc[i]] = index_of[err[loc[i]]]; err[loc[i]] = index_of[err[loc[i]]];
q = 0; /* form denominator of error term */ q = 0; /* form denominator of error term */
for (j = 0; j < l[u]; j++) for (j = 0; j < l[u]; j++)
if (j != i) if (j != i)
q += index_of[1 q += index_of[1 ^ alpha_to[(loc[j] + root[i]) % NN]];
^ alpha_to[(loc[j] + root[i]) % NN]];
q = q % NN; q = q % NN;
err[loc[i]] = alpha_to[(err[loc[i]] - q + NN) % 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 */ recd[loc[i]] ^= err[loc[i]]; /*recd[i] must be in polynomial form */
} }
} }
} } else {
else {
/* no. roots != degree of elp => >tt errors and cannot solve */ /* no. roots != degree of elp => >tt errors and cannot solve */
irrecoverable_error = 1; irrecoverable_error = 1;
} }
} } else {
else {
/* elp has degree >tt hence cannot solve */ /* elp has degree >tt hence cannot solve */
irrecoverable_error = 1; irrecoverable_error = 1;
} }
} } else {
else {
/* no non-zero syndromes => no errors: output received codeword */ /* no non-zero syndromes => no errors: output received codeword */
for (i = 0; i < NN; i++) for (i = 0; i < NN; i++)
if (recd[i] != -1) /* convert recd[] to polynomial form */ if (recd[i] != -1) /* convert recd[] to polynomial form */
@ -411,10 +449,6 @@ public:
} }
}; };
/**
* Convenience class that does a Reed-Solomon (36,20,17) error correction adapting input and output to
* the DSD data format: hex words packed as char arrays.
*/
class CRS362017 : public CReedSolomon63<8> class CRS362017 : public CReedSolomon63<8>
{ {
public: public:
@ -428,10 +462,6 @@ public:
private: private:
}; };
/**
* Convenience class that does a Reed-Solomon (24,12,13) error correction adapting input and output to
* the DSD data format: hex words packed as char arrays.
*/
class CRS241213 : public CReedSolomon63<6> class CRS241213 : public CReedSolomon63<6>
{ {
public: public:
@ -445,10 +475,6 @@ public:
private: private:
}; };
/**
* Convenience class that does a Reed-Solomon (24,16,9) error correction adapting input and output to
* the DSD data format: hex words packed as char arrays.
*/
class CRS24169 : public CReedSolomon63<4> class CRS24169 : public CReedSolomon63<4>
{ {
public: public: