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