IT++ Logo

rec_syst_conv_code.cpp

Go to the documentation of this file.
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
SourceForge Logo

Generated on Thu Apr 23 20:06:47 2009 for IT++ by Doxygen 1.5.8