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

Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4