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