00001 00030 #include <itpp/comm/punct_convcode.h> 00031 00032 00033 namespace itpp 00034 { 00035 00036 // --------------------- Punctured_Conv_Code -------------------------------- 00037 00038 //------- Protected functions ----------------------- 00039 int Punctured_Convolutional_Code::weight(int state, int input, int time) 00040 { 00041 int i, j, temp, out, w = 0, shiftreg = state; 00042 00043 shiftreg = shiftreg | (int(input) << m); 00044 for (j = 0; j < n; j++) { 00045 if (puncture_matrix(j, time) == bin(1)) { 00046 out = 0; 00047 temp = shiftreg & gen_pol(j); 00048 for (i = 0; i < K; i++) { 00049 out ^= (temp & 1); 00050 temp = temp >> 1; 00051 } 00052 w += out; 00053 } 00054 } 00055 return w; 00056 } 00057 00058 int Punctured_Convolutional_Code::weight_reverse(int state, int input, int time) 00059 { 00060 int i, j, temp, out, w = 0, shiftreg = state; 00061 00062 shiftreg = shiftreg | (int(input) << m); 00063 for (j = 0; j < n; j++) { 00064 if (puncture_matrix(j, Period - 1 - time) == bin(1)) { 00065 out = 0; 00066 temp = shiftreg & gen_pol_rev(j); 00067 for (i = 0; i < K; i++) { 00068 out ^= (temp & 1); 00069 temp = temp >> 1; 00070 } 00071 w += out; 00072 } 00073 } 00074 return w; 00075 } 00076 00077 void Punctured_Convolutional_Code::weight(int state, int &w0, int &w1, int time) 00078 { 00079 int i, j, temp, out, shiftreg = state; 00080 w0 = 0; 00081 w1 = 0; 00082 00083 shiftreg = shiftreg | (1 << m); 00084 for (j = 0; j < n; j++) { 00085 if (puncture_matrix(j, time) == bin(1)) { 00086 out = 0; 00087 temp = shiftreg & gen_pol(j); 00088 for (i = 0; i < m; i++) { 00089 out ^= (temp & 1); 00090 temp = temp >> 1; 00091 } 00092 w0 += out; 00093 w1 += out ^(temp & 1); 00094 } 00095 } 00096 } 00097 00098 void Punctured_Convolutional_Code::weight_reverse(int state, int &w0, int &w1, int time) 00099 { 00100 int i, j, temp, out, shiftreg = state; 00101 w0 = 0; 00102 w1 = 0; 00103 00104 shiftreg = shiftreg | (1 << m); 00105 for (j = 0; j < n; j++) { 00106 if (puncture_matrix(j, Period - 1 - time) == bin(1)) { 00107 out = 0; 00108 temp = shiftreg & gen_pol_rev(j); 00109 for (i = 0; i < m; i++) { 00110 out ^= (temp & 1); 00111 temp = temp >> 1; 00112 } 00113 w0 += out; 00114 w1 += out ^(temp & 1); 00115 } 00116 } 00117 } 00118 00119 //------- Public functions ----------------------- 00120 00121 void Punctured_Convolutional_Code::set_puncture_matrix(const bmat &pmatrix) 00122 { 00123 it_error_if((pmatrix.rows() != n) || (pmatrix.cols() == 0), "Wrong size of puncture matrix"); 00124 puncture_matrix = pmatrix; 00125 Period = puncture_matrix.cols(); 00126 00127 int p, j; 00128 total = 0; 00129 00130 for (j = 0; j < n; j++) { 00131 for (p = 0; p < Period; p++) 00132 total = total + (int)(puncture_matrix(j, p)); 00133 } 00134 rate = (double)Period / total; 00135 } 00136 00137 void Punctured_Convolutional_Code::encode(const bvec &input, bvec &output) 00138 { 00139 switch (cc_method) { 00140 case Trunc: 00141 encode_trunc(input, output); 00142 break; 00143 case Tail: 00144 encode_tail(input, output); 00145 break; 00146 case Tailbite: 00147 encode_tailbite(input, output); 00148 break; 00149 default: 00150 encode_tail(input, output); 00151 break; 00152 }; 00153 } 00154 00155 void Punctured_Convolutional_Code::encode_trunc(const bvec &input, bvec &output) 00156 { 00157 Convolutional_Code::encode_trunc(input, output); 00158 00159 int nn = 0, i, p = 0, j; 00160 00161 for (i = 0; i < int(output.size() / n); i++) { 00162 for (j = 0; j < n; j++) { 00163 if (puncture_matrix(j, p) == bin(1)) { 00164 output(nn) = output(i * n + j); 00165 nn++; 00166 } 00167 } 00168 p = (p + 1) % Period; 00169 } 00170 output.set_size(nn, true); 00171 } 00172 00173 void Punctured_Convolutional_Code::encode_tail(const bvec &input, bvec &output) 00174 { 00175 Convolutional_Code::encode_tail(input, output); 00176 00177 int nn = 0, i, p = 0, j; 00178 00179 for (i = 0; i < int(output.size() / n); i++) { 00180 for (j = 0; j < n; j++) { 00181 if (puncture_matrix(j, p) == bin(1)) { 00182 output(nn) = output(i * n + j); 00183 nn++; 00184 } 00185 } 00186 p = (p + 1) % Period; 00187 } 00188 output.set_size(nn, true); 00189 } 00190 00191 void Punctured_Convolutional_Code::encode_tailbite(const bvec &input, bvec &output) 00192 { 00193 Convolutional_Code::encode_tailbite(input, output); 00194 00195 int nn = 0, i, p = 0, j; 00196 00197 for (i = 0; i < int(output.size() / n); i++) { 00198 for (j = 0; j < n; j++) { 00199 if (puncture_matrix(j, p) == bin(1)) { 00200 output(nn) = output(i * n + j); 00201 nn++; 00202 } 00203 } 00204 p = (p + 1) % Period; 00205 } 00206 output.set_size(nn, true); 00207 } 00208 00209 00210 // --------------- Hard-decision decoding is not implemented -------------------------------- 00211 void Punctured_Convolutional_Code::decode(const bvec &, bvec &) 00212 { 00213 it_error("Punctured_Convolutional_Code::decode(bvec, bvec); hard-decision decoding is not implemented"); 00214 } 00215 00216 bvec Punctured_Convolutional_Code::decode(const bvec &) 00217 { 00218 it_error("Punctured_Convolutional_Code::decode(bvec, bvec); hard-decision decoding is not implemented"); 00219 return bvec(); 00220 } 00221 00222 /* 00223 Decode a block of encoded data using specified method 00224 */ 00225 void Punctured_Convolutional_Code::decode(const vec &received_signal, bvec &output) 00226 { 00227 switch (cc_method) { 00228 case Trunc: 00229 decode_trunc(received_signal, output); 00230 break; 00231 case Tail: 00232 decode_tail(received_signal, output); 00233 break; 00234 case Tailbite: 00235 decode_tailbite(received_signal, output); 00236 break; 00237 default: 00238 decode_tail(received_signal, output); 00239 break; 00240 }; 00241 } 00242 00243 00244 // Viterbi decoder using TruncLength (=5*K if not specified) 00245 void Punctured_Convolutional_Code::decode_trunc(const vec &received_signal, bvec &output) 00246 { 00247 int nn = 0, i = 0, p = received_signal.size() / total, j; 00248 00249 int temp_size = p * Period * n; 00250 // number of punctured bits in a fraction of the puncture matrix 00251 // correcponding to the end of the received_signal vector 00252 p = received_signal.size() - p * total; 00253 // precise calculation of the number of unpunctured bits 00254 // (in the above fraction of the puncture matrix) 00255 while (p > 0) { 00256 for (j = 0; j < n; j++) { 00257 if (puncture_matrix(j, nn) == bin(1)) 00258 p--; 00259 } 00260 nn++; 00261 } 00262 temp_size += n * nn; 00263 if (p != 0) { 00264 it_warning("Punctured_Convolutional_Code::decode(): Improper length of " 00265 "the received punctured block, dummy bits have been added"); 00266 } 00267 00268 vec temp(temp_size); 00269 nn = 0; 00270 j = 0; 00271 p = 0; 00272 00273 while (nn < temp.size()) { 00274 if ((puncture_matrix(j, p) == bin(1)) && (i < received_signal.size())) { 00275 temp(nn) = received_signal(i); 00276 i++; 00277 } 00278 else { // insert dummy symbols with the same contribution for 0 and 1 00279 temp(nn) = 0; 00280 } 00281 00282 nn++; 00283 j++; 00284 00285 if (j == n) { 00286 j = 0; 00287 p = (p + 1) % Period; 00288 } 00289 } // while 00290 00291 Convolutional_Code::decode_trunc(temp, output); 00292 } 00293 00294 // Viterbi decoder using TruncLength (=5*K if not specified) 00295 void Punctured_Convolutional_Code::decode_tail(const vec &received_signal, bvec &output) 00296 { 00297 int nn = 0, i = 0, p = received_signal.size() / total, j; 00298 00299 int temp_size = p * Period * n; 00300 // number of punctured bits in a fraction of the puncture matrix 00301 // correcponding to the end of the received_signal vector 00302 p = received_signal.size() - p * total; 00303 // precise calculation of the number of unpunctured bits 00304 // (in the above fraction of the puncture matrix) 00305 while (p > 0) { 00306 for (j = 0; j < n; j++) { 00307 if (puncture_matrix(j, nn) == bin(1)) 00308 p--; 00309 } 00310 nn++; 00311 } 00312 temp_size += n * nn; 00313 if (p != 0) { 00314 it_warning("Punctured_Convolutional_Code::decode_tail(): Improper length " 00315 "of the received punctured block, dummy bits have been added"); 00316 } 00317 00318 vec temp(temp_size); 00319 00320 nn = 0; 00321 j = 0; 00322 p = 0; 00323 00324 while (nn < temp.size()) { 00325 if ((puncture_matrix(j, p) == bin(1)) && (i < received_signal.size())) { 00326 temp(nn) = received_signal(i); 00327 i++; 00328 } 00329 else { // insert dummy symbols with same contribution for 0 and 1 00330 temp(nn) = 0; 00331 } 00332 00333 nn++; 00334 j++; 00335 00336 if (j == n) { 00337 j = 0; 00338 p = (p + 1) % Period; 00339 } 00340 } // while 00341 00342 Convolutional_Code::decode_tail(temp, output); 00343 } 00344 00345 // Decode a block of encoded data where encode_tailbite has been used. Tries all start states. 00346 void Punctured_Convolutional_Code::decode_tailbite(const vec &received_signal, bvec &output) 00347 { 00348 int nn = 0, i = 0, p = received_signal.size() / total, j; 00349 00350 int temp_size = p * Period * n; 00351 // number of punctured bits in a fraction of the puncture matrix 00352 // correcponding to the end of the received_signal vector 00353 p = received_signal.size() - p * total; 00354 // precise calculation of the number of unpunctured bits 00355 // (in the above fraction of the puncture matrix) 00356 while (p > 0) { 00357 for (j = 0; j < n; j++) { 00358 if (puncture_matrix(j, nn) == bin(1)) 00359 p--; 00360 } 00361 nn++; 00362 } 00363 temp_size += n * nn; 00364 if (p != 0) { 00365 it_warning("Punctured_Convolutional_Code::decode_tailbite(): Improper " 00366 "length of the received punctured block, dummy bits have been " 00367 "added"); 00368 } 00369 00370 vec temp(temp_size); 00371 00372 nn = 0; 00373 j = 0; 00374 p = 0; 00375 00376 while (nn < temp.size()) { 00377 if ((puncture_matrix(j, p) == bin(1)) && (i < received_signal.size())) { 00378 temp(nn) = received_signal(i); 00379 i++; 00380 } 00381 else { // insert dummy symbols with same contribution for 0 and 1 00382 temp(nn) = 0; 00383 } 00384 00385 nn++; 00386 j++; 00387 00388 if (j == n) { 00389 j = 0; 00390 p = (p + 1) % Period; 00391 } 00392 } // while 00393 00394 Convolutional_Code::decode_tailbite(temp, output); 00395 } 00396 00397 00398 /* 00399 Calculate the inverse sequence 00400 00401 Assumes that encode_tail is used in the encoding process. Returns false if there is an 00402 error in the coded sequence (not a valid codeword). Does not check that the tail forces 00403 the encoder into the zeroth state. 00404 */ 00405 bool Punctured_Convolutional_Code::inverse_tail(const bvec coded_sequence, bvec &input) 00406 { 00407 int state = 0, zero_state, one_state, zero_temp, one_temp, i, j, p = 0, prev_pos = 0, no_symbols; 00408 int block_length = 0; 00409 bvec zero_output(n), one_output(n), temp_output(n); 00410 00411 block_length = coded_sequence.size() * Period / total - m; 00412 00413 it_error_if(block_length <= 0, "The input sequence is to short"); 00414 input.set_length(block_length, false); // not include the tail in the output 00415 00416 p = 0; 00417 00418 for (i = 0; i < block_length; i++) { 00419 zero_state = state; 00420 one_state = state | (1 << m); 00421 no_symbols = 0; 00422 for (j = 0; j < n; j++) { 00423 if (puncture_matrix(j, p) == bin(1)) { 00424 zero_temp = zero_state & gen_pol(j); 00425 one_temp = one_state & gen_pol(j); 00426 zero_output(no_symbols) = xor_int_table(zero_temp); 00427 one_output(no_symbols) = xor_int_table(one_temp); 00428 no_symbols++; 00429 } 00430 } 00431 if (coded_sequence.mid(prev_pos, no_symbols) == zero_output.left(no_symbols)) { 00432 input(i) = bin(0); 00433 state = zero_state >> 1; 00434 } 00435 else if (coded_sequence.mid(prev_pos, no_symbols) == one_output.left(no_symbols)) { 00436 input(i) = bin(1); 00437 state = one_state >> 1; 00438 } 00439 else 00440 return false; 00441 00442 prev_pos += no_symbols; 00443 p = (p + 1) % Period; 00444 } 00445 00446 return true; 00447 } 00448 00449 00450 00451 00452 /* 00453 It is possible to do this more efficiently by registrating all (inputs,states,Periods) 00454 that has zero metric (just and with the generators). After that build paths from a input=1 state. 00455 */ 00456 bool Punctured_Convolutional_Code::catastrophic(void) 00457 { 00458 int max_stack_size = 50000; 00459 ivec S_stack(max_stack_size), t_stack(max_stack_size); 00460 int start, S, W0, W1, S0, S1, pos, t = 0; 00461 int stack_pos = -1; 00462 00463 for (pos = 0; pos < Period; pos++) { // test all start positions 00464 for (start = 0; start < (1 << m); start++) { 00465 stack_pos = -1; 00466 S = start; 00467 t = 0; 00468 00469 node1: 00470 if (t > (1 << m) * Period) { return true; } 00471 S0 = next_state(S, 0); 00472 S1 = next_state(S, 1); 00473 weight(S, W0, W1, (pos + t) % Period); 00474 if (W1 > 0) { goto node0; } 00475 if ((W0 == 0) && (S0 == start) && (((pos + t + 1) % Period) == pos)) { return true; } 00476 if ((W0 == 0) && (S0 == 0) && (S != 0)) { return true; } 00477 if ((S0 != 0) && (W0 == 0)) { 00478 stack_pos++; 00479 if (stack_pos >= max_stack_size) { 00480 max_stack_size = 2 * max_stack_size; 00481 S_stack.set_size(max_stack_size, true); 00482 t_stack.set_size(max_stack_size, true); 00483 } 00484 S_stack(stack_pos) = S0; 00485 t_stack(stack_pos) = t + 1; 00486 } 00487 if ((W1 == 0) && (S1 == start) && (((pos + t + 1) % Period) == pos)) { return true; } 00488 S = S1; 00489 t++; 00490 goto node1; 00491 00492 node0: 00493 if (W0 > 0) goto stack; 00494 if ((W0 == 0) && (S0 == start) && (((pos + t + 1) % Period) == pos)) { return true; } 00495 if ((W0 == 0) && (S0 == 0) && (S != 0)) { return true; } 00496 if (S0 != 0) { 00497 S = S0; 00498 t++; 00499 goto node1; 00500 } 00501 else { 00502 goto stack; 00503 } 00504 00505 stack: 00506 if (stack_pos >= 0) { 00507 S = S_stack(stack_pos); 00508 t = t_stack(stack_pos); 00509 stack_pos--; 00510 goto node1; 00511 } 00512 } 00513 } 00514 return false; 00515 } 00516 00517 void Punctured_Convolutional_Code::distance_profile(ivec &dist_prof, int start_time, int dmax, bool reverse) 00518 { 00519 int max_stack_size = 50000; 00520 ivec S_stack(max_stack_size), W_stack(max_stack_size), t_stack(max_stack_size); 00521 00522 int stack_pos = -1, t, S, W, W0, w0, w1; 00523 00524 00525 dist_prof.zeros(); 00526 dist_prof += dmax; // just select a big number! 00527 if (reverse) 00528 W = weight_reverse(0, 1, start_time); 00529 else 00530 W = weight(0, 1, start_time); 00531 00532 S = next_state(0, 1); // start from zero state and a one input 00533 t = 0; 00534 dist_prof(0) = W; 00535 00536 node1: 00537 if (reverse) 00538 weight_reverse(S, w0, w1, (start_time + t + 1) % Period); 00539 else 00540 weight(S, w0, w1, (start_time + t + 1) % Period); 00541 00542 if (t < m) { 00543 W0 = W + w0; 00544 if (W0 < dist_prof(m)) { // store node0 (S, W0, and t+1) 00545 stack_pos++; 00546 if (stack_pos >= max_stack_size) { 00547 max_stack_size = 2 * max_stack_size; 00548 S_stack.set_size(max_stack_size, true); 00549 W_stack.set_size(max_stack_size, true); 00550 t_stack.set_size(max_stack_size, true); 00551 } 00552 S_stack(stack_pos) = next_state(S, 0); 00553 W_stack(stack_pos) = W0; 00554 t_stack(stack_pos) = t + 1; 00555 } 00556 } 00557 else goto stack; 00558 00559 W += w1; 00560 if (W > dist_prof(m)) 00561 goto stack; 00562 00563 t++; 00564 S = next_state(S, 1); 00565 00566 if (W < dist_prof(t)) 00567 dist_prof(t) = W; 00568 00569 if (t == m) goto stack; 00570 goto node1; 00571 00572 00573 stack: 00574 if (stack_pos >= 0) { 00575 // get root node from stack 00576 S = S_stack(stack_pos); 00577 W = W_stack(stack_pos); 00578 t = t_stack(stack_pos); 00579 stack_pos--; 00580 00581 if (W < dist_prof(t)) 00582 dist_prof(t) = W; 00583 00584 if (t == m) goto stack; 00585 00586 goto node1; 00587 } 00588 } 00589 00590 int Punctured_Convolutional_Code::fast(Array<ivec> &spectrum, int start_time, int dfree, int no_terms, 00591 int d_best_so_far, bool test_catastrophic) 00592 { 00593 int cat_treshold = (1 << m) * Period; 00594 int i, j, t = 0; 00595 ivec dist_prof(K), dist_prof_rev(K), dist_prof_temp(K); 00596 00597 //calculate distance profile 00598 distance_profile(dist_prof, start_time, dfree); 00599 distance_profile(dist_prof_rev, 0, dfree, true); // for the reverse code 00600 00601 int dist_prof_rev0 = dist_prof_rev(0); 00602 00603 bool reverse = true; // true = use reverse dist_prof 00604 00605 // choose the lowest dist_prof over all periods 00606 for (i = 1; i < Period; i++) { 00607 distance_profile(dist_prof_temp, i, dfree, true); 00608 // switch if new is lower 00609 for (j = 0; j < K; j++) { 00610 if (dist_prof_temp(j) < dist_prof_rev(j)) { 00611 dist_prof_rev(j) = dist_prof_temp(j); 00612 } 00613 } 00614 } 00615 00616 dist_prof_rev0 = dist_prof(0); 00617 dist_prof = dist_prof_rev; 00618 00619 int d = dfree + no_terms - 1; 00620 int max_stack_size = 50000; 00621 ivec S_stack(max_stack_size), W_stack(max_stack_size), b_stack(max_stack_size); 00622 ivec c_stack(max_stack_size), t_stack(max_stack_size); 00623 int stack_pos = -1; 00624 00625 // F1. 00626 int mf = 1, b = 1; 00627 spectrum.set_size(2); 00628 spectrum(0).set_size(dfree + no_terms, 0); // ad 00629 spectrum(1).set_size(dfree + no_terms, 0); // cd 00630 spectrum(0).zeros(); 00631 spectrum(1).zeros(); 00632 int S, S0, S1, W0, W1, w0, w1, c = 0; 00633 S = next_state(0, 1); // start in zero state and one input 00634 int W = d - dist_prof_rev0; 00635 t = 1; 00636 00637 F2: 00638 S0 = next_state(S, 0); 00639 S1 = next_state(S, 1); 00640 00641 if (reverse) { 00642 weight(S, w0, w1, (start_time + t) % Period); 00643 } 00644 else { 00645 weight_reverse(S, w0, w1, (start_time + t) % Period); 00646 } 00647 W0 = W - w0; 00648 W1 = W - w1; 00649 00650 if (mf < m) goto F6; 00651 00652 //F3: 00653 if (W0 >= 0) { 00654 spectrum(0)(d - W0)++; 00655 spectrum(1)(d - W0) += b; 00656 } 00657 //Test on d_best_so_far 00658 if ((d - W0) < d_best_so_far) { return -1; } 00659 00660 F4: 00661 if ((W1 < dist_prof(m - 1)) || (W < dist_prof(m))) goto F5; 00662 // select node 1 00663 if (test_catastrophic && (W == W1)) { 00664 c++; 00665 if (c > cat_treshold) 00666 return 0; 00667 } 00668 else { 00669 c = 0; 00670 } 00671 00672 W = W1; 00673 S = S1; 00674 t++; 00675 mf = 1; 00676 b++; 00677 goto F2; 00678 00679 F5: 00680 if (stack_pos == -1) goto end; 00681 // get node from stack 00682 S = S_stack(stack_pos); 00683 W = W_stack(stack_pos); 00684 b = b_stack(stack_pos); 00685 c = c_stack(stack_pos); 00686 t = t_stack(stack_pos); 00687 stack_pos--; 00688 mf = 1; 00689 goto F2; 00690 00691 F6: 00692 if (W0 < dist_prof(m - mf - 1)) goto F4; 00693 00694 //F7: 00695 if ((W1 >= dist_prof(m - 1)) && (W >= dist_prof(m))) { 00696 // save node 1 00697 if (test_catastrophic && (stack_pos > 40000)) 00698 return 0; 00699 00700 stack_pos++; 00701 if (stack_pos >= max_stack_size) { 00702 max_stack_size = 2 * max_stack_size; 00703 S_stack.set_size(max_stack_size, true); 00704 W_stack.set_size(max_stack_size, true); 00705 b_stack.set_size(max_stack_size, true); 00706 c_stack.set_size(max_stack_size, true); 00707 t_stack.set_size(max_stack_size, true); 00708 } 00709 S_stack(stack_pos) = S1; 00710 W_stack(stack_pos) = W1; 00711 b_stack(stack_pos) = b + 1; // because gone to node 1 00712 c_stack(stack_pos) = c; 00713 t_stack(stack_pos) = t + 1; 00714 } 00715 // select node 0 00716 S = S0; 00717 t++; 00718 if (test_catastrophic && (W == W0)) { 00719 c++; 00720 if (c > cat_treshold) 00721 return false; 00722 } 00723 else { 00724 c = 0; 00725 } 00726 00727 W = W0; 00728 mf++; 00729 goto F2; 00730 00731 end: 00732 return 1; 00733 } 00734 00735 void Punctured_Convolutional_Code::calculate_spectrum(Array<ivec> &spectrum, int dmax, int no_terms) 00736 { 00737 Array<ivec> temp_spectra(2); 00738 spectrum.set_size(2); 00739 spectrum(0).set_size(dmax + no_terms, false); 00740 spectrum(1).set_size(dmax + no_terms, false); 00741 spectrum(0).zeros(); 00742 spectrum(1).zeros(); 00743 00744 for (int pos = 0; pos < Period; pos++) { 00745 calculate_spectrum(temp_spectra, pos, dmax, no_terms); 00746 spectrum(0) += temp_spectra(0); 00747 spectrum(1) += temp_spectra(1); 00748 } 00749 } 00750 00751 void Punctured_Convolutional_Code::calculate_spectrum(Array<ivec> &spectrum, int time, int dmax, int no_terms, int block_length) 00752 { 00753 imat Ad_states(1 << (K - 1), dmax + no_terms), Cd_states(1 << m, dmax + no_terms); 00754 imat Ad_temp(1 << (K - 1), dmax + no_terms), Cd_temp(1 << m, dmax + no_terms); 00755 ivec mindist(1 << (K - 1)), mindist_temp(1 << m); 00756 00757 spectrum.set_size(2); 00758 spectrum(0).set_size(dmax + no_terms, false); 00759 spectrum(1).set_size(dmax + no_terms, false); 00760 spectrum(0).zeros(); 00761 spectrum(1).zeros(); 00762 Ad_states.zeros(); 00763 Cd_states.zeros(); 00764 mindist.zeros(); 00765 int wmax = dmax + no_terms; 00766 ivec visited_states(1 << m), visited_states_temp(1 << m); 00767 bool proceede, expand_s1; 00768 int t, d, w0, w1, s, s0, s1 = 0, s_start; 00769 00770 // initial phase. Evolv trellis up to time K. 00771 visited_states.zeros(); // 0 = false 00772 00773 s1 = next_state(0, 1); // Start in state 0 and 1 input 00774 visited_states(s1) = 1; // 1 = true. 00775 w1 = weight(0, 1, time); 00776 Ad_states(s1, w1) = 1; 00777 Cd_states(s1, w1) = 1; 00778 mindist(s1) = w1; 00779 00780 if (block_length > 0) { 00781 s0 = next_state(0, 0); 00782 visited_states(s0) = 1; // 1 = true. Expand also the zero-path 00783 w0 = weight(0, 0, time); 00784 Ad_states(s0, w0) = 1; 00785 Cd_states(s0, w0) = 0; 00786 mindist(s0) = w0; 00787 s_start = 0; 00788 } 00789 else { 00790 s_start = 1; 00791 } 00792 00793 t = 1; 00794 proceede = true; 00795 while (proceede) { 00796 proceede = false; 00797 Ad_temp.zeros(); 00798 Cd_temp.zeros(); 00799 mindist_temp.zeros(); 00800 visited_states_temp.zeros(); //false 00801 00802 for (s = s_start; s < (1 << m); s++) { 00803 00804 if (visited_states(s) == 1) { 00805 if ((mindist(s) >= 0) && (mindist(s) < wmax)) { 00806 proceede = true; 00807 weight(s, w0, w1, (time + t) % Period); 00808 00809 s0 = next_state(s, 0); 00810 for (d = mindist(s); d < (wmax - w0); d++) { 00811 Ad_temp(s0, d + w0) += Ad_states(s, d); 00812 Cd_temp(s0, d + w0) += Cd_states(s, d); 00813 } 00814 if (visited_states_temp(s0) == 0) { mindist_temp(s0) = mindist(s) + w0; } 00815 else { mindist_temp(s0) = ((mindist(s) + w0) < mindist_temp(s0)) ? mindist(s) + w0 : mindist_temp(s0); } 00816 visited_states_temp(s0) = 1; //true 00817 00818 expand_s1 = false; 00819 if ((block_length == 0) || (t < (block_length - (K - 1)))) { 00820 expand_s1 = true; 00821 } 00822 00823 if (expand_s1) { 00824 s1 = next_state(s, 1); 00825 for (d = mindist(s); d < (wmax - w1); d++) { 00826 Ad_temp(s1, d + w1) += Ad_states(s, d); 00827 Cd_temp(s1, d + w1) += Cd_states(s, d) + Ad_states(s, d); 00828 } 00829 if (visited_states_temp(s1) == 0) { mindist_temp(s1) = mindist(s) + w1; } 00830 else { mindist_temp(s1) = ((mindist(s) + w1) < mindist_temp(s1)) ? mindist(s) + w1 : mindist_temp(s1); } 00831 visited_states_temp(s1) = 1; //true 00832 } 00833 } 00834 } 00835 } 00836 00837 Ad_states = Ad_temp; 00838 Cd_states = Cd_temp; 00839 if (block_length == 0) { 00840 spectrum(0) += Ad_temp.get_row(0); 00841 spectrum(1) += Cd_temp.get_row(0); 00842 } 00843 visited_states = visited_states_temp; 00844 mindist = elem_mult(mindist_temp, visited_states); 00845 t++; 00846 if ((t > block_length) && (block_length > 0)) { proceede = false; } 00847 } 00848 00849 if (block_length > 0) { 00850 spectrum(0) = Ad_states.get_row(0); 00851 spectrum(1) = Cd_states.get_row(0); 00852 } 00853 00854 } 00855 00856 } // namespace itpp
Generated on Tue Nov 23 08:47:56 2010 for IT++ by Doxygen 1.6.1