IT++ 4.3.1
reedsolomon.cpp
Go to the documentation of this file.
1
28
30#include <itpp/base/specmat.h>
32
33namespace 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.
55Reed_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
71void 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
114bvec Reed_Solomon::encode(const bvec &uncoded_bits)
115{
116 bvec coded_bits;
117 encode(uncoded_bits, coded_bits);
118 return coded_bits;
119}
120
121bool 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
278bool 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
284void 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
297bvec 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
306void Reed_Solomon::decode(const vec &, bvec &)
307{
308 it_error("Reed_Solomon::decode(): Soft-decision decoding not implemented");
309}
310
311bvec 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
Polynomials over GF(q)[x], where q=2^m, m=1,...,16.
Definition galois.h:176
int get_true_degree() const
Return true degree of GF(q)[x].
Definition galois.h:466
void set_degree(int indegree, bool copy=false)
Resize the polynomial to the given indegree. If copy is set to true, the old polynomial's coefficient...
Definition galois.h:459
int get_size() const
Return q.
Definition galois.h:449
void clear()
Set all coefficients to zero.
Definition galois.h:477
Galois Field GF(q).
Definition galois.h:75
void set(int qvalue, int inexp)
GF(q) equals alpha ^ inexp.
Definition galois.h:92
virtual void encode(const bvec &uncoded_bits, bvec &coded_bits)
Encoder function.
virtual bool decode(const bvec &coded_bits, const ivec &erasure_positions, bvec &decoded_message, bvec &cw_isvalid)
Decode the RS code bits. Return false if there has been any decoding failure.
Reed_Solomon(int in_m, int in_t, bool sys=false, int in_b=1)
GFX g
The generator polynomial of the RS code.
const bool systematic
Whether or not the code is systematic.
#define it_error(s)
Abort unconditionally.
Definition itassert.h:126
#define it_assert(t, s)
Abort if t is not true.
Definition itassert.h:94
int pow2i(int x)
Calculate two to the power of x (2^x); x is integer.
Definition log_exp.h:53
T max(const Vec< T > &v)
Maximum value of vector.
Definition min_max.h:45
ITPP_EXPORT ivec ones_i(int size)
A Int vector of ones.
ITPP_EXPORT bvec zeros_b(int size)
A Binary vector of zeros.
Logarithmic and exponenential functions - header file.
itpp namespace
Definition itmex.h:37
const Array< T > concat(const Array< T > &a, const T &e)
Append element e to the end of the Array a.
Definition array.h:486
GFX formal_derivate(const GFX &f)
Local help function.
GFX divgfx(const GFX &c, const GFX &g)
Division of two GFX (local help function)
Definition galois.cpp:157
int floor_i(double x)
The nearest smaller integer.
Definition converters.h:350
GFX modgfx(const GFX &a, const GFX &b)
Modulo function of two GFX (local help function)
Definition galois.cpp:183
Definitions of a Reed-Solomon codec class.
Definitions of special vectors and matrices.