IT++ Logo

convcode.cpp

Go to the documentation of this file.
00001 
00030 #include <itpp/comm/convcode.h>
00031 #include <itpp/base/binary.h>
00032 #include <itpp/base/matfunc.h>
00033 #include <limits>
00034 
00035 namespace itpp {
00036 
00037   // ----------------- Protected functions -----------------------------
00038 
00039   /*
00040     The weight of the transition from given state with the input given
00041   */
00042   int Convolutional_Code::weight(const int state, const int input)
00043   {
00044     int i, j, temp, out, w = 0, shiftreg = state;
00045 
00046     shiftreg = shiftreg | (input << m);
00047     for (j=0; j<n; j++) {
00048       out = 0;
00049       temp = shiftreg & gen_pol(j);
00050       for (i=0; i<K; i++) {
00051         out ^= (temp & 1);
00052         temp = temp >> 1;
00053       }
00054       w += out;
00055       //w += weight_int_table(temp);
00056     }
00057     return w;
00058   }
00059 
00060   /*
00061     The weight (of the reverse code) of the transition from given state with
00062     the input given
00063   */
00064   int Convolutional_Code::weight_reverse(const int state, const int input)
00065   {
00066     int i, j, temp, out, w = 0, shiftreg = state;
00067 
00068     shiftreg = shiftreg | (input << m);
00069     for (j=0; j<n; j++) {
00070       out = 0;
00071       temp = shiftreg & gen_pol_rev(j);
00072       for (i=0; i<K; i++) {
00073         out ^= (temp & 1);
00074         temp = temp >> 1;
00075       }
00076       w += out;
00077       //w += weight_int_table(temp);
00078     }
00079     return w;
00080   }
00081 
00082   /*
00083     The weight of the two paths (input 0 or 1) from given state
00084   */
00085   void Convolutional_Code::weight(const int state, int &w0, int &w1)
00086   {
00087     int i, j, temp, out, shiftreg = state;
00088     w0 = 0; w1 = 0;
00089 
00090     shiftreg = shiftreg | (1 << m);
00091     for (j=0; j<n; j++) {
00092       out = 0;
00093       temp = shiftreg & gen_pol(j);
00094 
00095       for (i=0; i<m; i++) {
00096         out ^= (temp & 1);
00097         temp = temp >> 1;
00098       }
00099       w0 += out;
00100       w1 += out ^ (temp & 1);
00101     }
00102   }
00103 
00104   /*
00105     The weight (of the reverse code) of the two paths (input 0 or 1) from
00106     given state
00107   */
00108   void Convolutional_Code::weight_reverse(const int state, int &w0, int &w1)
00109   {
00110     int i, j, temp, out, shiftreg = state;
00111     w0 = 0; w1 = 0;
00112 
00113     shiftreg = shiftreg | (1 << m);
00114     for (j=0; j<n; j++) {
00115       out = 0;
00116       temp = shiftreg & gen_pol_rev(j);
00117 
00118       for (i=0; i<m; i++) {
00119         out ^= (temp & 1);
00120         temp = temp >> 1;
00121       }
00122       w0 += out;
00123       w1 += out ^ (temp & 1);
00124     }
00125   }
00126 
00127   /*
00128     Output on transition (backwards) with input from state
00129   */
00130   bvec Convolutional_Code::output_reverse(const int state, const int input)
00131   {
00132     int j, temp, temp_state;
00133 
00134     bvec output(n);
00135 
00136     temp_state = (state<<1) | input;
00137     for (j=0; j<n; j++) {
00138       temp = temp_state & gen_pol(j);
00139       output(j) = xor_int_table(temp);
00140     }
00141 
00142     return output;
00143   }
00144 
00145   /*
00146     Output on transition (backwards) with input from state
00147   */
00148   void Convolutional_Code::output_reverse(const int state, bvec &zero_output,
00149                                           bvec &one_output)
00150   {
00151     int j, temp, temp_state;
00152     bin one_bit;
00153 
00154     temp_state = (state<<1) | 1;
00155     for (j=0; j<n; j++) {
00156       temp = temp_state & gen_pol(j);
00157       one_bit = temp & 1;
00158       temp = temp >> 1;
00159       one_output(j) = xor_int_table(temp) ^ one_bit;
00160       zero_output(j) = xor_int_table(temp);
00161     }
00162   }
00163 
00164   /*
00165     Output on transition (backwards) with input from state
00166   */
00167   void Convolutional_Code::output_reverse(const int state, int &zero_output,
00168                                           int &one_output)
00169   {
00170     int j, temp, temp_state;
00171     bin one_bit;
00172 
00173     zero_output=0, one_output=0;
00174     temp_state = (state<<1) | 1;
00175     for (j=0; j<n; j++) {
00176       temp = temp_state & gen_pol(j);
00177       one_bit = temp & 1;
00178       temp = temp >> 1;
00179 
00180       one_output = (one_output<<1) | int(xor_int_table(temp) ^ one_bit);
00181       zero_output = (zero_output<<1) | int(xor_int_table(temp));
00182     }
00183   }
00184 
00185   /*
00186     Output on transition (backwards) with input from state
00187   */
00188   void Convolutional_Code::calc_metric_reverse(int state,
00189                                                const vec &rx_codeword,
00190                                                double &zero_metric,
00191                                                double &one_metric)
00192   {
00193     int temp;
00194     bin one_bit;
00195     one_metric = zero_metric = 0;
00196 
00197     int temp_state = (state << 1) | 1;
00198     for (int j = 0; j < n; j++) {
00199       temp = temp_state & gen_pol(j);
00200       one_bit = temp & 1;
00201       temp >>= 1;
00202       one_metric += (2 * static_cast<int>(xor_int_table(temp) ^ one_bit) - 1)
00203         * rx_codeword(j);
00204       zero_metric += (2 * static_cast<int>(xor_int_table(temp)) - 1)
00205         * rx_codeword(j);
00206     }
00207   }
00208 
00209 
00210   // calculates metrics for all codewords (2^n of them) in natural order
00211   void Convolutional_Code::calc_metric(const vec &rx_codeword,
00212                                        vec &delta_metrics)
00213   {
00214     int no_outputs = pow2i(n), no_loop = pow2i(n - 1), mask = no_outputs - 1,
00215       temp;
00216     delta_metrics.set_size(no_outputs, false);
00217 
00218     if (no_outputs <= no_states) {
00219       for (int i = 0; i < no_loop; i++) { // all input possibilities
00220         delta_metrics(i) = 0;
00221         temp = i;
00222         for (int j = n - 1; j >= 0; j--) {
00223           if (temp & 1)
00224             delta_metrics(i) += rx_codeword(j);
00225           else
00226             delta_metrics(i) -= rx_codeword(j);
00227           temp >>= 1;
00228         }
00229         delta_metrics(i ^ mask) = -delta_metrics(i); // the inverse codeword
00230       }
00231     }
00232     else {
00233       double zero_metric, one_metric;
00234       int zero_output, one_output, temp_state;
00235       bin one_bit;
00236       for (int s = 0; s < no_states; s++) { // all states
00237         zero_metric = 0, one_metric = 0;
00238         zero_output = 0, one_output = 0;
00239         temp_state = (s << 1) | 1;
00240         for (int j = 0; j < n; j++) {
00241           temp = temp_state & gen_pol(j);
00242           one_bit = temp & 1;
00243           temp >>= 1;
00244           if (xor_int_table(temp)) {
00245             zero_metric += rx_codeword(j);
00246             one_metric -= rx_codeword(j);
00247           }
00248           else {
00249             zero_metric -= rx_codeword(j);
00250             one_metric += rx_codeword(j);
00251           }
00252           one_output = (one_output << 1)
00253             | static_cast<int>(xor_int_table(temp) ^ one_bit);
00254           zero_output = (zero_output << 1)
00255             | static_cast<int>(xor_int_table(temp));
00256         }
00257         delta_metrics(zero_output) = zero_metric;
00258         delta_metrics(one_output) = one_metric;
00259       }
00260     }
00261   }
00262 
00264 
00265   // MFD codes R=1/2
00266   int Conv_Code_MFD_2[15][2] = {
00267     {0,0},
00268     {0,0},
00269     {0,0},
00270     {05,07},
00271     {015,017},
00272     {023,035},
00273     {053,075},
00274     {0133,0171},
00275     {0247,0371},
00276     {0561,0753},
00277     {01167,01545},
00278     {02335,03661},
00279     {04335,05723},
00280     {010533,017661},
00281     {021675,027123},
00282   };
00283 
00284   // MFD codes R=1/3
00285   int Conv_Code_MFD_3[15][3] = {
00286     {0,0,0},
00287     {0,0,0},
00288     {0,0,0},
00289     {05,07,07},
00290     {013,015,017},
00291     {025,033,037},
00292     {047,053,075},
00293     {0133,0145,0175},
00294     {0225,0331,0367},
00295     {0557,0663,0711},
00296     {0117,01365,01633},
00297     {02353,02671,03175},
00298     {04767,05723,06265},
00299     {010533,010675,017661},
00300     {021645,035661,037133},
00301   };
00302 
00303   // MFD codes R=1/4
00304   int Conv_Code_MFD_4[15][4] = {
00305     {0,0,0,0},
00306     {0,0,0,0},
00307     {0,0,0,0},
00308     {05,07,07,07},
00309     {013,015,015,017},
00310     {025,027,033,037},
00311     {053,067,071,075},
00312     {0135,0135,0147,0163},
00313     {0235,0275,0313,0357},
00314     {0463,0535,0733,0745},
00315     {0117,01365,01633,01653},
00316     {02327,02353,02671,03175},
00317     {04767,05723,06265,07455},
00318     {011145,012477,015573,016727},
00319     {021113,023175,035527,035537},
00320   };
00321 
00322 
00323   // MFD codes R=1/5
00324   int Conv_Code_MFD_5[9][5] = {
00325     {0,0,0,0,0},
00326     {0,0,0,0,0},
00327     {0,0,0,0,0},
00328     {07,07,07,05,05},
00329     {017,017,013,015,015},
00330     {037,027,033,025,035},
00331     {075,071,073,065,057},
00332     {0175,0131,0135,0135,0147},
00333     {0257,0233,0323,0271,0357},
00334   };
00335 
00336   // MFD codes R=1/6
00337   int Conv_Code_MFD_6[9][6] = {
00338     {0,0,0,0,0,0},
00339     {0,0,0,0,0,0},
00340     {0,0,0,0,0,0},
00341     {07,07,07,07,05,05},
00342     {017,017,013,013,015,015},
00343     {037,035,027,033,025,035},
00344     {073,075,055,065,047,057},
00345     {0173,0151,0135,0135,0163,0137},
00346     {0253,0375,0331,0235,0313,0357},
00347   };
00348 
00349   // MFD codes R=1/7
00350   int Conv_Code_MFD_7[9][7] = {
00351     {0,0,0,0,0,0,0},
00352     {0,0,0,0,0,0,0},
00353     {0,0,0,0,0,0,0},
00354     {07,07,07,07,05,05,05},
00355     {017,017,013,013,013,015,015},
00356     {035,027,025,027,033,035,037},
00357     {053,075,065,075,047,067,057},
00358     {0165,0145,0173,0135,0135,0147,0137},
00359     {0275,0253,0375,0331,0235,0313,0357},
00360   };
00361 
00362   // MFD codes R=1/8
00363   int Conv_Code_MFD_8[9][8] = {
00364     {0,0,0,0,0,0,0,0},
00365     {0,0,0,0,0,0,0,0},
00366     {0,0,0,0,0,0,0,0},
00367     {07,07,05,05,05,07,07,07},
00368     {017,017,013,013,013,015,015,017},
00369     {037,033,025,025,035,033,027,037},
00370     {057,073,051,065,075,047,067,057},
00371     {0153,0111,0165,0173,0135,0135,0147,0137},
00372     {0275,0275,0253,0371,0331,0235,0313,0357},
00373   };
00374 
00375   int maxK_Conv_Code_MFD[9] = {0,0,14,14,14,8,8,8,8};
00376 
00377   void get_MFD_gen_pol(int n, int K, ivec & gen)
00378   {
00379     gen.set_size(n);
00380 
00381     switch (n) {
00382     case 2: // Rate 1/2
00383       it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[2], "This convolutional code doesn't exist in the tables");
00384       gen(0) = Conv_Code_MFD_2[K][0];
00385       gen(1) = Conv_Code_MFD_2[K][1];
00386       break;
00387     case 3: // Rate 1/3
00388       it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[3], "This convolutional code doesn't exist in the tables");
00389       gen(0) = Conv_Code_MFD_3[K][0];
00390       gen(1) = Conv_Code_MFD_3[K][1];
00391       gen(2) = Conv_Code_MFD_3[K][2];
00392       break;
00393     case 4: // Rate 1/4
00394       it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[4], "This convolutional code doesn't exist in the tables");
00395       gen(0) = Conv_Code_MFD_4[K][0];
00396       gen(1) = Conv_Code_MFD_4[K][1];
00397       gen(2) = Conv_Code_MFD_4[K][2];
00398       gen(3) = Conv_Code_MFD_4[K][3];
00399       break;
00400     case 5: // Rate 1/5
00401       it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[5], "This convolutional code doesn't exist in the tables");
00402       gen(0) = Conv_Code_MFD_5[K][0];
00403       gen(1) = Conv_Code_MFD_5[K][1];
00404       gen(2) = Conv_Code_MFD_5[K][2];
00405       gen(3) = Conv_Code_MFD_5[K][3];
00406       gen(4) = Conv_Code_MFD_5[K][4];
00407       break;
00408     case 6: // Rate 1/6
00409       it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[6], "This convolutional code doesn't exist in the tables");
00410       gen(0) = Conv_Code_MFD_6[K][0];
00411       gen(1) = Conv_Code_MFD_6[K][1];
00412       gen(2) = Conv_Code_MFD_6[K][2];
00413       gen(3) = Conv_Code_MFD_6[K][3];
00414       gen(4) = Conv_Code_MFD_6[K][4];
00415       gen(5) = Conv_Code_MFD_6[K][5];
00416       break;
00417     case 7: // Rate 1/7
00418       it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[7], "This convolutional code doesn't exist in the tables");
00419       gen(0) = Conv_Code_MFD_7[K][0];
00420       gen(1) = Conv_Code_MFD_7[K][1];
00421       gen(2) = Conv_Code_MFD_7[K][2];
00422       gen(3) = Conv_Code_MFD_7[K][3];
00423       gen(4) = Conv_Code_MFD_7[K][4];
00424       gen(5) = Conv_Code_MFD_7[K][5];
00425       gen(6) = Conv_Code_MFD_7[K][6];
00426       break;
00427     case 8: // Rate 1/8
00428       it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[8], "This convolutional code doesn't exist in the tables");
00429       gen(0) = Conv_Code_MFD_8[K][0];
00430       gen(1) = Conv_Code_MFD_8[K][1];
00431       gen(2) = Conv_Code_MFD_8[K][2];
00432       gen(3) = Conv_Code_MFD_8[K][3];
00433       gen(4) = Conv_Code_MFD_8[K][4];
00434       gen(5) = Conv_Code_MFD_8[K][5];
00435       gen(6) = Conv_Code_MFD_8[K][6];
00436       gen(7) = Conv_Code_MFD_8[K][7];
00437       break;
00438     default: // wrong rate
00439       it_assert(false, "This convolutional code doesn't exist in the tables");
00440     }
00441   }
00442 
00443   // ODS codes R=1/2
00444   int Conv_Code_ODS_2[17][2] = {
00445     {0,0},
00446     {0,0},
00447     {0,0},
00448     {05,07},
00449     {015,017},
00450     {023,035},
00451     {053,075},
00452     {0133,0171},
00453     {0247,0371},
00454     {0561,0753},
00455     {01151,01753},
00456     {03345,03613},
00457     {05261,07173},
00458     {012767,016461},
00459     {027251,037363},
00460     {063057,044735},
00461     {0126723,0152711},
00462   };
00463 
00464   // ODS codes R=1/3
00465   int Conv_Code_ODS_3[14][3] = {
00466     {0,0,0},
00467     {0,0,0},
00468     {0,0,0},
00469     {05,07,07},
00470     {013,015,017},
00471     {025,033,037},
00472     {047,053,075},
00473     {0133,0165,0171},
00474     {0225,0331,0367},
00475     {0575,0623,0727},
00476     {01233,01375,01671},
00477     {02335,02531,03477},
00478     {05745,06471,07553},
00479     {013261,015167,017451},
00480   };
00481 
00482   // ODS codes R=1/4
00483   int Conv_Code_ODS_4[12][4] = {
00484     {0,0,0,0},
00485     {0,0,0,0},
00486     {0,0,0,0},
00487     {05,05,07,07},
00488     {013,015,015,017},
00489     {025,027,033,037},
00490     {051,055,067,077},
00491     {0117,0127,0155,0171},
00492     {0231,0273,0327,0375},
00493     {0473,0513,0671,0765},
00494     {01173,01325,01467,01751},
00495     {02565,02747,03311,03723},
00496   };
00497 
00498   int maxK_Conv_Code_ODS[5] = {0,0,16,13,11};
00499 
00500   void get_ODS_gen_pol(int n, int K, ivec & gen)
00501   {
00502     gen.set_size(n);
00503 
00504     switch (n) {
00505     case 2: // Rate 1/2
00506       it_assert(K >= 3 && K <= maxK_Conv_Code_ODS[2], "This convolutional code doesn't exist in the tables");
00507       gen(0) = Conv_Code_ODS_2[K][0];
00508       gen(1) = Conv_Code_ODS_2[K][1];
00509       break;
00510     case 3: // Rate 1/3
00511       it_assert(K >= 3 && K <= maxK_Conv_Code_ODS[3], "This convolutional code doesn't exist in the tables");
00512       gen(0) = Conv_Code_ODS_3[K][0];
00513       gen(1) = Conv_Code_ODS_3[K][1];
00514       gen(2) = Conv_Code_ODS_3[K][2];
00515       break;
00516     case 4: // Rate 1/4
00517       it_assert(K >= 3 && K <= maxK_Conv_Code_ODS[4], "This convolutional code doesn't exist in the tables");
00518       gen(0) = Conv_Code_ODS_4[K][0];
00519       gen(1) = Conv_Code_ODS_4[K][1];
00520       gen(2) = Conv_Code_ODS_4[K][2];
00521       gen(3) = Conv_Code_ODS_4[K][3];
00522       break;
00523     default: // wrong rate
00524       it_assert(false, "This convolutional code doesn't exist in the tables");
00525     }
00526   }
00527 
00529 
00530   // --------------- Public functions -------------------------
00531 
00532   void Convolutional_Code::set_code(const CONVOLUTIONAL_CODE_TYPE type_of_code,
00533                                     int inverse_rate, int constraint_length)
00534   {
00535     ivec gen;
00536 
00537     if (type_of_code == MFD)
00538       get_MFD_gen_pol(inverse_rate, constraint_length, gen);
00539     else if (type_of_code == ODS)
00540       get_ODS_gen_pol(inverse_rate, constraint_length, gen);
00541     else
00542       it_assert(false, "This convolutional code doesn't exist in the tables");
00543 
00544     set_generator_polynomials(gen, constraint_length);
00545   }
00546 
00547   /*
00548     Set generator polynomials. Given in Proakis integer form
00549   */
00550   void Convolutional_Code::set_generator_polynomials(const ivec &gen,
00551                                                      int constraint_length)
00552   {
00553     it_error_if(constraint_length <= 0, "Convolutional_Code::set_generator_polynomials(): Constraint length out of range");
00554     gen_pol = gen;
00555     n = gen.size();
00556     it_error_if(n <= 0, "Convolutional_Code::set_generator_polynomials(): Invalid code rate");
00557 
00558     // Generate table look-up of weight of integers of size K bits
00559     if (constraint_length != K) {
00560       K = constraint_length;
00561       xor_int_table.set_size(pow2i(K), false);
00562       for (int i = 0; i < pow2i(K); i++) {
00563         xor_int_table(i) = (weight_int(K, i) & 1);
00564       }
00565     }
00566 
00567     trunc_length = 5 * K;
00568     m = K - 1;
00569     no_states = pow2i(m);
00570     gen_pol_rev.set_size(n, false);
00571     rate = 1.0 / n;
00572 
00573     for (int i = 0; i < n; i++) {
00574       gen_pol_rev(i) = reverse_int(K, gen_pol(i));
00575     }
00576 
00577     int zero_output, one_output;
00578     output_reverse_int.set_size(no_states, 2, false);
00579 
00580     for (int i = 0; i < no_states; i++) {
00581       output_reverse(i, zero_output, one_output);
00582       output_reverse_int(i, 0) = zero_output;
00583       output_reverse_int(i, 1) = one_output;
00584     }
00585 
00586     // initialise memory structures
00587     visited_state.set_size(no_states);
00588     visited_state = false;
00589     visited_state(start_state) = true; // mark start state
00590 
00591     sum_metric.set_size(no_states);
00592     sum_metric.clear();
00593 
00594     trunc_ptr = 0;
00595     trunc_state = 0;
00596 
00597   }
00598 
00599   // Reset encoder and decoder states
00600   void Convolutional_Code::reset()
00601   {
00602     init_encoder();
00603 
00604     visited_state = false;
00605     visited_state(start_state) = true; // mark start state
00606 
00607     sum_metric.clear();
00608 
00609     trunc_ptr = 0;
00610     trunc_state = 0;
00611   }
00612 
00613 
00614   /*
00615     Encode a binary vector of inputs using specified method
00616   */
00617   void Convolutional_Code::encode(const bvec &input, bvec &output)
00618   {
00619     switch (cc_method) {
00620     case Trunc:
00621       encode_trunc(input, output);
00622       break;
00623     case Tailbite:
00624       encode_tailbite(input, output);
00625       break;
00626     case Tail:
00627     default:
00628       encode_tail(input, output);
00629       break;
00630     };
00631   }
00632 
00633   /*
00634     Encode a binary vector of inputs starting from state set by the
00635     set_state function
00636   */
00637   void Convolutional_Code::encode_trunc(const bvec &input, bvec &output)
00638   {
00639     int temp;
00640     output.set_size(input.size() * n, false);
00641 
00642     for (int i = 0; i < input.size(); i++) {
00643       encoder_state |=  static_cast<int>(input(i)) << m;
00644       for (int j = 0; j < n; j++) {
00645         temp = encoder_state & gen_pol(j);
00646         output(i * n + j) = xor_int_table(temp);
00647       }
00648       encoder_state >>= 1;
00649     }
00650   }
00651 
00652   /*
00653     Encode a binary vector of inputs starting from zero state and also adds
00654     a tail of K-1 zeros to force the encoder into the zero state. Well
00655     suited for packet transmission.
00656   */
00657   void Convolutional_Code::encode_tail(const bvec &input, bvec &output)
00658   {
00659     int temp;
00660     output.set_size((input.size() + m) * n, false);
00661 
00662     // always start from state 0
00663     encoder_state = 0;
00664 
00665     for (int i = 0; i < input.size(); i++) {
00666       encoder_state |=  static_cast<int>(input(i)) << m;
00667       for (int j = 0; j < n; j++) {
00668         temp = encoder_state & gen_pol(j);
00669         output(i * n + j) = xor_int_table(temp);
00670       }
00671       encoder_state >>= 1;
00672     }
00673 
00674     // add tail of m = K-1 zeros
00675     for (int i = input.size(); i < input.size() + m; i++) {
00676       for (int j = 0; j < n; j++) {
00677         temp = encoder_state & gen_pol(j);
00678         output(i * n + j) = xor_int_table(temp);
00679       }
00680       encoder_state >>= 1;
00681     }
00682   }
00683 
00684   /*
00685     Encode a binary vector of inputs starting using tail biting
00686   */
00687   void Convolutional_Code::encode_tailbite(const bvec &input, bvec &output)
00688   {
00689     int temp;
00690     output.set_size(input.size() * n, false);
00691 
00692     // Set the start state equal to the end state:
00693     encoder_state = 0;
00694     bvec last_bits = input.right(m);
00695     for (int i = 0; i < m; i++) {
00696       encoder_state |= static_cast<int>(last_bits(i)) << m;
00697       encoder_state >>= 1;
00698     }
00699 
00700     for (int i = 0; i < input.size(); i++) {
00701       encoder_state |= static_cast<int>(input(i)) << m;
00702       for (int j = 0; j < n; j++) {
00703         temp = encoder_state & gen_pol(j);
00704         output(i * n + j) = xor_int_table(temp);
00705       }
00706       encoder_state >>= 1;
00707     }
00708   }
00709 
00710   /*
00711     Encode a binary bit starting from the internal encoder state.
00712     To initialize the encoder state use set_start_state() and init_encoder()
00713   */
00714   void Convolutional_Code::encode_bit(const bin &input, bvec &output)
00715   {
00716     int temp;
00717     output.set_size(n, false);
00718 
00719     encoder_state |= static_cast<int>(input) << m;
00720     for (int j = 0; j < n; j++) {
00721       temp = encoder_state & gen_pol(j);
00722       output(j) = xor_int_table(temp);
00723     }
00724     encoder_state >>= 1;
00725   }
00726 
00727 
00728   // --------------- Hard-decision decoding is not implemented -----------------
00729 
00730   void Convolutional_Code::decode(const bvec &coded_bits, bvec &output)
00731   {
00732     it_error("Convolutional_Code::decode(): Hard-decision decoding not implemented");
00733   }
00734 
00735   bvec Convolutional_Code::decode(const bvec &coded_bits)
00736   {
00737     it_error("Convolutional_Code::decode(): Hard-decision decoding not implemented");
00738     return bvec();
00739   }
00740 
00741 
00742   /*
00743     Decode a block of encoded data using specified method
00744   */
00745   void Convolutional_Code::decode(const vec &received_signal, bvec &output)
00746   {
00747     switch (cc_method) {
00748     case Trunc:
00749       decode_trunc(received_signal, output);
00750       break;
00751     case Tailbite:
00752       decode_tailbite(received_signal, output);
00753       break;
00754     case Tail:
00755     default:
00756       decode_tail(received_signal, output);
00757       break;
00758     };
00759   }
00760 
00761   /*
00762     Decode a block of encoded data where encode_tail has been used.
00763     Thus is assumes a decoder start state of zero and that a tail of
00764     K-1 zeros has been added. No memory truncation.
00765   */
00766   void Convolutional_Code::decode_tail(const vec &received_signal, bvec &output)
00767   {
00768     int block_length = received_signal.size() / n; // no input symbols
00769     it_error_if(block_length - m <= 0,
00770                 "Convolutional_Code::decode_tail(): Input sequence to short");
00771     int S0, S1;
00772     vec temp_sum_metric(no_states), temp_rec(n), delta_metrics;
00773     Array<bool> temp_visited_state(no_states);
00774     double temp_metric_zero, temp_metric_one;
00775 
00776     path_memory.set_size(no_states, block_length, false);
00777     output.set_size(block_length - m, false);    // no tail in the output
00778 
00779     // clear visited states
00780     visited_state = false;
00781     temp_visited_state = visited_state;
00782     visited_state(0) = true; // starts in the zero state
00783 
00784     // clear accumulated metrics
00785     sum_metric.clear();
00786 
00787     // evolve up to m where not all states visited
00788     for (int l = 0; l < m; l++) { // all transitions including the tail
00789       temp_rec = received_signal.mid(l * n, n);
00790 
00791       // calculate all metrics for all codewords at the same time
00792       calc_metric(temp_rec, delta_metrics);
00793 
00794       for (int s = 0; s < no_states; s++) { // all states
00795         // S0 and S1 are the states that expanded end at state s
00796         previous_state(s, S0, S1);
00797         if (visited_state(S0)) { // expand trellis if state S0 is visited
00798           temp_metric_zero = sum_metric(S0)
00799             + delta_metrics(output_reverse_int(s, 0));
00800           temp_visited_state(s) = true;
00801         }
00802         else {
00803           temp_metric_zero = std::numeric_limits<double>::max();
00804         }
00805         if (visited_state(S1)) { // expand trellis if state S0 is visited
00806           temp_metric_one = sum_metric(S1)
00807             + delta_metrics(output_reverse_int(s, 1));
00808           temp_visited_state(s) = true;
00809         }
00810         else {
00811           temp_metric_one = std::numeric_limits<double>::max();
00812         }
00813         if (temp_metric_zero < temp_metric_one) { // path zero survives
00814           temp_sum_metric(s) = temp_metric_zero;
00815           path_memory(s, l) = 0;
00816         }
00817         else { // path one survives
00818           temp_sum_metric(s) = temp_metric_one;
00819           path_memory(s, l) = 1;
00820         }
00821       } // all states, s
00822       sum_metric = temp_sum_metric;
00823       visited_state = temp_visited_state;
00824     } // all transitions, l
00825 
00826     // evolve from m and to the end of the block
00827     for (int l = m; l < block_length; l++) { // all transitions except the tail
00828       temp_rec = received_signal.mid(l * n, n);
00829 
00830       // calculate all metrics for all codewords at the same time
00831       calc_metric(temp_rec, delta_metrics);
00832 
00833       for (int s = 0; s < no_states; s++) { // all states
00834         // S0 and S1 are the states that expanded end at state s
00835         previous_state(s, S0, S1);
00836         temp_metric_zero = sum_metric(S0)
00837           + delta_metrics(output_reverse_int(s, 0));
00838         temp_metric_one = sum_metric(S1)
00839           + delta_metrics(output_reverse_int(s, 1));
00840         if (temp_metric_zero < temp_metric_one) { // path zero survives
00841           temp_sum_metric(s) = temp_metric_zero;
00842           path_memory(s, l) = 0;
00843         }
00844         else { // path one survives
00845           temp_sum_metric(s) = temp_metric_one;
00846           path_memory(s, l) = 1;
00847         }
00848       } // all states, s
00849       sum_metric = temp_sum_metric;
00850     } // all transitions, l
00851 
00852     // minimum metric is the zeroth state due to the tail
00853     int min_metric_state = 0;
00854     // trace back to remove tail of zeros
00855     for (int l = block_length - 1; l > block_length - 1 - m; l--) {
00856       // previous state calculation
00857       min_metric_state = previous_state(min_metric_state,
00858                                         path_memory(min_metric_state, l));
00859     }
00860 
00861     // trace back to calculate output sequence
00862     for (int l = block_length - 1 - m; l >= 0; l--) {
00863       output(l) = get_input(min_metric_state);
00864       // previous state calculation
00865       min_metric_state = previous_state(min_metric_state,
00866                                         path_memory(min_metric_state, l));
00867     }
00868   }
00869 
00870 
00871   /*
00872     Decode a block of encoded data where encode_tailbite has been used.
00873   */
00874   void Convolutional_Code::decode_tailbite(const vec &received_signal,
00875                                            bvec &output)
00876   {
00877     int block_length = received_signal.size() / n; // no input symbols
00878     it_error_if(block_length <= 0,
00879                 "Convolutional_Code::decode_tailbite(): Input sequence to short");
00880     int S0, S1;
00881     vec temp_sum_metric(no_states), temp_rec(n), delta_metrics;
00882     Array<bool> temp_visited_state(no_states);
00883     double temp_metric_zero, temp_metric_one;
00884     double best_metric = std::numeric_limits<double>::max();
00885     bvec best_output(block_length), temp_output(block_length);
00886 
00887     path_memory.set_size(no_states, block_length, false);
00888     output.set_size(block_length, false);
00889 
00890 
00891     // Try all start states ss
00892     for (int ss = 0; ss < no_states; ss++) {
00893       //Clear the visited_state vector:
00894       visited_state = false;
00895       temp_visited_state = visited_state;
00896       visited_state(ss) = true; // starts in the state s
00897 
00898       // clear accumulated metrics
00899       sum_metric.zeros();
00900 
00901       for (int l = 0; l < block_length; l++) { // all transitions
00902         temp_rec = received_signal.mid(l * n, n);
00903         // calculate all metrics for all codewords at the same time
00904         calc_metric(temp_rec, delta_metrics);
00905 
00906         for (int s = 0; s < no_states; s++) { // all states
00907           // S0 and S1 are the states that expanded end at state s
00908           previous_state(s, S0, S1);
00909           if (visited_state(S0)) { // expand trellis if state S0 is visited
00910             temp_metric_zero = sum_metric(S0)
00911               + delta_metrics(output_reverse_int(s, 0));
00912             temp_visited_state(s) = true;
00913           }
00914           else {
00915             temp_metric_zero = std::numeric_limits<double>::max();
00916           }
00917           if (visited_state(S1)) { // expand trellis if state S0 is visited
00918             temp_metric_one = sum_metric(S1)
00919               + delta_metrics(output_reverse_int(s, 1));
00920             temp_visited_state(s) = true;
00921           }
00922           else {
00923             temp_metric_one = std::numeric_limits<double>::max();
00924           }
00925           if (temp_metric_zero < temp_metric_one) { // path zero survives
00926             temp_sum_metric(s) = temp_metric_zero;
00927             path_memory(s, l) = 0;
00928           }
00929           else { // path one survives
00930             temp_sum_metric(s) = temp_metric_one;
00931             path_memory(s, l) = 1;
00932           }
00933         } // all states, s
00934         sum_metric = temp_sum_metric;
00935         visited_state = temp_visited_state;
00936       } // all transitions, l
00937 
00938       // minimum metric is the ss state due to the tailbite
00939       int min_metric_state = ss;
00940 
00941       // trace back to calculate output sequence
00942       for (int l = block_length - 1; l >= 0; l--) {
00943         temp_output(l) = get_input(min_metric_state);
00944         // previous state calculation
00945         min_metric_state = previous_state(min_metric_state,
00946                                           path_memory(min_metric_state, l));
00947       }
00948       if (sum_metric(ss) < best_metric) {
00949         best_metric = sum_metric(ss);
00950         best_output = temp_output;
00951       }
00952     } // all start states ss
00953     output = best_output;
00954   }
00955 
00956 
00957   /*
00958     Viterbi decoding using truncation of memory (default = 5*K)
00959   */
00960   void Convolutional_Code::decode_trunc(const vec &received_signal,
00961                                         bvec &output)
00962   {
00963     int block_length = received_signal.size() / n; // no input symbols
00964     it_error_if(block_length <= 0,
00965                 "Convolutional_Code::decode_trunc(): Input sequence to short");
00966     int S0, S1;
00967     vec temp_sum_metric(no_states), temp_rec(n), delta_metrics;
00968     Array<bool> temp_visited_state(no_states);
00969     double temp_metric_zero, temp_metric_one;
00970 
00971     path_memory.set_size(no_states, trunc_length, false);
00972     output.set_size(0);
00973 
00974     // copy visited states
00975     temp_visited_state = visited_state;
00976 
00977     for (int i = 0; i < block_length; i++) {
00978       // update path memory pointer
00979       trunc_ptr = (trunc_ptr + 1) % trunc_length;
00980 
00981       temp_rec = received_signal.mid(i * n, n);
00982       // calculate all metrics for all codewords at the same time
00983       calc_metric(temp_rec, delta_metrics);
00984 
00985       for (int s = 0; s < no_states; s++) { // all states
00986         // the states that expanded end at state s
00987         previous_state(s, S0, S1);
00988         if (visited_state(S0)) { // expand trellis if state S0 is visited
00989           temp_metric_zero = sum_metric(S0)
00990             + delta_metrics(output_reverse_int(s, 0));
00991           temp_visited_state(s) = true;
00992         }
00993         else {
00994           temp_metric_zero = std::numeric_limits<double>::max();
00995         }
00996         if (visited_state(S1)) { // expand trellis if state S0 is visited
00997           temp_metric_one = sum_metric(S1)
00998             + delta_metrics(output_reverse_int(s, 1));
00999           temp_visited_state(s) = true;
01000         }
01001         else {
01002           temp_metric_one = std::numeric_limits<double>::max();
01003         }
01004         if (temp_metric_zero < temp_metric_one) { // path zero survives
01005           temp_sum_metric(s) = temp_metric_zero;
01006           path_memory(s, trunc_ptr) = 0;
01007         }
01008         else { // path one survives
01009           temp_sum_metric(s) = temp_metric_one;
01010           path_memory(s, trunc_ptr) = 1;
01011         }
01012       } // all states, s
01013       sum_metric = temp_sum_metric;
01014       visited_state = temp_visited_state;
01015 
01016       // find minimum metric
01017       int min_metric_state = min_index(sum_metric);
01018 
01019       // normalise accumulated metrics
01020       sum_metric -= sum_metric(min_metric_state);
01021 
01022       // check if we had enough metrics to generate output
01023       if (trunc_state >= trunc_length) {
01024         // trace back to calculate the output symbol
01025         for (int j = trunc_length; j > 0; j--) {
01026           // previous state calculation
01027           min_metric_state =
01028             previous_state(min_metric_state,
01029                            path_memory(min_metric_state,
01030                                        (trunc_ptr + j) % trunc_length));
01031         }
01032         output.ins(output.size(), get_input(min_metric_state));
01033       }
01034       else { // if not increment trunc_state counter
01035         trunc_state++;
01036       }
01037     } // end for (int i = 0; i < block_length; l++)
01038   }
01039 
01040 
01041   /*
01042     Calculate the inverse sequence
01043 
01044     Assumes that encode_tail is used in the encoding process. Returns false
01045     if there is an error in the coded sequence (not a valid codeword). Do
01046     not check that the tail forces the encoder into the zeroth state.
01047   */
01048   bool Convolutional_Code::inverse_tail(const bvec coded_sequence, bvec &input)
01049   {
01050     int state = 0, zero_state, one_state, zero_temp, one_temp, i, j;
01051     bvec zero_output(n), one_output(n);
01052 
01053     int block_length = coded_sequence.size()/n - m; // no input symbols
01054     it_error_if(block_length<=0, "The input sequence is to short");
01055     input.set_length(block_length, false); // not include the tail in the output
01056 
01057 
01058     for (i=0; i<block_length; i++) {
01059       zero_state = state;
01060       one_state = state | (1 << m);
01061       for (j=0; j<n; j++) {
01062         zero_temp = zero_state & gen_pol(j);
01063         one_temp = one_state & gen_pol(j);
01064         zero_output(j) = xor_int_table(zero_temp);
01065         one_output(j) = xor_int_table(one_temp);
01066       }
01067       if (coded_sequence.mid(i*n, n) == zero_output) {
01068         input(i) = bin(0);
01069         state = zero_state >> 1;
01070       } else if (coded_sequence.mid(i*n, n) == one_output) {
01071         input(i) = bin(1);
01072         state = one_state >> 1;
01073       } else
01074         return false;
01075     }
01076 
01077     return true;
01078   }
01079 
01080   /*
01081     Check if catastrophic. Returns true if catastrophic
01082   */
01083   bool Convolutional_Code::catastrophic(void)
01084   {
01085     int start, S, W0, W1, S0, S1;
01086     bvec visited(1<<m);
01087 
01088     for (start=1; start<no_states; start++) {
01089       visited.zeros();
01090       S = start;
01091       visited(S) = 1;
01092 
01093     node1:
01094       S0 = next_state(S, 0);
01095       S1 = next_state(S, 1);
01096       weight(S, W0, W1);
01097       if (S1 < start) goto node0;
01098       if (W1 > 0) goto node0;
01099 
01100       if (visited(S1) == bin(1))
01101         return true;
01102       S = S1; // goto node1
01103       visited(S) = 1;
01104       goto node1;
01105 
01106     node0:
01107       if (S0 < start) continue; // goto end;
01108       if (W0 > 0) continue; // goto end;
01109 
01110       if (visited(S0) == bin(1))
01111         return true;
01112       S = S0;
01113       visited(S) = 1;
01114       goto node1;
01115 
01116       //end:
01117       //void;
01118     }
01119 
01120     return false;
01121   }
01122 
01123   /*
01124     Calculate distance profile. If reverse = true calculate for the reverse code instead.
01125   */
01126   void Convolutional_Code::distance_profile(ivec &dist_prof, int dmax, bool reverse)
01127   {
01128     int max_stack_size = 50000;
01129     ivec S_stack(max_stack_size), W_stack(max_stack_size), t_stack(max_stack_size);
01130 
01131     int stack_pos=-1, t, S, W, W0, w0, w1;
01132 
01133     dist_prof.set_size(K, false);
01134     dist_prof.zeros();
01135     dist_prof += dmax; // just select a big number!
01136     if (reverse)
01137       W = weight_reverse(0, 1);
01138     else
01139       W = weight(0, 1);
01140 
01141     S = next_state(0, 1); // first state 0 and one as input, S = 1<<(m-1);
01142     t = 0;
01143     dist_prof(0) = W;
01144 
01145   node1:
01146     if (reverse)
01147       weight_reverse(S, w0, w1);
01148     else
01149       weight(S, w0, w1);
01150 
01151     if (t < m) {
01152       W0 = W + w0;
01153       if (W0 < dist_prof(m)) { // store node0 (S, W0, and t+1)
01154         stack_pos++;
01155         if (stack_pos >= max_stack_size) {
01156           max_stack_size = 2*max_stack_size;
01157           S_stack.set_size(max_stack_size, true);
01158           W_stack.set_size(max_stack_size, true);
01159           t_stack.set_size(max_stack_size, true);
01160         }
01161 
01162         S_stack(stack_pos) = next_state(S, 0); //S>>1;
01163         W_stack(stack_pos) = W0;
01164         t_stack(stack_pos) = t+1;
01165       }
01166     } else goto stack;
01167 
01168     W += w1;
01169     if (W > dist_prof(m))
01170       goto stack;
01171 
01172     t++;
01173     S = next_state(S, 1); //S = (S>>1)|(1<<(m-1));
01174 
01175     if (W < dist_prof(t))
01176       dist_prof(t) = W;
01177 
01178     if(t == m) goto stack;
01179     goto node1;
01180 
01181 
01182   stack:
01183     if (stack_pos >= 0) {
01184       // get root node from stack
01185       S = S_stack(stack_pos);
01186       W = W_stack(stack_pos);
01187       t = t_stack(stack_pos);
01188       stack_pos--;
01189 
01190       if (W < dist_prof(t))
01191         dist_prof(t) = W;
01192 
01193       if (t == m) goto stack;
01194 
01195       goto node1;
01196     }
01197   }
01198 
01199   /*
01200     Calculate spectrum
01201 
01202     Calculates both the weight spectrum (Ad) and the information weight spectrum (Cd) and
01203     returns it as ivec:s in the 0:th and 1:st component of spectrum, respectively. Suitable
01204     for calculating many terms in the spectra (uses an breadth first algorithm). It is assumed
01205     that the code is non-catastrophic or else it is a possibility for an eternal loop.
01206     dmax = an upper bound on the free distance
01207     no_terms = no_terms including the dmax term that should be calculated
01208   */
01209   void Convolutional_Code::calculate_spectrum(Array<ivec> &spectrum, int dmax, int no_terms)
01210   {
01211     imat Ad_states(no_states, dmax+no_terms), Cd_states(no_states, dmax+no_terms);
01212     imat Ad_temp(no_states, dmax+no_terms), Cd_temp(no_states, dmax+no_terms);
01213     ivec mindist(no_states),  mindist_temp(1<<m);
01214 
01215     spectrum.set_size(2);
01216     spectrum(0).set_size(dmax+no_terms, false);
01217     spectrum(1).set_size(dmax+no_terms, false);
01218     spectrum(0).zeros();
01219     spectrum(1).zeros();
01220     Ad_states.zeros();
01221     Cd_states.zeros();
01222     mindist.zeros();
01223     int wmax = dmax+no_terms;
01224     ivec visited_states(no_states), visited_states_temp(no_states);
01225     bool proceede;
01226     int d, w0, w1, s, s0, s1;
01227 
01228     visited_states.zeros(); // 0 = false
01229     s = next_state(0, 1); // Start in state from 0 with an one input.
01230     visited_states(s) = 1;  // 1 = true. Start in state 2^(m-1).
01231     w1 = weight(0, 1);
01232     Ad_states(s, w1) = 1;
01233     Cd_states(s, w1) = 1;
01234     mindist(s) = w1;
01235     proceede = true;
01236 
01237     while (proceede) {
01238       proceede = false;
01239       Ad_temp.zeros();
01240       Cd_temp.zeros();
01241       mindist_temp.zeros();
01242       visited_states_temp.zeros(); //false
01243       for (s=1; s<no_states; s++) {
01244         if ((mindist(s)>0) && (mindist(s)<wmax)) {
01245           proceede = true;
01246           weight(s,w0,w1);
01247           s0 = next_state(s, 0);
01248           for (d=mindist(s); d<(wmax-w0); d++) {
01249             Ad_temp(s0,d+w0) += Ad_states(s,d);
01250             Cd_temp(s0,d+w0) += Cd_states(s,d);
01251             visited_states_temp(s0) = 1; //true
01252           }
01253 
01254           s1 = next_state(s, 1);
01255           for (d=mindist(s); d<(wmax-w1); d++) {
01256             Ad_temp(s1,d+w1) += Ad_states(s,d);
01257             Cd_temp(s1,d+w1) += Cd_states(s,d) + Ad_states(s,d);
01258             visited_states_temp(s1) = 1; //true
01259           }
01260           if (mindist_temp(s0)>0) { mindist_temp(s0) = ( mindist(s)+w0 ) < mindist_temp(s0) ? mindist(s)+w0 :  mindist_temp(s0); }
01261           else { mindist_temp(s0) = mindist(s)+w0; }
01262           if (mindist_temp(s1)>0) { mindist_temp(s1) = ( mindist(s)+w1 ) < mindist_temp(s1) ? mindist(s)+w1 :  mindist_temp(s1); }
01263           else { mindist_temp(s1) = mindist(s)+w1; }
01264 
01265         }
01266       }
01267       Ad_states = Ad_temp;
01268       Cd_states = Cd_temp;
01269       spectrum(0) += Ad_temp.get_row(0);
01270       spectrum(1) += Cd_temp.get_row(0);
01271       visited_states = visited_states_temp;
01272       mindist = elem_mult(mindist_temp, visited_states);
01273     }
01274   }
01275 
01276   /*
01277     Cederwall's fast algorithm
01278 
01279     See IT No. 6, pp. 1146-1159, Nov. 1989 for details.
01280     Calculates both the weight spectrum (Ad) and the information weight spectrum (Cd) and
01281     returns it as ivec:s in the 0:th and 1:st component of spectrum, respectively. The FAST algorithm
01282     is good for calculating only a few terms in the spectrum. If many terms are desired, use calc_spectrum instead.
01283     The algorithm returns -1 if the code tested is worse that the input dfree and Cdfree.
01284     It returns 0 if the code MAY be catastrophic (assuming that test_catastrophic is true),
01285     and returns 1 if everything went right.
01286 
01287     \arg \c dfree the free distance of the code (or an upper bound)
01288     \arg \c no_terms including the dfree term that should be calculated
01289     \ar \c Cdfree is the best value of information weight spectrum found so far
01290   */
01291   int Convolutional_Code::fast(Array<ivec> &spectrum, const int dfree, const int no_terms, const int Cdfree, const bool test_catastrophic)
01292   {
01293     int cat_treshold = 7*K; // just a big number, but not to big!
01294     int i;
01295     ivec dist_prof(K), dist_prof_rev(K);
01296     distance_profile(dist_prof, dfree);
01297     distance_profile(dist_prof_rev, dfree, true); // for the reverse code
01298 
01299     int dist_prof_rev0 = dist_prof_rev(0);
01300 
01301     bool reverse = false; // false = use dist_prof
01302 
01303     // is the reverse distance profile better?
01304     for (i=0; i<K; i++) {
01305       if (dist_prof_rev(i) > dist_prof(i)) {
01306         reverse = true;
01307         dist_prof_rev0 = dist_prof(0);
01308         dist_prof = dist_prof_rev;
01309         break;
01310       } else if (dist_prof_rev(i) < dist_prof(i)) {
01311         break;
01312       }
01313     }
01314 
01315     int d = dfree + no_terms - 1;
01316     int max_stack_size = 50000;
01317     ivec S_stack(max_stack_size), W_stack(max_stack_size), b_stack(max_stack_size), c_stack(max_stack_size);
01318     int stack_pos=-1;
01319 
01320     // F1.
01321     int mf = 1, b = 1;
01322     spectrum.set_size(2);
01323     spectrum(0).set_size(dfree+no_terms, false); // ad
01324     spectrum(1).set_size(dfree+no_terms, false); // cd
01325     spectrum(0).zeros();
01326     spectrum(1).zeros();
01327     int S, S0, S1, W0, W1, w0, w1, c = 0;
01328     S = next_state(0, 1); //first state zero and one as input
01329     int W = d - dist_prof_rev0;
01330 
01331 
01332   F2:
01333     S0 = next_state(S, 0);
01334     S1 = next_state(S, 1);
01335 
01336     if (reverse) {
01337       weight(S, w0, w1);
01338     } else {
01339       weight_reverse(S, w0, w1);
01340     }
01341     W0 = W - w0;
01342     W1 = W - w1;
01343     if(mf < m) goto F6;
01344 
01345     //F3:
01346     if (W0 >= 0) {
01347       spectrum(0)(d-W0)++;
01348       spectrum(1)(d-W0) += b;
01349       // The code is worse than the best found.
01350       if ( ((d-W0) < dfree) || ( ((d-W0) == dfree) && (spectrum(1)(d-W0) > Cdfree) ) )
01351         return -1;
01352     }
01353 
01354 
01355   F4:
01356     if ( (W1 < dist_prof(m-1)) || (W < dist_prof(m)) ) goto F5;
01357     // select node 1
01358     if (test_catastrophic && W == W1) {
01359       c++;
01360       if (c>cat_treshold)
01361         return 0;
01362     } else {
01363       c = 0;
01364       W = W1;
01365     }
01366 
01367     S = S1;
01368     mf = 1;
01369     b++;
01370     goto F2;
01371 
01372   F5:
01373     //if (stack_pos == -1) return n_values;
01374     if (stack_pos == -1) goto end;
01375     // get node from stack
01376     S = S_stack(stack_pos);
01377     W = W_stack(stack_pos);
01378     b = b_stack(stack_pos);
01379     c = c_stack(stack_pos);
01380     stack_pos--;
01381     mf = 1;
01382     goto F2;
01383 
01384   F6:
01385     if (W0 < dist_prof(m-mf-1)) goto F4;
01386 
01387     //F7:
01388     if ( (W1 >= dist_prof(m-1)) && (W >= dist_prof(m)) ) {
01389       // save node 1
01390       if (test_catastrophic && stack_pos > 10000)
01391         return 0;
01392 
01393       stack_pos++;
01394       if (stack_pos >= max_stack_size) {
01395         max_stack_size = 2*max_stack_size;
01396         S_stack.set_size(max_stack_size, true);
01397         W_stack.set_size(max_stack_size, true);
01398         b_stack.set_size(max_stack_size, true);
01399         c_stack.set_size(max_stack_size, true);
01400       }
01401       S_stack(stack_pos) = S1;
01402       W_stack(stack_pos) = W1;
01403       b_stack(stack_pos) = b + 1; // because gone to node 1
01404       c_stack(stack_pos) = c; // because gone to node 1
01405     }
01406     // select node 0
01407     S = S0;
01408     if (test_catastrophic && W == W0) {
01409       c++;
01410       if (c>cat_treshold)
01411         return 0;
01412     } else {
01413       c = 0;
01414       W = W0;
01415     }
01416 
01417 
01418     mf++;
01419     goto F2;
01420 
01421   end:
01422     return 1;
01423   }
01424 
01425   //----------- These functions should be moved into some other place -------
01426 
01430   int reverse_int(int length, int in)
01431   {
01432     int i, j, out = 0;
01433 
01434     for (i=0; i<(length>>1); i++) {
01435       out = out | ( (in & (1<<i)) << (length-1-(i<<1)) );
01436     }
01437     for (j=0; j<(length-i); j++) {
01438       out = out | ( (in & (1<<(j+i))) >> ((j<<1)-(length&1)+1) );
01439       //out = out | ( (in & (1<<j+i)) >> ((j<<1)-(length&1)+1) ); old version with preecedence problems in MSC
01440 
01441     }
01442     return out;
01443   }
01444 
01448   int weight_int(int length, int in)
01449   {
01450     int i, w=0;
01451     for (i=0; i<length; i++) {
01452       w += (in & (1<<i)) >> i;
01453     }
01454     return w;
01455   }
01456 
01460   int compare_spectra(ivec v1, ivec v2)
01461   {
01462     it_assert_debug(v1.size() == v2.size(), "compare_spectra: wrong sizes");
01463 
01464     for (int i=0; i<v1.size(); i++) {
01465       if (v1(i) < v2(i)) {
01466         return 1;
01467       } else if (v1(i) > v2(i)) {
01468         return 0;
01469       }
01470     }
01471     return -1;
01472   }
01473 
01479   int compare_spectra(ivec v1, ivec v2, vec weight_profile)
01480   {
01481     double t1=0, t2=0;
01482     for (int i=0; i<v1.size(); i++) {
01483       t1 += (double)v1(i) * weight_profile(i);
01484       t2 += (double)v2(i) * weight_profile(i);
01485     }
01486     if (t1 < t2) return 1;
01487     else if (t1 > t2) return 0;
01488     else return -1;
01489   }
01490 
01491 } // namespace itpp
SourceForge Logo

Generated on Sun Apr 20 12:40:06 2008 for IT++ by Doxygen 1.5.5