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

Generated on Sat Jul 9 2011 15:21:31 for IT++ by Doxygen 1.7.4