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