IT++ Logo Newcom Logo

punct_convcode.cpp

Go to the documentation of this file.
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
SourceForge Logo

Generated on Wed Mar 21 12:21:44 2007 for IT++ by Doxygen 1.4.7