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.
bch.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_BCH_H
18 #define LEANSDR_BCH_H
19 
20 #include "leansdr/discrmath.h"
21 
22 namespace leansdr
23 {
24 
25 // Interface to hide the template parameters
26 
28 {
29  virtual void encode(const uint8_t *msg, size_t msgbytes, uint8_t *out) = 0;
30  virtual int decode(uint8_t *cw, size_t cwbytes) = 0;
31 }; // bch_interface
32 
33 // BCH error correction.
34 // T: Unsigned type for packing binary polynomials.
35 // N: Number of parity bits.
36 // NP: Width of the polynomials supplied.
37 // DP: Actual degree of the minimum polynomials (all must be same).
38 // TGF: Unsigned type for syndromes (must be wider than DP).
39 // GFTRUNCGEN: Generator polynomial for GF(2^DP), with X^DP omitted.
40 
41 template <typename T, int N, int NP, int DP, typename TGF, int GFTRUNCGEN>
43 {
44  bch_engine(const bitvect<T, NP> *polys, int _npolys)
45  : npolys(_npolys)
46  {
47  // Build the generator polynomial (product of polys[]).
48  g = 1;
49  for (int i = 0; i < npolys; ++i)
50  g = g * polys[i];
51  // Convert the polynomials to truncated representation
52  // (with X^DP omitted) for use with divmod().
53  truncpolys = new bitvect<T, DP>[npolys];
54  for (int i = 0; i < npolys; ++i)
55  truncpolys[i].copy(polys[i]);
56  // Check which polynomial contains each root.
57  // Note: The DVB-S2 polynomials are numbered so that
58  // syndpoly[2*i]==i, but we don't use that property.
59  syndpolys = new int[2 * npolys];
60  for (int i = 0; i < 2 * npolys; ++i)
61  {
62  int j;
63  for (j = 0; j < npolys; ++j)
64  if (!eval_poly(truncpolys[j], true, 1 + i))
65  break;
66  if (j == npolys)
67  fail("Bad polynomials/root");
68  syndpolys[i] = j;
69  }
70  }
71 
72  // Generate BCH parity bits.
73 
74  void encode(const uint8_t *msg, size_t msgbytes, uint8_t *out)
75  {
76  bitvect<T, N> parity = shiftdivmod(msg, msgbytes, g);
77  // Output as bytes, coefficient of highest degree first
78  for (int i = N / 8; i--; ++out)
79  *out = parity.v[i / sizeof(T)] >> ((i & (sizeof(T) - 1)) * 8);
80  }
81 
82  // Decode BCH.
83  // Return number of bits corrected, or -1 on failure.
84 
85  int decode(uint8_t *cw, size_t cwbytes)
86  {
87  //again:
88  bool corrupted = false;
89  // Divide by individual polynomials.
90  // TBD Maybe do in parallel, scanning cw only once.
91  bitvect<T, DP> *rem = new bitvect<T, DP>[npolys]; // npolys is not static hence 'bitvect<T, DP> rem[npolys]' does not compile in all compilers
92  for (int j = 0; j < npolys; ++j)
93  {
94  rem[j] = divmod(cw, cwbytes, truncpolys[j]);
95  }
96  // Compute syndromes.
97  TGF *S = new TGF[2 * npolys]; // npolys is not static hence 'TGF S[2 * npolys]' does not compile in all compilers
98  for (int i = 0; i < 2 * npolys; ++i)
99  {
100  // Compute R(alpha^(1+i)), exploiting the fact that
101  // R(x)=Q(x)g_j(X)+rem_j(X) and g_j(alpha^(1+i))=0
102  // for some j that we already determined.
103  // TBD Compute even exponents using conjugates.
104  S[i] = eval_poly(rem[syndpolys[i]], false, 1 + i);
105  if (S[i])
106  corrupted = true;
107  }
108  if (!corrupted)
109  {
110  delete[] S;
111  delete[] rem;
112  return 0;
113  }
114 #if 0
115  fprintf(stderr, "synd:");
116  for ( int i=0; i<2*npolys; ++i ) fprintf(stderr, " %04x", S[i]);
117  fprintf(stderr, "\n");
118 #endif
119 
120  // S_j = R(alpha_j) = 0+E(alpha_j) = sum(l=1..L)((alpha^j)^i_l)
121  // where i_1 .. i_L are the degrees of the non-zero coefficients of E.
122  // S_j = sum(l=1..L)((alpha^i_l)^j) = sum(l=1..L)(X_l^j)
123 
124  // Berlekamp - Massey
125  // http://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm
126  // TBD More efficient to work with logs of syndromes ?
127 
128  int NN = 2 * npolys;
129  TGF *C = new TGF[NN];
130  C[0] = 1;
131  TGF *B = new TGF[NN];
132  B[0] = 1;
133  // TGF C[NN] = { crap code
134  // 1,
135  // },
136  // B[NN] = {
137  // 1,
138  // };
139  int L = 0, m = 1;
140  TGF b = 1;
141  for (int n = 0; n < NN; ++n)
142  {
143  TGF d = S[n];
144  for (int i = 1; i <= L; ++i)
145  d = GF.add(d, GF.mul(C[i], S[n - i]));
146  if (d == 0)
147  ++m;
148  else
149  {
150  TGF d_div_b = GF.mul(d, GF.inv(b));
151  if (2 * L <= n)
152  {
153  TGF *tmp = new TGF[NN]; // replaced crap code
154  memcpy(tmp, C, sizeof(tmp));
155  for (int i = 0; i < NN - m; ++i)
156  C[m + i] = GF.sub(C[m + i], GF.mul(d_div_b, B[i]));
157  L = n + 1 - L;
158  memcpy(B, tmp, sizeof(B));
159  b = d;
160  m = 1;
161  delete[] tmp;
162  }
163  else
164  {
165  for (int i = 0; i < NN - m; ++i)
166  C[m + i] = GF.sub(C[m + i], GF.mul(d_div_b, B[i]));
167  ++m;
168  }
169  }
170  }
171  // L is the number of errors.
172  // C of degree L is the error locator polynomial (Lambda).
173  // C(X) = sum(l=1..L)(1-X_l*X).
174 #if 0
175  fprintf(stderr, "C[%d]=", L);
176  for ( int i=0; i<NN; ++i ) fprintf(stderr, " %04x", C[i]);
177  fprintf(stderr, "\n");
178 #endif
179 
180  // Forney
181  // http://en.wikipedia.org/wiki/Forney_algorithm
182  // Simplified because coefficients are in GF(2).
183 
184  // Find zeroes of C by exhaustive search.
185  // TODO Chien method
186  int roots_found = 0;
187  for (int i = 0; i < (1 << DP) - 1; ++i)
188  {
189  // Candidate root ALPHA^i
190  TGF v = eval_poly(C, L, i);
191  if (!v)
192  {
193  // ALPHA^i is a root of C, i.e. the inverse of an X_l.
194  int loc = (i ? (1 << DP) - 1 - i : 0); // exponent of inverse
195  // Reverse because cw[0..cwbytes-1] is stored MSB first
196  int rloc = cwbytes * 8 - 1 - loc;
197  if (rloc < 0)
198  {
199  // This may happen if the code is used truncated.
200  delete[] C;
201  delete[] B;
202  delete[] S;
203  delete[] rem;
204  return -1;
205  }
206  cw[rloc / 8] ^= 128 >> (rloc & 7);
207  ++roots_found;
208  if (roots_found == L)
209  break;
210  }
211  }
212 
213  delete[] C;
214  delete[] B;
215  delete[] S;
216  delete[] rem;
217 
218  if (roots_found != L)
219  return -1;
220  return L;
221  }
222 
223  private:
224  // Eval a GF(2)[X] polynomial at a power of ALPHA.
225 
226  TGF eval_poly(const bitvect<T, DP> &poly, bool is_trunc, int rootexp)
227  {
228  TGF acc = 0;
229  int re = 0;
230  for (int i = 0; i < DP; ++i)
231  {
232  if (poly[i])
233  acc = GF.add(acc, GF.exp(re));
234  re += rootexp;
235  if (re >= (1 << DP) - 1)
236  re -= (1 << DP) - 1; // mod 2^DP-1 incrementally
237  }
238  if (is_trunc)
239  acc = GF.add(acc, GF.exp(re));
240  return acc;
241  }
242 
243  // Eval a GF(2^16)[X] polynomial at a power of ALPHA.
244 
245  TGF eval_poly(const TGF *poly, int deg, int rootexp)
246  {
247  TGF acc = 0;
248  int re = 0;
249  for (int i = 0; i <= deg; ++i)
250  {
251  acc = GF.add(acc, GF.mul(poly[i], GF.exp(re)));
252  re += rootexp;
253  if (re >= (1 << DP) - 1)
254  re -= (1 << DP) - 1; // mod 2^DP-1 incrementally
255  }
256  return acc;
257  }
258 
260  int npolys;
261  int *syndpolys;
263  // Finite group for syndrome calculations
265 }; // bch_engine
266 
267 } // namespace leansdr
268 
269 #endif // LEANSDR_BCH_H
virtual int decode(uint8_t *cw, size_t cwbytes)=0
int * syndpolys
Definition: bch.h:261
bitvect< T, DP > * truncpolys
Definition: bch.h:259
int decode(uint8_t *cw, size_t cwbytes)
Definition: bch.h:85
bitvect< T, N > g
Definition: bch.h:262
void encode(const uint8_t *msg, size_t msgbytes, uint8_t *out)
Definition: bch.h:74
virtual void encode(const uint8_t *msg, size_t msgbytes, uint8_t *out)=0
unsigned char uint8_t
Definition: rtptypes_win.h:42
void fail(const char *s)
Definition: framework.cpp:11
int32_t i
Definition: decimators.h:244
bitvect< T, N > divmod(const Tm *m, size_t nm, const bitvect< T, N > &p)
Definition: discrmath.h:111
gf2n< TGF, DP, 2, GFTRUNCGEN > GF
Definition: bch.h:264
TGF eval_poly(const bitvect< T, DP > &poly, bool is_trunc, int rootexp)
Definition: bch.h:226
unsigned char parity(uint8_t x)
Definition: math.cpp:27
TGF eval_poly(const TGF *poly, int deg, int rootexp)
Definition: bch.h:245
bch_engine(const bitvect< T, NP > *polys, int _npolys)
Definition: bch.h:44
bitvect< T, N > shiftdivmod(const Tm *m, size_t nm, const bitvect< T, N > &p, T init=0)
Definition: discrmath.h:138