IT++ Logo

reedsolomon.cpp

Go to the documentation of this file.
00001 
00030 #include <itpp/comm/reedsolomon.h>
00031 #include <itpp/base/math/log_exp.h>
00032 
00033 namespace itpp
00034 {
00035 
00036 //-------------------- Help Function ----------------------------
00037 
00039 GFX formal_derivate(const GFX &f)
00040 {
00041   int degree = f.get_true_degree();
00042   int q = f.get_size();
00043   int i;
00044   GFX fprim(q, degree);
00045   fprim.clear();
00046   for (i = 0; i <= degree - 1; i += 2) {
00047     fprim[i] = f[i+1];
00048   }
00049   return fprim;
00050 }
00051 
00052 //-------------------- Reed-Solomon ----------------------------
00053 //A Reed-Solomon code is a q^m-ary BCH code of length n = pow(q,m)-1.
00054 //k = pow(q,m)-1-t. This class works for q==2.
00055 Reed_Solomon::Reed_Solomon(int in_m, int in_t, bool sys):
00056     m(in_m), t(in_t), systematic(sys)
00057 {
00058   n = pow2i(m) - 1;
00059   k = pow2i(m) - 1 - 2 * t;
00060   q = pow2i(m);
00061   GFX x(q, (char *)"-1 0");
00062   ivec alphapow(1);
00063   g.set(q, (char *)"0");
00064   for (int i = 1; i <= 2*t; i++) {
00065     alphapow(0) = i;
00066     g *= (x - GFX(q, alphapow));
00067   }
00068 }
00069 
00070 void Reed_Solomon::encode(const bvec &uncoded_bits, bvec &coded_bits)
00071 {
00072   int i, j, itterations = floor_i(static_cast<double>(uncoded_bits.length())
00073                                   / (k * m));
00074   GFX mx(q, k), cx(q, n);
00075   GFX r(n + 1, n - k);
00076   GFX uncoded_shifted(n + 1, n);
00077   GF mpow;
00078   bvec mbit(k*m), cbit(m);
00079 
00080   coded_bits.set_size(itterations*n*m, false);
00081 
00082   if (systematic)
00083     for (i = 0; i < n - k; i++)
00084       uncoded_shifted[i] = GF(n + 1, -1);
00085 
00086   for (i = 0; i < itterations; i++) {
00087     //Fix the message polynom m(x).
00088     for (j = 0; j < k; j++) {
00089       mpow.set(q, uncoded_bits.mid((i*m*k) + (j*m), m));
00090       mx[j] = mpow;
00091       if (systematic) {
00092         cx[j] = mx[j];
00093         uncoded_shifted[j+n-k] = mx[j];
00094       }
00095     }
00096     //Fix the outputbits cbit.
00097     if (systematic) {
00098       r = modgfx(uncoded_shifted, g);
00099       for (j = k; j < n; j++) {
00100         cx[j] = r[j-k];
00101       }
00102     }
00103     else {
00104       cx = g * mx;
00105     }
00106     for (j = 0; j < n; j++) {
00107       cbit = cx[j].get_vectorspace();
00108       coded_bits.replace_mid((i*n*m) + (j*m), cbit);
00109     }
00110   }
00111 }
00112 
00113 bvec Reed_Solomon::encode(const bvec &uncoded_bits)
00114 {
00115   bvec coded_bits;
00116   encode(uncoded_bits, coded_bits);
00117   return coded_bits;
00118 }
00119 
00120 void Reed_Solomon::decode(const bvec &coded_bits, bvec &decoded_bits)
00121 {
00122   int j, i, kk, l, L, foundzeros, decoderfailure,
00123   itterations = floor_i(static_cast<double>(coded_bits.length()) / (n * m));
00124   bvec mbit(m*k);
00125   decoded_bits.set_size(itterations*k*m, false);
00126 
00127   GFX rx(q, n - 1), cx(q, n - 1), mx(q, k - 1), ex(q, n - 1), S(q, 2*t), Lambda(q),
00128   Lambdaprim(q), OldLambda(q), T(q), Ohmega(q);
00129   GFX dummy(q), One(q, (char*)"0"), Ohmegatemp(q);
00130   GF delta(q), tempsum(q), rtemp(q), temp(q), Xk(q), Xkinv(q);
00131   ivec errorpos;
00132 
00133   for (i = 0; i < itterations; i++) {
00134     decoderfailure = false;
00135     //Fix the received polynomial r(x)
00136     for (j = 0; j < n; j++) {
00137       rtemp.set(q, coded_bits.mid(i*n*m + j*m, m));
00138       rx[j] = rtemp;
00139     }
00140     //Fix the syndrome polynomial S(x).
00141     S.clear();
00142     for (j = 1; j <= 2*t; j++) {
00143       S[j] = rx(GF(q, j));
00144     }
00145     if (S.get_true_degree() == 0) {
00146       cx = rx;
00147       decoderfailure = false;
00148     }
00149     else {//Errors in the received word
00150       //Itterate to find Lambda(x).
00151       kk = 0;
00152       Lambda = GFX(q, (char*)"0");
00153       L = 0;
00154       T = GFX(q, (char*)"-1 0");
00155       while (kk < 2*t) {
00156         kk = kk + 1;
00157         tempsum = GF(q, -1);
00158         for (l = 1; l <= L; l++) {
00159           tempsum += Lambda[l] * S[kk-l];
00160         }
00161         delta = S[kk] - tempsum;
00162         if (delta != GF(q, -1)) {
00163           OldLambda = Lambda;
00164           Lambda -= delta * T;
00165           if (2*L < kk) {
00166             L = kk - L;
00167             T = OldLambda / delta;
00168           }
00169         }
00170         T = GFX(q, (char*)"-1 0") * T;
00171       }
00172       //Find the zeros to Lambda(x).
00173       errorpos.set_size(Lambda.get_true_degree(), false);
00174       errorpos.clear();
00175       foundzeros = 0;
00176       for (j = q - 2; j >= 0; j--) {
00177         temp = Lambda(GF(q, j));
00178         if (Lambda(GF(q, j)) == GF(q, -1)) {
00179           errorpos(foundzeros) = (n - j) % n;
00180           foundzeros += 1;
00181           if (foundzeros >= Lambda.get_true_degree()) {
00182             break;
00183           }
00184         }
00185       }
00186       if (foundzeros != Lambda.get_true_degree()) {
00187         decoderfailure = false;
00188       }
00189       else {
00190         //Compute Ohmega(x) using the key equation for RS-decoding
00191         Ohmega.set_degree(2*t);
00192         Ohmegatemp = Lambda * (One + S);
00193         for (j = 0; j <= 2*t; j++) {
00194           Ohmega[j] = Ohmegatemp[j];
00195         }
00196         Lambdaprim = formal_derivate(Lambda);
00197         //Find the error polynomial
00198         ex.clear();
00199         for (j = 0; j < foundzeros; j++) {
00200           Xk = GF(q, errorpos(j));
00201           Xkinv = GF(q, 0) / Xk;
00202           ex[errorpos(j)] = (Xk * Ohmega(Xkinv)) / Lambdaprim(Xkinv);
00203         }
00204         //Reconstruct the corrected codeword.
00205         cx = rx + ex;
00206         //Code word validation
00207         S.clear();
00208         for (j = 1; j <= 2*t; j++) {
00209           S[j] = rx(GF(q, j));
00210         }
00211         if (S.get_true_degree() >= 1) {
00212           decoderfailure = false;
00213         }
00214       }
00215     }
00216     //Find the message polynomial
00217     mbit.clear();
00218     if (decoderfailure == false) {
00219       if (cx.get_true_degree() >= 1) {// A nonzero codeword was transmitted
00220         if (systematic) {
00221           for (j = 0; j < k; j++) {
00222             mx[j] = cx[j];
00223           }
00224         }
00225         else {
00226           mx = divgfx(cx, g);
00227         }
00228         for (j = 0; j <= mx.get_true_degree(); j++) {
00229           mbit.replace_mid(j*m, mx[j].get_vectorspace());
00230         }
00231       }
00232     }
00233     decoded_bits.replace_mid(i*m*k, mbit);
00234   }
00235 }
00236 
00237 bvec Reed_Solomon::decode(const bvec &coded_bits)
00238 {
00239   bvec decoded_bits;
00240   decode(coded_bits, decoded_bits);
00241   return decoded_bits;
00242 }
00243 
00244 // --- Soft-decision decoding is not implemented ---
00245 
00246 void Reed_Solomon::decode(const vec &, bvec &)
00247 {
00248   it_error("Reed_Solomon::decode(): Soft-decision decoding not implemented");
00249 }
00250 
00251 bvec Reed_Solomon::decode(const vec &)
00252 {
00253   it_error("Reed_Solomon::decode(): Soft-decision decoding not implemented");
00254   return bvec();
00255 }
00256 
00257 } // namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jan 20 23:03:05 2010 for IT++ by Doxygen 1.6.2