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 
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
SourceForge Logo

Generated on Sun Dec 9 17:31:03 2007 for IT++ by Doxygen 1.5.4