IT++ Logo Newcom Logo

rec_syst_conv_code.cpp

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

Generated on Sat Aug 25 23:37:28 2007 for IT++ by Doxygen 1.5.2