IT++ Logo Newcom Logo

convcode.cpp

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

Generated on Wed Mar 21 12:21:41 2007 for IT++ by Doxygen 1.4.7