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

Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4