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 
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
SourceForge Logo

Generated on Sun Jul 26 08:54:56 2009 for IT++ by Doxygen 1.5.9