IT++ Logo

punct_convcode.cpp

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

Generated on Sun Sep 14 18:57:04 2008 for IT++ by Doxygen 1.5.6