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