SDRAngel  4.11.5
Developer docs for <a href="https://github.com/f4exb/sdrangel">SDRangel<\a>, an Open Source Qt5 / OpenGL 3.0+ SDR and signal analyzer frontend to various hardware.
rs.h
Go to the documentation of this file.
1 // This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
2 // See the toplevel README for more information.
3 //
4 // This program is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation, either version 3 of the License, or
7 // (at your option) any later version.
8 //
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with this program. If not, see <http://www.gnu.org/licenses/>.
16 
17 #ifndef LEANSDR_RS_H
18 #define LEANSDR_RS_H
19 
20 #include "leansdr/math.h"
21 
22 #define DEBUG_RS 0
23 
24 namespace leansdr
25 {
26 
27 // Finite group GF(2^N).
28 
29 // GF(2) is the ring ({0,1},+,*).
30 // GF(2)[X] is the ring of polynomials with coefficients in GF(2).
31 // P(X) is an irreducible polynomial of GF(2)[X].
32 // N is the degree of P(x).
33 // P is the bitfield representation of P(X), with degree 0 at LSB.
34 // GF(2)[X]/(P) is GF(2)[X] modulo P(X).
35 // (GF(2)[X]/(P), +) is a group with 2^N elements.
36 // (GF(2)[X]/(P)*, *) is a group with 2^N-1 elements.
37 // (GF(2)[X]/(P), +, *) is a field with 2^N elements, noted GF(2^N).
38 // Te is a C++ integer type for representing elements of GF(2^N).
39 // "0" is 0
40 // "1" is 1
41 // "2" is X
42 // "3" is X+1
43 // "4" is X^2
44 // Tp is a C++ integer type for representing P(X) (1 bit larger than Te).
45 // ALPHA is a primitive element of GF(2^N). Usually "2"=[X] is chosen.
46 
47 template <typename Te, typename Tp, Tp P, int N, Te ALPHA>
48 struct gf2x_p
49 {
51  {
52  if (ALPHA != 2)
53  fail("alpha!=2 not implemented");
54  // Precompute log and exp tables.
55  Tp alpha_i = 1;
56  for (Tp i = 0; i < (1 << N); ++i)
57  {
58  lut_exp[i] = alpha_i;
59  lut_exp[((1 << N) - 1) + i] = alpha_i;
60  lut_log[alpha_i] = i;
61  alpha_i <<= 1; // Multiply by alpha=[X] i.e. increase degrees
62  if (alpha_i & (1 << N))
63  alpha_i ^= P; // Modulo P iteratively
64  }
65  }
66  static const Te alpha = ALPHA;
67  inline Te add(Te x, Te y) { return x ^ y; } // Addition modulo 2
68  inline Te sub(Te x, Te y) { return x ^ y; } // Subtraction modulo 2
69  inline Te mul(Te x, Te y)
70  {
71  if (!x || !y)
72  return 0;
73  return lut_exp[lut_log[x] + lut_log[y]];
74  }
75  inline Te div(Te x, Te y)
76  {
77  //if ( ! y ) fail("div"); // TODO
78  if (!x)
79  return 0;
80  return lut_exp[lut_log[x] + ((1 << N) - 1) - lut_log[y]];
81  }
82  inline Te inv(Te x)
83  {
84  // if ( ! x ) fail("inv");
85  return lut_exp[((1 << N) - 1) - lut_log[x]];
86  }
87  inline Te exp(Te x) { return lut_exp[x]; }
88  inline Te log(Te x) { return lut_log[x]; }
89 
90  private:
91  Te lut_exp[(1 << N) * 2]; // Wrap to avoid indexing modulo 2^N-1
92  Te lut_log[1 << N];
93 };
94 
95 // Reed-Solomon for RS(204,188) shortened from RS(255,239).
96 
97 struct rs_engine
98 {
99  // EN 300 421, section 4.4.2, Field Generator Polynomial
100  // p(X) = X^8 + X^4 + X^3 + X^2 + 1
102 
103  u8 G[17]; // { G_16, ..., G_0 }
104 
106  {
107  // EN 300 421, section 4.4.2, Code Generator Polynomial
108  // G(X) = (X-alpha^0)*...*(X-alpha^15)
109  for (int i = 0; i <= 16; ++i)
110  G[i] = (i == 16) ? 1 : 0; // Init G=1
111  for (int d = 0; d < 16; ++d)
112  {
113  // Multiply by (X-alpha^d)
114  // G := X*G - alpha^d*G
115  for (int i = 0; i <= 16; ++i)
116  G[i] = gf.sub((i == 16) ? 0 : G[i + 1], gf.mul(gf.exp(d), G[i]));
117  }
118 #if DEBUG_RS
119  fprintf(stderr, "RS generator:");
120  for (int i = 0; i <= 16; ++i)
121  fprintf(stderr, " %02x", G[i]);
122  fprintf(stderr, "\n");
123 #endif
124  }
125 
126  // RS-encoded messages are interpreted as coefficients in
127  // GF(256) of a polynomial P.
128  // The syndromes are synd[i] = P(alpha^i).
129  // By convention coefficients are listed by decreasing degree here,
130  // so we can evaluate syndromes of the shortened code without
131  // prepending with 51 zeroes.
132  bool syndromes(const u8 *poly, u8 *synd)
133  {
134  bool corrupted = false;
135  for (int i = 0; i < 16; ++i)
136  {
137  synd[i] = eval_poly_rev(poly, 204, gf.exp(i));
138  if (synd[i])
139  corrupted = true;
140  }
141  return corrupted;
142  }
143  u8 eval_poly_rev(const u8 *poly, int n, u8 x)
144  {
145  // poly[0]*x^(n-1) + .. + poly[n-1]*x^0 with Hörner method.
146  u8 acc = 0;
147  for (int i = 0; i < n; ++i)
148  acc = gf.add(gf.mul(acc, x), poly[i]);
149  return acc;
150  }
151 
152  // Evaluation with coefficients listed by increasing degree.
153  u8 eval_poly(const u8 *poly, int deg, u8 x)
154  {
155  // poly[0]*x^0 + .. + poly[deg]*x^deg with Hörner method.
156  u8 acc = 0;
157  for (; deg >= 0; --deg)
158  acc = gf.add(gf.mul(acc, x), poly[deg]);
159  return acc;
160  }
161 
162  // Append parity symbols
163 
164  void encode(u8 msg[204])
165  {
166  // TBD Avoid copying
167  u8 p[204];
168  memcpy(p, msg, 188);
169  memset(p + 188, 0, 16);
170  // p = msg*X^16
171 #if DEBUG_RS
172  fprintf(stderr, "uncoded:");
173  for (int i = 0; i < 204; ++i)
174  fprintf(stderr, " %d", p[i]);
175  fprintf(stderr, "\n");
176 #endif
177  // Compute remainder modulo G
178  for (int d = 0; d < 188; ++d)
179  {
180  // Clear monomial of degree d
181  if (!p[d])
182  continue;
183  u8 k = gf.div(p[d], G[0]);
184  // p(X) := p(X) - k*G(X)*X^(188-d)
185  for (int i = 0; i <= 16; ++i)
186  p[d + i] = gf.sub(p[d + i], gf.mul(k, G[i]));
187  }
188 #if DEBUG_RS
189  fprintf(stderr, "coded:");
190  for (int i = 0; i < 204; ++i)
191  fprintf(stderr, " %d", p[i]);
192  fprintf(stderr, "\n");
193 #endif
194  memcpy(msg + 188, p + 188, 16);
195  }
196 
197  // Try to fix errors in pout[].
198  // If pin[] is provided, errors will be fixed in the original
199  // message too and syndromes will be updated.
200 
201  bool correct(u8 synd[16], u8 pout[188],
202  u8 pin[204] = NULL, int *bits_corrected = NULL)
203  {
204  // Berlekamp - Massey
205  // http://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm#Code_sample
206  u8 C[16] = {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; // Max degree is L
207  u8 B[16] = {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
208  int L = 0;
209  int m = 1;
210  u8 b = 1;
211  for (int n = 0; n < 16; ++n)
212  {
213  u8 d = synd[n];
214  for (int i = 1; i <= L; ++i)
215  d ^= gf.mul(C[i], synd[n - i]);
216  if (!d)
217  {
218  ++m;
219  }
220  else if (2 * L <= n)
221  {
222  u8 T[16];
223  memcpy(T, C, sizeof(T));
224  for (int i = 0; i < 16 - m; ++i)
225  C[m + i] ^= gf.mul(d, gf.mul(gf.inv(b), B[i]));
226  L = n + 1 - L;
227  memcpy(B, T, sizeof(B));
228  b = d;
229  m = 1;
230  }
231  else
232  {
233  for (int i = 0; i < 16 - m; ++i)
234  C[m + i] ^= gf.mul(d, gf.mul(gf.inv(b), B[i]));
235  ++m;
236  }
237  }
238  // L is the number of errors
239  // C of degree L is now the error locator polynomial (Lambda)
240 #if DEBUG_RS
241  fprintf(stderr, "[L=%d C=", L);
242  for (int i = 0; i < 16; ++i)
243  fprintf(stderr, " %d", C[i]);
244  fprintf(stderr, "]\n");
245  fprintf(stderr, "[S=");
246  for (int i = 0; i < 16; ++i)
247  fprintf(stderr, " %d", synd[i]);
248  fprintf(stderr, "]\n");
249 #endif
250 
251  // Forney
252  // http://en.wikipedia.org/wiki/Forney_algorithm (2t=16)
253 
254  // Compute Omega
255  u8 omega[16];
256  memset(omega, 0, sizeof(omega));
257  // TODO loops
258  for (int i = 0; i < 16; ++i)
259  for (int j = 0; j < 16; ++j)
260  if (i + j < 16)
261  omega[i + j] ^= gf.mul(synd[i], C[j]);
262 #if DEBUG_RS
263  fprintf(stderr, "omega=");
264  for (int i = 0; i < 16; ++i)
265  fprintf(stderr, " %d", omega[i]);
266  fprintf(stderr, "\n");
267 #endif
268 
269  // Compute Lambda'
270  u8 Cprime[15];
271  for (int i = 0; i < 15; ++i)
272  Cprime[i] = (i & 1) ? 0 : C[i + 1];
273 #if DEBUG_RS
274  fprintf(stderr, "Cprime=");
275  for (int i = 0; i < 15; ++i)
276  fprintf(stderr, " %d", Cprime[i]);
277  fprintf(stderr, "\n");
278 #endif
279 
280  // Find zeroes of C by exhaustive search?
281  // TODO Chien method
282  int roots_found = 0;
283  for (int i = 0; i < 255; ++i)
284  {
285  u8 r = gf.exp(i); // Candidate root alpha^0..alpha^254
286  u8 v = eval_poly(C, L, r);
287  if (!v)
288  {
289  // r is a root X_k^-1 of the error locator polynomial.
290  u8 xk = gf.inv(r);
291  int loc = (255 - i) % 255; // == log(xk)
292 #if DEBUG_RS
293  fprintf(stderr, "found root=%d, inv=%d, loc=%d\n", r, xk, loc);
294 #endif
295  if (loc < 204)
296  {
297  // Evaluate e_k
298  u8 num = gf.mul(xk, eval_poly(omega, L, r));
299  u8 den = eval_poly(Cprime, 14, r);
300  u8 e = gf.div(num, den);
301  // Subtract e from coefficient of degree loc.
302  // Note: Coeffients listed by decreasing degree in pin[] and pout[].
303  if (bits_corrected)
304  *bits_corrected += hamming_weight(e);
305  if (loc >= 16)
306  pout[203 - loc] ^= e;
307  if (pin)
308  pin[203 - loc] ^= e;
309  }
310  if (++roots_found == L)
311  break;
312  }
313  }
314  if (pin)
315  return syndromes(pin, synd);
316  else
317  return false;
318  }
319 };
320 
321 } // namespace leansdr
322 
323 #endif // LEANSDR_RS_H
void encode(u8 msg[204])
Definition: rs.h:164
gf2x_p< unsigned char, unsigned short, 0x11d, 8, 2 > gf
Definition: rs.h:101
Te log(Te x)
Definition: rs.h:88
Te mul(Te x, Te y)
Definition: rs.h:69
Te add(Te x, Te y)
Definition: rs.h:67
unsigned char u8
Definition: framework.h:453
int hamming_weight(uint8_t x)
Definition: math.cpp:6
bool correct(u8 synd[16], u8 pout[188], u8 pin[204]=NULL, int *bits_corrected=NULL)
Definition: rs.h:201
Te sub(Te x, Te y)
Definition: rs.h:68
Te lut_exp[(1<< N) *2]
Definition: rs.h:91
Te inv(Te x)
Definition: rs.h:82
static const Te alpha
Definition: rs.h:66
void fail(const char *s)
Definition: framework.cpp:11
u8 eval_poly(const u8 *poly, int deg, u8 x)
Definition: rs.h:153
u8 eval_poly_rev(const u8 *poly, int n, u8 x)
Definition: rs.h:143
int32_t i
Definition: decimators.h:244
Te div(Te x, Te y)
Definition: rs.h:75
Te exp(Te x)
Definition: rs.h:87
bool syndromes(const u8 *poly, u8 *synd)
Definition: rs.h:132
Te lut_log[1<< N]
Definition: rs.h:92
gf2x_p()
Definition: rs.h:50