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

Generated on Sat Jul 9 2011 15:21:31 for IT++ by Doxygen 1.7.4