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