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 Sat May 3 16:10:43 2008 for IT++ by Doxygen 1.5.5