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

Generated on Thu Apr 24 13:39:00 2008 for IT++ by Doxygen 1.5.5