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