00001 00029 #include <itpp/comm/rec_syst_conv_code.h> 00030 00031 00032 namespace itpp 00033 { 00034 00036 double(*com_log)(double, double) = NULL; 00037 00039 // This wrapper is because "com_log = std::max;" below caused an error 00040 inline double max(double x, double y) { return std::max(x, y); } 00042 00043 // ----------------- Protected functions ----------------------------- 00044 00045 int Rec_Syst_Conv_Code::calc_state_transition(const int instate, const int input, ivec &parity) 00046 { 00047 int i, j, in = 0, temp = (gen_pol_rev(0) & (instate << 1)), parity_temp, parity_bit; 00048 00049 for (i = 0; i < K; i++) { 00050 in = (temp & 1) ^ in; 00051 temp = temp >> 1; 00052 } 00053 in = in ^ input; 00054 00055 parity.set_size(n - 1, false); 00056 for (j = 0; j < (n - 1); j++) { 00057 parity_temp = ((instate << 1) + in) & gen_pol_rev(j + 1); 00058 parity_bit = 0; 00059 for (i = 0; i < K; i++) { 00060 parity_bit = (parity_temp & 1) ^ parity_bit; 00061 parity_temp = parity_temp >> 1; 00062 } 00063 parity(j) = parity_bit; 00064 } 00065 return in + ((instate << 1) & ((1 << m) - 1)); 00066 } 00067 00068 // --------------- Public functions ------------------------- 00069 void Rec_Syst_Conv_Code::set_generator_polynomials(const ivec &gen, int constraint_length) 00070 { 00071 int j; 00072 gen_pol = gen; 00073 n = gen.size(); 00074 K = constraint_length; 00075 m = K - 1; 00076 rate = 1.0 / n; 00077 00078 gen_pol_rev.set_size(n, false); 00079 for (int i = 0; i < n; i++) { 00080 gen_pol_rev(i) = reverse_int(K, gen_pol(i)); 00081 } 00082 00083 Nstates = (1 << m); 00084 state_trans.set_size(Nstates, 2, false); 00085 rev_state_trans.set_size(Nstates, 2, false); 00086 output_parity.set_size(Nstates, 2*(n - 1), false); 00087 rev_output_parity.set_size(Nstates, 2*(n - 1), false); 00088 int s0, s1, s_prim; 00089 ivec p0, p1; 00090 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00091 s0 = calc_state_transition(s_prim, 0, p0); 00092 state_trans(s_prim, 0) = s0; 00093 rev_state_trans(s0, 0) = s_prim; 00094 for (j = 0; j < (n - 1); j++) { 00095 output_parity(s_prim, 2*j + 0) = p0(j); 00096 rev_output_parity(s0, 2*j + 0) = p0(j); 00097 } 00098 00099 s1 = calc_state_transition(s_prim, 1, p1); 00100 state_trans(s_prim, 1) = s1; 00101 rev_state_trans(s1, 1) = s_prim; 00102 for (j = 0; j < (n - 1); j++) { 00103 output_parity(s_prim, 2*j + 1) = p1(j); 00104 rev_output_parity(s1, 2*j + 1) = p1(j); 00105 } 00106 } 00107 00108 ln2 = std::log(2.0); 00109 00110 //The default value of Lc is 1: 00111 Lc = 1.0; 00112 } 00113 00114 void Rec_Syst_Conv_Code::set_awgn_channel_parameters(double Ec, double N0) 00115 { 00116 Lc = 4.0 * std::sqrt(Ec) / N0; 00117 } 00118 00119 void Rec_Syst_Conv_Code::set_scaling_factor(double in_Lc) 00120 { 00121 Lc = in_Lc; 00122 } 00123 00124 void Rec_Syst_Conv_Code::encode_tail(const bvec &input, bvec &tail, bmat &parity_bits) 00125 { 00126 int i, j, length = input.size(), target_state; 00127 parity_bits.set_size(length + m, n - 1, false); 00128 tail.set_size(m, false); 00129 00130 encoder_state = 0; 00131 for (i = 0; i < length; i++) { 00132 for (j = 0; j < (n - 1); j++) { 00133 parity_bits(i, j) = output_parity(encoder_state, 2 * j + int(input(i))); 00134 } 00135 encoder_state = state_trans(encoder_state, int(input(i))); 00136 } 00137 00138 // add tail of m=K-1 zeros 00139 for (i = 0; i < m; i++) { 00140 target_state = (encoder_state << 1) & ((1 << m) - 1); 00141 if (state_trans(encoder_state, 0) == target_state) { tail(i) = bin(0); } 00142 else { tail(i) = bin(1); } 00143 for (j = 0; j < (n - 1); j++) { 00144 parity_bits(length + i, j) = output_parity(encoder_state, 2 * j + int(tail(i))); 00145 } 00146 encoder_state = target_state; 00147 } 00148 terminated = true; 00149 } 00150 00151 void Rec_Syst_Conv_Code::encode(const bvec &input, bmat &parity_bits) 00152 { 00153 int i, j, length = input.size(); 00154 parity_bits.set_size(length, n - 1, false); 00155 00156 encoder_state = 0; 00157 for (i = 0; i < length; i++) { 00158 for (j = 0; j < (n - 1); j++) { 00159 parity_bits(i, j) = output_parity(encoder_state, 2 * j + int(input(i))); 00160 } 00161 encoder_state = state_trans(encoder_state, int(input(i))); 00162 } 00163 terminated = false; 00164 } 00165 00166 void Rec_Syst_Conv_Code::map_decode(const vec &rec_systematic, const mat &rec_parity, const vec &extrinsic_input, 00167 vec &extrinsic_output, bool in_terminated) 00168 { 00169 double gamma_k_e, nom, den, temp0, temp1, exp_temp0, exp_temp1; 00170 int j, s0, s1, k, kk, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length(); 00171 ivec p0, p1; 00172 00173 alpha.set_size(Nstates, block_length + 1, false); 00174 beta.set_size(Nstates, block_length + 1, false); 00175 gamma.set_size(2*Nstates, block_length + 1, false); 00176 denom.set_size(block_length + 1, false); 00177 denom.clear(); 00178 extrinsic_output.set_size(block_length, false); 00179 00180 if (in_terminated) { terminated = true; } 00181 00182 //Calculate gamma 00183 for (k = 1; k <= block_length; k++) { 00184 kk = k - 1; 00185 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00186 exp_temp0 = 0.0; 00187 exp_temp1 = 0.0; 00188 for (j = 0; j < (n - 1); j++) { 00189 exp_temp0 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 0)); /* Is this OK? */ 00190 exp_temp1 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 1)); 00191 } 00192 // gamma(2*s_prim+0,k) = std::exp( 0.5*(extrinsic_input(kk) + Lc*rec_systematic(kk))) * std::exp( exp_temp0 ); 00193 // gamma(2*s_prim+1,k) = std::exp(-0.5*(extrinsic_input(kk) + Lc*rec_systematic(kk))) * std::exp( exp_temp1 ); 00194 /* == Changed to trunc_exp() to address bug 1088420 which 00195 described a numerical instability problem in map_decode() 00196 at high SNR. This should be regarded as a temporary fix and 00197 it is not necessarily a waterproof one: multiplication of 00198 probabilities still can result in values out of 00199 range. (Range checking for the multiplication operator was 00200 not implemented as it was felt that this would sacrifice 00201 too much runtime efficiency. Some margin was added to the 00202 numerical hardlimits below to reflect this. The hardlimit 00203 values below were taken as the minimum range that a 00204 "double" should support reduced by a few orders of 00205 magnitude to make sure multiplication of several values 00206 does not exceed the limits.) 00207 00208 It is suggested to use the QLLR based log-domain() decoders 00209 instead of map_decode() as they are much faster and more 00210 numerically stable. 00211 00212 EGL 8/06. == */ 00213 gamma(2*s_prim + 0, k) = trunc_exp(0.5 * (extrinsic_input(kk) + Lc * rec_systematic(kk)) + exp_temp0); 00214 gamma(2*s_prim + 1, k) = trunc_exp(-0.5 * (extrinsic_input(kk) + Lc * rec_systematic(kk)) + exp_temp1); 00215 } 00216 } 00217 00218 //Initiate alpha 00219 alpha.set_col(0, zeros(Nstates)); 00220 alpha(0, 0) = 1.0; 00221 00222 //Calculate alpha and denom going forward through the trellis 00223 for (k = 1; k <= block_length; k++) { 00224 for (s = 0; s < Nstates; s++) { 00225 s_prim0 = rev_state_trans(s, 0); 00226 s_prim1 = rev_state_trans(s, 1); 00227 temp0 = alpha(s_prim0, k - 1) * gamma(2 * s_prim0 + 0, k); 00228 temp1 = alpha(s_prim1, k - 1) * gamma(2 * s_prim1 + 1, k); 00229 alpha(s, k) = temp0 + temp1; 00230 denom(k) += temp0 + temp1; 00231 } 00232 alpha.set_col(k, alpha.get_col(k) / denom(k)); 00233 } 00234 00235 //Initiate beta 00236 if (terminated) { 00237 beta.set_col(block_length, zeros(Nstates)); 00238 beta(0, block_length) = 1.0; 00239 } 00240 else { 00241 beta.set_col(block_length, alpha.get_col(block_length)); 00242 } 00243 00244 //Calculate beta going backward in the trellis 00245 for (k = block_length; k >= 2; k--) { 00246 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00247 s0 = state_trans(s_prim, 0); 00248 s1 = state_trans(s_prim, 1); 00249 beta(s_prim, k - 1) = (beta(s0, k) * gamma(2 * s_prim + 0, k)) + (beta(s1, k) * gamma(2 * s_prim + 1, k)); 00250 } 00251 beta.set_col(k - 1, beta.get_col(k - 1) / denom(k)); 00252 } 00253 00254 //Calculate extrinsic output for each bit 00255 for (k = 1; k <= block_length; k++) { 00256 kk = k - 1; 00257 nom = 0; 00258 den = 0; 00259 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00260 s0 = state_trans(s_prim, 0); 00261 s1 = state_trans(s_prim, 1); 00262 exp_temp0 = 0.0; 00263 exp_temp1 = 0.0; 00264 for (j = 0; j < (n - 1); j++) { 00265 exp_temp0 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 0)); 00266 exp_temp1 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 1)); 00267 } 00268 // gamma_k_e = std::exp( exp_temp0 ); 00269 gamma_k_e = trunc_exp(exp_temp0); 00270 nom += alpha(s_prim, k - 1) * gamma_k_e * beta(s0, k); 00271 00272 // gamma_k_e = std::exp( exp_temp1 ); 00273 gamma_k_e = trunc_exp(exp_temp1); 00274 den += alpha(s_prim, k - 1) * gamma_k_e * beta(s1, k); 00275 } 00276 // extrinsic_output(kk) = std::log(nom/den); 00277 extrinsic_output(kk) = trunc_log(nom / den); 00278 } 00279 00280 } 00281 00282 void Rec_Syst_Conv_Code::log_decode(const vec &rec_systematic, const mat &rec_parity, 00283 const vec &extrinsic_input, vec &extrinsic_output, bool in_terminated, std::string metric) 00284 { 00285 if (metric == "TABLE") { 00286 /* Use the QLLR decoder. This can probably be done more 00287 efficiently since it converts floating point vectors to QLLR. 00288 However we have to live with this for the time being. */ 00289 QLLRvec rec_systematic_q = llrcalc.to_qllr(rec_systematic); 00290 QLLRmat rec_parity_q = llrcalc.to_qllr(rec_parity); 00291 QLLRvec extrinsic_input_q = llrcalc.to_qllr(extrinsic_input); 00292 QLLRvec extrinsic_output_q(length(extrinsic_output)); 00293 log_decode(rec_systematic_q, rec_parity_q, extrinsic_input_q, 00294 extrinsic_output_q, in_terminated); 00295 extrinsic_output = llrcalc.to_double(extrinsic_output_q); 00296 return; 00297 } 00298 00299 double nom, den, exp_temp0, exp_temp1, rp, temp0, temp1; 00300 int i, j, s0, s1, k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length(); 00301 ivec p0, p1; 00302 00303 //Set the internal metric: 00304 if (metric == "LOGMAX") { com_log = max; } 00305 else if (metric == "LOGMAP") { com_log = log_add; } 00306 else { 00307 it_error("Rec_Syst_Conv_Code::log_decode: Illegal metric parameter"); 00308 } 00309 00310 alpha.set_size(Nstates, block_length + 1, false); 00311 beta.set_size(Nstates, block_length + 1, false); 00312 gamma.set_size(2*Nstates, block_length + 1, false); 00313 extrinsic_output.set_size(block_length, false); 00314 denom.set_size(block_length + 1, false); 00315 for (k = 0; k <= block_length; k++) { denom(k) = -infinity; } 00316 00317 if (in_terminated) { terminated = true; } 00318 00319 //Check that Lc = 1.0 00320 it_assert(Lc == 1.0, 00321 "Rec_Syst_Conv_Code::log_decode: This function assumes that Lc = 1.0. Please use proper scaling of the input data"); 00322 00323 //Calculate gamma 00324 for (k = 1; k <= block_length; k++) { 00325 kk = k - 1; 00326 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00327 exp_temp0 = 0.0; 00328 exp_temp1 = 0.0; 00329 for (j = 0; j < (n - 1); j++) { 00330 rp = rec_parity(kk, j); 00331 if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; } 00332 else { exp_temp0 -= rp; } 00333 if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; } 00334 else { exp_temp1 -= rp; } 00335 } 00336 gamma(2*s_prim + 0, k) = 0.5 * ((extrinsic_input(kk) + rec_systematic(kk)) + exp_temp0); 00337 gamma(2*s_prim + 1, k) = -0.5 * ((extrinsic_input(kk) + rec_systematic(kk)) - exp_temp1); 00338 } 00339 } 00340 00341 //Initiate alpha 00342 for (j = 1; j < Nstates; j++) { alpha(j, 0) = -infinity; } 00343 alpha(0, 0) = 0.0; 00344 00345 //Calculate alpha, going forward through the trellis 00346 for (k = 1; k <= block_length; k++) { 00347 for (s = 0; s < Nstates; s++) { 00348 s_prim0 = rev_state_trans(s, 0); 00349 s_prim1 = rev_state_trans(s, 1); 00350 temp0 = alpha(s_prim0, k - 1) + gamma(2 * s_prim0 + 0, k); 00351 temp1 = alpha(s_prim1, k - 1) + gamma(2 * s_prim1 + 1, k); 00352 alpha(s, k) = com_log(temp0, temp1); 00353 denom(k) = com_log(alpha(s, k), denom(k)); 00354 } 00355 //Normalization of alpha 00356 for (l = 0; l < Nstates; l++) { alpha(l, k) -= denom(k); } 00357 } 00358 00359 //Initiate beta 00360 if (terminated) { 00361 for (i = 1; i < Nstates; i++) { beta(i, block_length) = -infinity; } 00362 beta(0, block_length) = 0.0; 00363 } 00364 else { 00365 beta.set_col(block_length, alpha.get_col(block_length)); 00366 } 00367 00368 //Calculate beta going backward in the trellis 00369 for (k = block_length; k >= 1; k--) { 00370 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00371 s0 = state_trans(s_prim, 0); 00372 s1 = state_trans(s_prim, 1); 00373 beta(s_prim, k - 1) = com_log(beta(s0, k) + gamma(2 * s_prim + 0, k) , beta(s1, k) + gamma(2 * s_prim + 1, k)); 00374 } 00375 //Normalization of beta 00376 for (l = 0; l < Nstates; l++) { beta(l, k - 1) -= denom(k); } 00377 } 00378 00379 //Calculate extrinsic output for each bit 00380 for (k = 1; k <= block_length; k++) { 00381 kk = k - 1; 00382 nom = -infinity; 00383 den = -infinity; 00384 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00385 s0 = state_trans(s_prim, 0); 00386 s1 = state_trans(s_prim, 1); 00387 exp_temp0 = 0.0; 00388 exp_temp1 = 0.0; 00389 for (j = 0; j < (n - 1); j++) { 00390 rp = rec_parity(kk, j); 00391 if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; } 00392 else { exp_temp0 -= rp; } 00393 if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; } 00394 else { exp_temp1 -= rp; } 00395 } 00396 nom = com_log(nom, alpha(s_prim, kk) + 0.5 * exp_temp0 + beta(s0, k)); 00397 den = com_log(den, alpha(s_prim, kk) + 0.5 * exp_temp1 + beta(s1, k)); 00398 } 00399 extrinsic_output(kk) = nom - den; 00400 } 00401 00402 } 00403 00404 void Rec_Syst_Conv_Code::log_decode_n2(const vec &rec_systematic, const vec &rec_parity, 00405 const vec &extrinsic_input, vec &extrinsic_output, bool in_terminated, std::string metric) 00406 { 00407 if (metric == "TABLE") { // use the QLLR decoder; also see comment under log_decode() 00408 QLLRvec rec_systematic_q = llrcalc.to_qllr(rec_systematic); 00409 QLLRvec rec_parity_q = llrcalc.to_qllr(rec_parity); 00410 QLLRvec extrinsic_input_q = llrcalc.to_qllr(extrinsic_input); 00411 QLLRvec extrinsic_output_q(length(extrinsic_output)); 00412 log_decode_n2(rec_systematic_q, rec_parity_q, extrinsic_input_q, 00413 extrinsic_output_q, in_terminated); 00414 extrinsic_output = llrcalc.to_double(extrinsic_output_q); 00415 return; 00416 } 00417 00418 // const double INF = 10e300; // replaced by DEFINE to be file-wide in scope 00419 double nom, den, exp_temp0, exp_temp1, rp; 00420 int k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length(); 00421 int ext_info_length = extrinsic_input.length(); 00422 ivec p0, p1; 00423 double ex, norm; 00424 00425 //Set the internal metric: 00426 if (metric == "LOGMAX") { com_log = max; } 00427 else if (metric == "LOGMAP") { com_log = log_add; } 00428 else { 00429 it_error("Rec_Syst_Conv_Code::log_decode_n2: Illegal metric parameter"); 00430 } 00431 00432 alpha.set_size(Nstates, block_length + 1, false); 00433 beta.set_size(Nstates, block_length + 1, false); 00434 gamma.set_size(2*Nstates, block_length + 1, false); 00435 extrinsic_output.set_size(ext_info_length, false); 00436 //denom.set_size(block_length+1,false); for (k=0; k<=block_length; k++) { denom(k) = -infinity; } 00437 00438 if (in_terminated) { terminated = true; } 00439 00440 //Check that Lc = 1.0 00441 it_assert(Lc == 1.0, 00442 "Rec_Syst_Conv_Code::log_decode_n2: This function assumes that Lc = 1.0. Please use proper scaling of the input data"); 00443 00444 //Initiate alpha 00445 for (s = 1; s < Nstates; s++) { alpha(s, 0) = -infinity; } 00446 alpha(0, 0) = 0.0; 00447 00448 //Calculate alpha and gamma going forward through the trellis 00449 for (k = 1; k <= block_length; k++) { 00450 kk = k - 1; 00451 if (kk < ext_info_length) { 00452 ex = 0.5 * (extrinsic_input(kk) + rec_systematic(kk)); 00453 } 00454 else { 00455 ex = 0.5 * rec_systematic(kk); 00456 } 00457 rp = 0.5 * rec_parity(kk); 00458 for (s = 0; s < Nstates; s++) { 00459 s_prim0 = rev_state_trans(s, 0); 00460 s_prim1 = rev_state_trans(s, 1); 00461 if (output_parity(s_prim0 , 0)) { exp_temp0 = -rp; } 00462 else { exp_temp0 = rp; } 00463 if (output_parity(s_prim1 , 1)) { exp_temp1 = -rp; } 00464 else { exp_temp1 = rp; } 00465 gamma(2*s_prim0 , k) = ex + exp_temp0; 00466 gamma(2*s_prim1 + 1, k) = -ex + exp_temp1; 00467 alpha(s, k) = com_log(alpha(s_prim0, kk) + gamma(2 * s_prim0 , k), 00468 alpha(s_prim1, kk) + gamma(2 * s_prim1 + 1, k)); 00469 //denom(k) = com_log( alpha(s,k), denom(k) ); 00470 } 00471 norm = alpha(0, k); //norm = denom(k); 00472 for (l = 0; l < Nstates; l++) { alpha(l, k) -= norm; } 00473 } 00474 00475 //Initiate beta 00476 if (terminated) { 00477 for (s = 1; s < Nstates; s++) { beta(s, block_length) = -infinity; } 00478 beta(0, block_length) = 0.0; 00479 } 00480 else { 00481 beta.set_col(block_length, alpha.get_col(block_length)); 00482 } 00483 00484 //Calculate beta going backward in the trellis 00485 for (k = block_length; k >= 1; k--) { 00486 kk = k - 1; 00487 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00488 beta(s_prim, kk) = com_log(beta(state_trans(s_prim, 0), k) + gamma(2 * s_prim, k), 00489 beta(state_trans(s_prim, 1), k) + gamma(2 * s_prim + 1, k)); 00490 } 00491 norm = beta(0, k); //norm = denom(k); 00492 for (l = 0; l < Nstates; l++) { beta(l, k) -= norm; } 00493 } 00494 00495 //Calculate extrinsic output for each bit 00496 for (k = 1; k <= block_length; k++) { 00497 kk = k - 1; 00498 if (kk < ext_info_length) { 00499 nom = -infinity; 00500 den = -infinity; 00501 rp = 0.5 * rec_parity(kk); 00502 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00503 if (output_parity(s_prim , 0)) { exp_temp0 = -rp; } 00504 else { exp_temp0 = rp; } 00505 if (output_parity(s_prim , 1)) { exp_temp1 = -rp; } 00506 else { exp_temp1 = rp; } 00507 nom = com_log(nom, alpha(s_prim, kk) + exp_temp0 + beta(state_trans(s_prim, 0), k)); 00508 den = com_log(den, alpha(s_prim, kk) + exp_temp1 + beta(state_trans(s_prim, 1), k)); 00509 } 00510 extrinsic_output(kk) = nom - den; 00511 } 00512 } 00513 } 00514 00515 // === Below new decoder functions by EGL, using QLLR arithmetics =========== 00516 00517 void Rec_Syst_Conv_Code::log_decode(const QLLRvec &rec_systematic, const QLLRmat &rec_parity, 00518 const QLLRvec &extrinsic_input, 00519 QLLRvec &extrinsic_output, bool in_terminated) 00520 { 00521 00522 int nom, den, exp_temp0, exp_temp1, rp, temp0, temp1; 00523 int i, j, s0, s1, k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length(); 00524 // ivec p0, p1; 00525 00526 alpha_q.set_size(Nstates, block_length + 1, false); 00527 beta_q.set_size(Nstates, block_length + 1, false); 00528 gamma_q.set_size(2*Nstates, block_length + 1, false); 00529 extrinsic_output.set_size(block_length, false); 00530 denom_q.set_size(block_length + 1, false); 00531 for (k = 0; k <= block_length; k++) { denom_q(k) = -QLLR_MAX; } 00532 00533 if (in_terminated) { terminated = true; } 00534 00535 //Check that Lc = 1.0 00536 it_assert(Lc == 1.0, 00537 "Rec_Syst_Conv_Code::log_decode: This function assumes that Lc = 1.0. Please use proper scaling of the input data"); 00538 00539 //Calculate gamma_q 00540 for (k = 1; k <= block_length; k++) { 00541 kk = k - 1; 00542 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00543 exp_temp0 = 0; 00544 exp_temp1 = 0; 00545 for (j = 0; j < (n - 1); j++) { 00546 rp = rec_parity(kk, j); 00547 if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; } 00548 else { exp_temp0 -= rp; } 00549 if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; } 00550 else { exp_temp1 -= rp; } 00551 } 00552 // right shift cannot be used due to implementation dependancy of how sign is handled? 00553 gamma_q(2*s_prim + 0, k) = ((extrinsic_input(kk) + rec_systematic(kk)) + exp_temp0) / 2; 00554 gamma_q(2*s_prim + 1, k) = - ((extrinsic_input(kk) + rec_systematic(kk)) - exp_temp1) / 2; 00555 } 00556 } 00557 00558 //Initiate alpha_q 00559 for (j = 1; j < Nstates; j++) { alpha_q(j, 0) = -QLLR_MAX; } 00560 alpha_q(0, 0) = 0; 00561 00562 //Calculate alpha_q, going forward through the trellis 00563 for (k = 1; k <= block_length; k++) { 00564 for (s = 0; s < Nstates; s++) { 00565 s_prim0 = rev_state_trans(s, 0); 00566 s_prim1 = rev_state_trans(s, 1); 00567 temp0 = alpha_q(s_prim0, k - 1) + gamma_q(2 * s_prim0 + 0, k); 00568 temp1 = alpha_q(s_prim1, k - 1) + gamma_q(2 * s_prim1 + 1, k); 00569 // alpha_q(s,k) = com_log( temp0, temp1 ); 00570 // denom_q(k) = com_log( alpha_q(s,k), denom_q(k) ); 00571 alpha_q(s, k) = llrcalc.jaclog(temp0, temp1); 00572 denom_q(k) = llrcalc.jaclog(alpha_q(s, k), denom_q(k)); 00573 } 00574 //Normalization of alpha_q 00575 for (l = 0; l < Nstates; l++) { alpha_q(l, k) -= denom_q(k); } 00576 } 00577 00578 //Initiate beta_q 00579 if (terminated) { 00580 for (i = 1; i < Nstates; i++) { beta_q(i, block_length) = -QLLR_MAX; } 00581 beta_q(0, block_length) = 0; 00582 } 00583 else { 00584 beta_q.set_col(block_length, alpha_q.get_col(block_length)); 00585 } 00586 00587 //Calculate beta_q going backward in the trellis 00588 for (k = block_length; k >= 1; k--) { 00589 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00590 s0 = state_trans(s_prim, 0); 00591 s1 = state_trans(s_prim, 1); 00592 // beta_q(s_prim,k-1) = com_log( beta_q(s0,k) + gamma_q(2*s_prim+0,k) , beta_q(s1,k) + gamma_q(2*s_prim+1,k) ); 00593 beta_q(s_prim, k - 1) = llrcalc.jaclog(beta_q(s0, k) + gamma_q(2 * s_prim + 0, k) , beta_q(s1, k) + gamma_q(2 * s_prim + 1, k)); 00594 } 00595 //Normalization of beta_q 00596 for (l = 0; l < Nstates; l++) { beta_q(l, k - 1) -= denom_q(k); } 00597 } 00598 00599 //Calculate extrinsic output for each bit 00600 for (k = 1; k <= block_length; k++) { 00601 kk = k - 1; 00602 nom = -QLLR_MAX; 00603 den = -QLLR_MAX; 00604 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00605 s0 = state_trans(s_prim, 0); 00606 s1 = state_trans(s_prim, 1); 00607 exp_temp0 = 0; 00608 exp_temp1 = 0; 00609 for (j = 0; j < (n - 1); j++) { 00610 rp = rec_parity(kk, j); 00611 if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; } 00612 else { exp_temp0 -= rp; } 00613 if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; } 00614 else { exp_temp1 -= rp; } 00615 } 00616 // nom = com_log(nom, alpha_q(s_prim,kk) + 0.5*exp_temp0 + beta_q(s0,k) ); 00617 // den = com_log(den, alpha_q(s_prim,kk) + 0.5*exp_temp1 + beta_q(s1,k) ); 00618 nom = llrcalc.jaclog(nom, alpha_q(s_prim, kk) + exp_temp0 / 2 + beta_q(s0, k)); 00619 den = llrcalc.jaclog(den, alpha_q(s_prim, kk) + exp_temp1 / 2 + beta_q(s1, k)); 00620 } 00621 extrinsic_output(kk) = nom - den; 00622 } 00623 00624 } 00625 00626 00627 00628 void Rec_Syst_Conv_Code::log_decode_n2(const QLLRvec &rec_systematic, 00629 const QLLRvec &rec_parity, 00630 const QLLRvec &extrinsic_input, 00631 QLLRvec &extrinsic_output, 00632 bool in_terminated) 00633 { 00634 int nom, den, exp_temp0, exp_temp1, rp; 00635 int k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length(); 00636 int ext_info_length = extrinsic_input.length(); 00637 ivec p0, p1; 00638 int ex, norm; 00639 00640 00641 alpha_q.set_size(Nstates, block_length + 1, false); 00642 beta_q.set_size(Nstates, block_length + 1, false); 00643 gamma_q.set_size(2*Nstates, block_length + 1, false); 00644 extrinsic_output.set_size(ext_info_length, false); 00645 //denom.set_size(block_length+1,false); for (k=0; k<=block_length; k++) { denom(k) = -infinity; } 00646 00647 if (in_terminated) { terminated = true; } 00648 00649 //Check that Lc = 1.0 00650 it_assert(Lc == 1.0, 00651 "Rec_Syst_Conv_Code::log_decode_n2: This function assumes that Lc = 1.0. Please use proper scaling of the input data"); 00652 00653 //Initiate alpha 00654 for (s = 1; s < Nstates; s++) { alpha_q(s, 0) = -QLLR_MAX; } 00655 alpha_q(0, 0) = 0; 00656 00657 //Calculate alpha and gamma going forward through the trellis 00658 for (k = 1; k <= block_length; k++) { 00659 kk = k - 1; 00660 if (kk < ext_info_length) { 00661 ex = (extrinsic_input(kk) + rec_systematic(kk)) / 2; 00662 } 00663 else { 00664 ex = rec_systematic(kk) / 2; 00665 } 00666 rp = rec_parity(kk) / 2; 00667 for (s = 0; s < Nstates; s++) { 00668 s_prim0 = rev_state_trans(s, 0); 00669 s_prim1 = rev_state_trans(s, 1); 00670 if (output_parity(s_prim0 , 0)) { exp_temp0 = -rp; } 00671 else { exp_temp0 = rp; } 00672 if (output_parity(s_prim1 , 1)) { exp_temp1 = -rp; } 00673 else { exp_temp1 = rp; } 00674 gamma_q(2*s_prim0 , k) = ex + exp_temp0; 00675 gamma_q(2*s_prim1 + 1, k) = -ex + exp_temp1; 00676 alpha_q(s, k) = llrcalc.jaclog(alpha_q(s_prim0, kk) + gamma_q(2 * s_prim0 , k), 00677 alpha_q(s_prim1, kk) + gamma_q(2 * s_prim1 + 1, k)); 00678 //denom(k) = com_log( alpha(s,k), denom(k) ); 00679 } 00680 norm = alpha_q(0, k); //norm = denom(k); 00681 for (l = 0; l < Nstates; l++) { alpha_q(l, k) -= norm; } 00682 } 00683 00684 //Initiate beta 00685 if (terminated) { 00686 for (s = 1; s < Nstates; s++) { beta_q(s, block_length) = -QLLR_MAX; } 00687 beta_q(0, block_length) = 0; 00688 } 00689 else { 00690 beta_q.set_col(block_length, alpha_q.get_col(block_length)); 00691 } 00692 00693 //Calculate beta going backward in the trellis 00694 for (k = block_length; k >= 1; k--) { 00695 kk = k - 1; 00696 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00697 beta_q(s_prim, kk) = llrcalc.jaclog(beta_q(state_trans(s_prim, 0), k) + gamma_q(2 * s_prim, k), 00698 beta_q(state_trans(s_prim, 1), k) + gamma_q(2 * s_prim + 1, k)); 00699 } 00700 norm = beta_q(0, k); //norm = denom(k); 00701 for (l = 0; l < Nstates; l++) { beta_q(l, k) -= norm; } 00702 } 00703 00704 //Calculate extrinsic output for each bit 00705 for (k = 1; k <= block_length; k++) { 00706 kk = k - 1; 00707 if (kk < ext_info_length) { 00708 nom = -QLLR_MAX; 00709 den = -QLLR_MAX; 00710 rp = rec_parity(kk) / 2; 00711 for (s_prim = 0; s_prim < Nstates; s_prim++) { 00712 if (output_parity(s_prim , 0)) { exp_temp0 = -rp; } 00713 else { exp_temp0 = rp; } 00714 if (output_parity(s_prim , 1)) { exp_temp1 = -rp; } 00715 else { exp_temp1 = rp; } 00716 nom = llrcalc.jaclog(nom, alpha_q(s_prim, kk) + exp_temp0 + beta_q(state_trans(s_prim, 0), k)); 00717 den = llrcalc.jaclog(den, alpha_q(s_prim, kk) + exp_temp1 + beta_q(state_trans(s_prim, 1), k)); 00718 } 00719 extrinsic_output(kk) = nom - den; 00720 } 00721 } 00722 } 00723 00724 void Rec_Syst_Conv_Code::set_llrcalc(LLR_calc_unit in_llrcalc) 00725 { 00726 llrcalc = in_llrcalc; 00727 } 00728 00729 00730 } // namespace itpp
Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4