IT++ Logo
reedsolomon.cpp
Go to the documentation of this file.
1 
29 #include <itpp/comm/reedsolomon.h>
30 #include <itpp/base/specmat.h>
31 #include <itpp/base/math/log_exp.h>
32 
33 namespace itpp
34 {
35 
36 //-------------------- Help Function ----------------------------
37 
40 {
41  int degree = f.get_true_degree();
42  int q = f.get_size();
43  int i;
44  GFX fprim(q, degree);
45  fprim.clear();
46  for (i = 0; i <= degree - 1; i += 2) {
47  fprim[i] = f[i + 1];
48  }
49  return fprim;
50 }
51 
52 //-------------------- Reed-Solomon ----------------------------
53 //A Reed-Solomon code is a q^m-ary BCH code of length n = pow(q,m)-1.
54 //k = pow(q,m)-1-t. This class works for q==2.
55 Reed_Solomon::Reed_Solomon(int in_m, int in_t, bool sys, int in_b):
56  m(in_m), t(in_t), b(in_b), systematic(sys)
57 {
58  n = pow2i(m) - 1;
59  k = pow2i(m) - 1 - 2 * t;
60  q = pow2i(m);
61  it_assert( (b >= 0) && (b < n), "Reed_Solomon::Reed_Solomon: narrow-sense parameter restricted to 0 <= b <= n.");
62  GFX x(q, (char *)"-1 0");
63  ivec alphapow(1);
64  g.set(q, (char *)"0");
65  for (int i = 1; i <= 2 * t; i++) {
66  alphapow(0) = b + i - 1;
67  g *= (x - GFX(q, alphapow));
68  }
69 }
70 
71 void Reed_Solomon::encode(const bvec &uncoded_bits, bvec &coded_bits)
72 {
73  int i, j, iterations = floor_i(static_cast<double>(uncoded_bits.length())
74  / (k * m));
75  GFX mx(q, k), cx(q, n);
76  GFX r(n + 1, n - k);
77  GFX uncoded_shifted(n + 1, n);
78  GF mpow;
79  bvec mbit(k * m), cbit(m);
80 
81  coded_bits.set_size(iterations * n * m, false);
82 
83  if (systematic)
84  for (i = 0; i < n - k; i++)
85  uncoded_shifted[i] = GF(n + 1, -1);
86 
87  for (i = 0; i < iterations; i++) {
88  //Fix the message polynom m(x).
89  for (j = 0; j < k; j++) {
90  mpow.set(q, uncoded_bits.mid((i * m * k) + (j * m), m));
91  mx[j] = mpow;
92  if (systematic) {
93  cx[j] = mx[j];
94  uncoded_shifted[j + n - k] = mx[j];
95  }
96  }
97  //Fix the outputbits cbit.
98  if (systematic) {
99  r = modgfx(uncoded_shifted, g);
100  for (j = k; j < n; j++) {
101  cx[j] = r[j - k];
102  }
103  }
104  else {
105  cx = g * mx;
106  }
107  for (j = 0; j < n; j++) {
108  cbit = cx[j].get_vectorspace();
109  coded_bits.replace_mid((i * n * m) + (j * m), cbit);
110  }
111  }
112 }
113 
114 bvec Reed_Solomon::encode(const bvec &uncoded_bits)
115 {
116  bvec coded_bits;
117  encode(uncoded_bits, coded_bits);
118  return coded_bits;
119 }
120 
121 bool Reed_Solomon::decode(const bvec &coded_bits, const ivec &erasure_positions, bvec &decoded_message, bvec &cw_isvalid)
122 {
123  bool decoderfailure, no_dec_failure;
124  int j, i, kk, l, L, foundzeros, iterations = floor_i(static_cast<double>(coded_bits.length()) / (n * m));
125  bvec mbit(m * k);
126  decoded_message.set_size(iterations * k * m, false);
127  cw_isvalid.set_length(iterations);
128 
129  GFX rx(q, n - 1), cx(q, n - 1), mx(q, k - 1), ex(q, n - 1), S(q, 2 * t), Xi(q, 2 * t), Gamma(q), Lambda(q),
130  Psiprime(q), OldLambda(q), T(q), Omega(q);
131  GFX dummy(q), One(q, (char*)"0"), Omegatemp(q);
132  GF delta(q), tempsum(q), rtemp(q), temp(q), Xk(q), Xkinv(q);
133  ivec errorpos;
134 
135  if ( erasure_positions.length() ) {
136  it_assert(max(erasure_positions) < iterations*n, "Reed_Solomon::decode: erasure position is invalid.");
137  }
138 
139  no_dec_failure = true;
140  for (i = 0; i < iterations; i++) {
141  decoderfailure = false;
142  //Fix the received polynomial r(x)
143  for (j = 0; j < n; j++) {
144  rtemp.set(q, coded_bits.mid(i * n * m + j * m, m));
145  rx[j] = rtemp;
146  }
147  // Fix the Erasure polynomial Gamma(x)
148  // and replace erased coordinates with zeros
149  rtemp.set(q, -1);
150  ivec alphapow = - ones_i(2);
151  Gamma = One;
152  for (j = 0; j < erasure_positions.length(); j++) {
153  rx[erasure_positions(j)] = rtemp;
154  alphapow(1) = erasure_positions(j);
155  Gamma *= (One - GFX(q, alphapow));
156  }
157  //Fix the syndrome polynomial S(x).
158  S.clear();
159  for (j = 1; j <= 2 * t; j++) {
160  S[j] = rx(GF(q, b + j - 1));
161  }
162  // calculate the modified syndrome polynomial Xi(x) = Gamma * (1+S) - 1
163  Xi = Gamma * (One + S) - One;
164  // Apply Berlekam-Massey algorithm
165  if (Xi.get_true_degree() >= 1) { //Errors in the received word
166  // Iterate to find Lambda(x), which hold all error locations
167  kk = 0;
168  Lambda = One;
169  L = 0;
170  T = GFX(q, (char*)"-1 0");
171  while (kk < 2 * t) {
172  kk = kk + 1;
173  tempsum = GF(q, -1);
174  for (l = 1; l <= L; l++) {
175  tempsum += Lambda[l] * Xi[kk - l];
176  }
177  delta = Xi[kk] - tempsum;
178  if (delta != GF(q, -1)) {
179  OldLambda = Lambda;
180  Lambda -= delta * T;
181  if (2 * L < kk) {
182  L = kk - L;
183  T = OldLambda / delta;
184  }
185  }
186  T = GFX(q, (char*)"-1 0") * T;
187  }
188  // Find the zeros to Lambda(x)
189  errorpos.set_size(Lambda.get_true_degree());
190  foundzeros = 0;
191  for (j = q - 2; j >= 0; j--) {
192  if (Lambda(GF(q, j)) == GF(q, -1)) {
193  errorpos(foundzeros) = (n - j) % n;
194  foundzeros += 1;
195  if (foundzeros >= Lambda.get_true_degree()) {
196  break;
197  }
198  }
199  }
200  if (foundzeros != Lambda.get_true_degree()) {
201  decoderfailure = true;
202  }
203  else { // Forney algorithm...
204  //Compute Omega(x) using the key equation for RS-decoding
205  Omega.set_degree(2 * t);
206  Omegatemp = Lambda * (One + Xi);
207  for (j = 0; j <= 2 * t; j++) {
208  Omega[j] = Omegatemp[j];
209  }
210  //Find the error/erasure magnitude polynomial by treating them the same
211  Psiprime = formal_derivate(Lambda*Gamma);
212  errorpos = concat(errorpos, erasure_positions);
213  ex.clear();
214  for (j = 0; j < errorpos.length(); j++) {
215  Xk = GF(q, errorpos(j));
216  Xkinv = GF(q, 0) / Xk;
217  // we calculate ex = - error polynomial, in order to avoid the
218  // subtraction when recunstructing the corrected codeword
219  ex[errorpos(j)] = (Xk * Omega(Xkinv)) / Psiprime(Xkinv);
220  if (b != 1) { // non-narrow-sense code needs corrected error magnitudes
221  int correction_exp = ( errorpos(j)*(1-b) ) % n;
222  ex[errorpos(j)] *= GF(q, correction_exp + ( (correction_exp < 0) ? n : 0 ));
223  }
224  }
225  //Reconstruct the corrected codeword.
226  // instead of subtracting the error/erasures, we calculated
227  // the negative error with 'ex' above
228  cx = rx + ex;
229  //Code word validation
230  S.clear();
231  for (j = 1; j <= 2 * t; j++) {
232  S[j] = cx(GF(q, b + j - 1));
233  }
234  if (S.get_true_degree() >= 1) {
235  decoderfailure = true;
236  }
237  }
238  }
239  else {
240  cx = rx;
241  decoderfailure = false;
242  }
243  //Find the message polynomial
244  mbit.clear();
245  if (decoderfailure == false) {
246  if (cx.get_true_degree() >= 1) { // A nonzero codeword was transmitted
247  if (systematic) {
248  for (j = 0; j < k; j++) {
249  mx[j] = cx[j];
250  }
251  }
252  else {
253  mx = divgfx(cx, g);
254  }
255  for (j = 0; j <= mx.get_true_degree(); j++) {
256  mbit.replace_mid(j * m, mx[j].get_vectorspace());
257  }
258  }
259  }
260  else { //Decoder failure.
261  // for a systematic code it is better to extract the undecoded message
262  // from the received code word, i.e. obtaining a bit error
263  // prob. p_b << 1/2, than setting all-zero (p_b = 1/2)
264  if (systematic) {
265  mbit = coded_bits.mid(i * n * m, k * m);
266  }
267  else {
268  mbit = zeros_b(k);
269  }
270  no_dec_failure = false;
271  }
272  decoded_message.replace_mid(i * m * k, mbit);
273  cw_isvalid(i) = (!decoderfailure);
274  }
275  return no_dec_failure;
276 }
277 
278 bool Reed_Solomon::decode(const bvec &coded_bits, bvec &decoded_message, bvec &cw_isvalid)
279 {
280  ivec erasures(0);
281  return decode(coded_bits, erasures, decoded_message, cw_isvalid);
282 }
283 
284 void Reed_Solomon::decode(const bvec &coded_bits, bvec &decoded_bits)
285 {
286  bvec cw_isvalid;
287  ivec erasures(0);
288  if (!decode(coded_bits, erasures, decoded_bits, cw_isvalid)) {
289  for (int i = 0; i < cw_isvalid.length(); i++) {
290  if (!cw_isvalid(i)) {
291  decoded_bits.replace_mid(i * k * m, zeros_b(k * m));
292  }
293  }
294  }
295 }
296 
297 bvec Reed_Solomon::decode(const bvec &coded_bits)
298 {
299  bvec decoded_bits;
300  decode(coded_bits, decoded_bits);
301  return decoded_bits;
302 }
303 
304 // --- Soft-decision decoding is not implemented ---
305 
306 void Reed_Solomon::decode(const vec &, bvec &)
307 {
308  it_error("Reed_Solomon::decode(): Soft-decision decoding not implemented");
309 }
310 
311 bvec Reed_Solomon::decode(const vec &)
312 {
313  it_error("Reed_Solomon::decode(): Soft-decision decoding not implemented");
314  return bvec();
315 }
316 
317 } // namespace itpp
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:23 for IT++ by Doxygen 1.8.2