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 Sun Dec 9 17:31:03 2007 for IT++ by Doxygen 1.5.4