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
Generated on Sat Jul 9 2011 15:21:31 for IT++ by Doxygen 1.7.4