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
Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4