00001 00033 #include <itpp/comm/reedsolomon.h> 00034 00035 00036 namespace itpp { 00037 00038 //-------------------- Help Function ---------------------------- 00039 00041 GFX formal_derivate(const GFX &f) { 00042 int degree = f.get_true_degree(); 00043 int q = f.get_size(); 00044 int i; 00045 GFX fprim(q,degree); 00046 fprim.clear(); 00047 for (i=0; i<=degree-1; i+=2) { 00048 fprim[i] = f[i+1]; 00049 } 00050 return fprim; 00051 } 00052 00053 //-------------------- Reed-Solomon ---------------------------- 00054 //A Reed-Solomon code is a q^m-ary BCH code of length n = pow(q,m)-1. 00055 //k = pow(q,m)-1-t. This class works for q==2. 00056 Reed_Solomon::Reed_Solomon(int in_m, int in_t) { 00057 m = in_m; 00058 t = in_t; 00059 n = pow2i(m) - 1; 00060 k = pow2i(m) - 1 - 2*t; 00061 q = pow2i(m); 00062 int i; 00063 GFX x(q,(char *)"-1 0"); 00064 ivec alphapow(1); 00065 g.set(q,(char *)"0"); 00066 for(i=1; i<=2*t; i++) { 00067 alphapow(0) = i; 00068 g *= (x-GFX(q,alphapow)); 00069 } 00070 } 00071 00072 void Reed_Solomon::encode(const bvec &uncoded_bits, bvec &coded_bits) 00073 { 00074 int i, j, itterations = (int)floor( (double)uncoded_bits.length() / (k*m) ); 00075 GFX mx(q,k), cx(q,n); 00076 GF mpow; 00077 bvec mbit(k*m), cbit(m); 00078 00079 coded_bits.set_size(itterations*n*m, false); 00080 00081 for (i=0; i<itterations; i++) { 00082 //Fix the message polynom m(x). 00083 for (j=0; j<k; j++) { 00084 mpow.set(q,uncoded_bits.mid((i*m*k)+(j*m),m)); 00085 mx[j] = mpow; 00086 } 00087 //Fix the outputbits cbit. 00088 cx = g*mx; 00089 for (j=0; j<n; j++) { 00090 cbit = cx[j].get_vectorspace(); 00091 coded_bits.replace_mid((i*n*m)+(j*m), cbit); 00092 } 00093 } 00094 } 00095 00096 bvec Reed_Solomon::encode(const bvec &uncoded_bits) 00097 { 00098 bvec coded_bits; 00099 encode(uncoded_bits, coded_bits); 00100 return coded_bits; 00101 } 00102 00103 void Reed_Solomon::decode(const bvec &coded_bits, bvec &decoded_bits) 00104 { 00105 int j, i, kk, l, L, foundzeros, itterations = (int)floor( (double)coded_bits.length() / (n*m) ), decoderfailure; 00106 bvec mbit(m*k); 00107 decoded_bits.set_size(itterations*k*m, false); 00108 00109 GFX rx(q,n-1), cx(q,n-1), mx(q,k-1), ex(q,n-1), S(q,2*t), Lambda(q), Lambdaprim(q), OldLambda(q), T(q), Ohmega(q); 00110 GFX dummy(q), One(q,(char*)"0"), Ohmegatemp(q); 00111 GF delta(q), tempsum(q), rtemp(q), temp(q), Xk(q), Xkinv(q); 00112 ivec errorpos; 00113 00114 for (i=0; i<itterations; i++) { 00115 decoderfailure = false; 00116 //Fix the received polynomial r(x) 00117 for (j=0; j<n; j++) { rtemp.set(q,coded_bits.mid(i*n*m + j*m, m)); rx[j] = rtemp; } 00118 //Fix the syndrome polynomial S(x). 00119 S.clear(); 00120 for (j=1; j<=2*t; j++) { S[j] = rx(GF(q,j)); } 00121 if (S.get_true_degree() == 0) {cx = rx; decoderfailure = false; } 00122 else {//Errors in the received word 00123 //Itterate to find Lambda(x). 00124 kk = 0; Lambda = GFX(q,(char*)"0"); L = 0; T = GFX(q,(char*)"-1 0"); 00125 while (kk<2*t) { 00126 kk = kk + 1; 00127 tempsum = GF(q,-1); 00128 for (l=1; l<=L; l++) { tempsum += Lambda[l] * S[kk-l]; } 00129 delta = S[kk] - tempsum; 00130 if (delta != GF(q,-1)) { 00131 OldLambda = Lambda; 00132 Lambda -= delta*T; 00133 if (2*L<kk) { L = kk - L; T = OldLambda/delta;} 00134 } 00135 T = GFX(q,(char*)"-1 0") * T; 00136 } 00137 //Find the zeros to Lambda(x). 00138 errorpos.set_size(Lambda.get_true_degree(), false); 00139 errorpos.clear(); 00140 foundzeros = 0; 00141 for (j=q-2; j>=0; j--) { 00142 temp = Lambda( GF(q,j) ); 00143 if (Lambda( GF(q,j) ) == GF(q,-1) ) { 00144 errorpos( foundzeros ) = (n-j) % n; 00145 foundzeros +=1; 00146 if (foundzeros >= Lambda.get_true_degree()) { break; } 00147 } 00148 } 00149 if (foundzeros != Lambda.get_true_degree()) { decoderfailure = false; } 00150 else { 00151 //Compute Ohmega(x) using the key equation for RS-decoding 00152 Ohmega.set_degree(2*t); 00153 Ohmegatemp = Lambda * (One +S); 00154 for (j=0; j<=2*t; j++) { Ohmega[j] = Ohmegatemp[j]; } 00155 Lambdaprim = formal_derivate(Lambda); 00156 //Find the error polynomial 00157 ex.clear(); 00158 for (j=0; j<foundzeros; j++) { 00159 Xk = GF(q,errorpos(j)); Xkinv = GF(q,0) / Xk; 00160 ex[errorpos(j)] = (Xk * Ohmega(Xkinv)) / Lambdaprim(Xkinv); 00161 } 00162 //Reconstruct the corrected codeword. 00163 cx = rx + ex; 00164 //Code word validation 00165 S.clear(); for (j=1; j<=2*t; j++) { S[j] = rx(GF(q,j)); } 00166 if (S.get_true_degree() >= 1) { decoderfailure=false; } 00167 } 00168 } 00169 //Find the message polynomial 00170 mbit.clear(); 00171 if (decoderfailure == false) { 00172 if (cx.get_true_degree() >= 1) {// A nonzero codeword was transmitted 00173 mx = divgfx(cx,g); 00174 for (j=0; j<=mx.get_true_degree(); j++) { mbit.replace_mid(j*m,mx[j].get_vectorspace()); } 00175 } 00176 } 00177 decoded_bits.replace_mid(i*m*k,mbit); 00178 } 00179 } 00180 00181 bvec Reed_Solomon::decode(const bvec &coded_bits) 00182 { 00183 bvec decoded_bits; 00184 decode(coded_bits, decoded_bits); 00185 return decoded_bits; 00186 } 00187 00188 // --------------- Soft-decision decoding is not implemented -------------------------------- 00189 void Reed_Solomon::decode(const vec &received_signal, bvec &output) 00190 { 00191 it_error("Reed_Solomon::decode(vec, bvec); soft-decision decoding is not implemented"); 00192 } 00193 00194 bvec Reed_Solomon::decode(const vec &received_signal) 00195 { 00196 it_error("Reed_Solomon::decode(vec, bvec); soft-decision decoding is not implemented"); 00197 return bvec(); 00198 } 00199 00200 00201 } // namespace itpp
Generated on Wed Apr 18 11:23:33 2007 for IT++ by Doxygen 1.5.2