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
Generated on Thu Apr 19 14:19:54 2007 for IT++ by Doxygen 1.4.6