IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
convcode.cpp
Go to the documentation of this file.
1 
29 #include <itpp/comm/convcode.h>
30 #include <itpp/base/binary.h>
31 #include <itpp/base/matfunc.h>
32 #include <limits>
33 
34 namespace itpp
35 {
36 
37 // ----------------- Protected functions -----------------------------
38 
39 /*
40  The weight of the transition from given state with the input given
41 */
42 int Convolutional_Code::weight(const int state, const int input)
43 {
44  int i, j, temp, out, w = 0, shiftreg = state;
45 
46  shiftreg = shiftreg | (input << m);
47  for (j = 0; j < n; j++) {
48  out = 0;
49  temp = shiftreg & gen_pol(j);
50  for (i = 0; i < K; i++) {
51  out ^= (temp & 1);
52  temp = temp >> 1;
53  }
54  w += out;
55  //w += weight_int_table(temp);
56  }
57  return w;
58 }
59 
60 /*
61  The weight (of the reverse code) of the transition from given state with
62  the input given
63 */
64 int Convolutional_Code::weight_reverse(const int state, const int input)
65 {
66  int i, j, temp, out, w = 0, shiftreg = state;
67 
68  shiftreg = shiftreg | (input << m);
69  for (j = 0; j < n; j++) {
70  out = 0;
71  temp = shiftreg & gen_pol_rev(j);
72  for (i = 0; i < K; i++) {
73  out ^= (temp & 1);
74  temp = temp >> 1;
75  }
76  w += out;
77  //w += weight_int_table(temp);
78  }
79  return w;
80 }
81 
82 /*
83  The weight of the two paths (input 0 or 1) from given state
84 */
85 void Convolutional_Code::weight(const int state, int &w0, int &w1)
86 {
87  int i, j, temp, out, shiftreg = state;
88  w0 = 0;
89  w1 = 0;
90 
91  shiftreg = shiftreg | (1 << m);
92  for (j = 0; j < n; j++) {
93  out = 0;
94  temp = shiftreg & gen_pol(j);
95 
96  for (i = 0; i < m; i++) {
97  out ^= (temp & 1);
98  temp = temp >> 1;
99  }
100  w0 += out;
101  w1 += out ^(temp & 1);
102  }
103 }
104 
105 /*
106  The weight (of the reverse code) of the two paths (input 0 or 1) from
107  given state
108 */
109 void Convolutional_Code::weight_reverse(const int state, int &w0, int &w1)
110 {
111  int i, j, temp, out, shiftreg = state;
112  w0 = 0;
113  w1 = 0;
114 
115  shiftreg = shiftreg | (1 << m);
116  for (j = 0; j < n; j++) {
117  out = 0;
118  temp = shiftreg & gen_pol_rev(j);
119 
120  for (i = 0; i < m; i++) {
121  out ^= (temp & 1);
122  temp = temp >> 1;
123  }
124  w0 += out;
125  w1 += out ^(temp & 1);
126  }
127 }
128 
129 /*
130  Output on transition (backwards) with input from state
131 */
132 bvec Convolutional_Code::output_reverse(const int state, const int input)
133 {
134  int j, temp, temp_state;
135 
136  bvec output(n);
137 
138  temp_state = (state << 1) | input;
139  for (j = 0; j < n; j++) {
140  temp = temp_state & gen_pol(j);
141  output(j) = xor_int_table(temp);
142  }
143 
144  return output;
145 }
146 
147 /*
148  Output on transition (backwards) with input from state
149 */
150 void Convolutional_Code::output_reverse(const int state, bvec &zero_output,
151  bvec &one_output)
152 {
153  int j, temp, temp_state;
154  bin one_bit;
155 
156  temp_state = (state << 1) | 1;
157  for (j = 0; j < n; j++) {
158  temp = temp_state & gen_pol(j);
159  one_bit = temp & 1;
160  temp = temp >> 1;
161  one_output(j) = xor_int_table(temp) ^ one_bit;
162  zero_output(j) = xor_int_table(temp);
163  }
164 }
165 
166 /*
167  Output on transition (backwards) with input from state
168 */
169 void Convolutional_Code::output_reverse(const int state, int &zero_output,
170  int &one_output)
171 {
172  int j, temp, temp_state;
173  bin one_bit;
174 
175  zero_output = 0, one_output = 0;
176  temp_state = (state << 1) | 1;
177  for (j = 0; j < n; j++) {
178  temp = temp_state & gen_pol(j);
179  one_bit = temp & 1;
180  temp = temp >> 1;
181 
182  one_output = (one_output << 1) | int(xor_int_table(temp) ^ one_bit);
183  zero_output = (zero_output << 1) | int(xor_int_table(temp));
184  }
185 }
186 
187 /*
188  Output on transition (backwards) with input from state
189 */
191  const vec &rx_codeword,
192  double &zero_metric,
193  double &one_metric)
194 {
195  int temp;
196  bin one_bit;
197  one_metric = zero_metric = 0;
198 
199  int temp_state = (state << 1) | 1;
200  for (int j = 0; j < n; j++) {
201  temp = temp_state & gen_pol(j);
202  one_bit = temp & 1;
203  temp >>= 1;
204  one_metric += (2 * static_cast<int>(xor_int_table(temp) ^ one_bit) - 1)
205  * rx_codeword(j);
206  zero_metric += (2 * static_cast<int>(xor_int_table(temp)) - 1)
207  * rx_codeword(j);
208  }
209 }
210 
211 
212 // calculates metrics for all codewords (2^n of them) in natural order
213 void Convolutional_Code::calc_metric(const vec &rx_codeword,
214  vec &delta_metrics)
215 {
216  int no_outputs = pow2i(n), no_loop = pow2i(n - 1), mask = no_outputs - 1,
217  temp;
218  delta_metrics.set_size(no_outputs, false);
219 
220  if (no_outputs <= no_states) {
221  for (int i = 0; i < no_loop; i++) { // all input possibilities
222  delta_metrics(i) = 0;
223  temp = i;
224  for (int j = n - 1; j >= 0; j--) {
225  if (temp & 1)
226  delta_metrics(i) += rx_codeword(j);
227  else
228  delta_metrics(i) -= rx_codeword(j);
229  temp >>= 1;
230  }
231  delta_metrics(i ^ mask) = -delta_metrics(i); // the inverse codeword
232  }
233  }
234  else {
235  double zero_metric, one_metric;
236  int zero_output, one_output, temp_state;
237  bin one_bit;
238  for (int s = 0; s < no_states; s++) { // all states
239  zero_metric = 0, one_metric = 0;
240  zero_output = 0, one_output = 0;
241  temp_state = (s << 1) | 1;
242  for (int j = 0; j < n; j++) {
243  temp = temp_state & gen_pol(j);
244  one_bit = temp & 1;
245  temp >>= 1;
246  if (xor_int_table(temp)) {
247  zero_metric += rx_codeword(j);
248  one_metric -= rx_codeword(j);
249  }
250  else {
251  zero_metric -= rx_codeword(j);
252  one_metric += rx_codeword(j);
253  }
254  one_output = (one_output << 1)
255  | static_cast<int>(xor_int_table(temp) ^ one_bit);
256  zero_output = (zero_output << 1)
257  | static_cast<int>(xor_int_table(temp));
258  }
259  delta_metrics(zero_output) = zero_metric;
260  delta_metrics(one_output) = one_metric;
261  }
262  }
263 }
264 
266 
267 // MFD codes R=1/2
268 int Conv_Code_MFD_2[15][2] = {
269  {0, 0},
270  {0, 0},
271  {0, 0},
272  {05, 07},
273  {015, 017},
274  {023, 035},
275  {053, 075},
276  {0133, 0171},
277  {0247, 0371},
278  {0561, 0753},
279  {01167, 01545},
280  {02335, 03661},
281  {04335, 05723},
282  {010533, 017661},
283  {021675, 027123},
284 };
285 
286 // MFD codes R=1/3
287 int Conv_Code_MFD_3[15][3] = {
288  {0, 0, 0},
289  {0, 0, 0},
290  {0, 0, 0},
291  {05, 07, 07},
292  {013, 015, 017},
293  {025, 033, 037},
294  {047, 053, 075},
295  {0133, 0145, 0175},
296  {0225, 0331, 0367},
297  {0557, 0663, 0711},
298  {0117, 01365, 01633},
299  {02353, 02671, 03175},
300  {04767, 05723, 06265},
301  {010533, 010675, 017661},
302  {021645, 035661, 037133},
303 };
304 
305 // MFD codes R=1/4
306 int Conv_Code_MFD_4[15][4] = {
307  {0, 0, 0, 0},
308  {0, 0, 0, 0},
309  {0, 0, 0, 0},
310  {05, 07, 07, 07},
311  {013, 015, 015, 017},
312  {025, 027, 033, 037},
313  {053, 067, 071, 075},
314  {0135, 0135, 0147, 0163},
315  {0235, 0275, 0313, 0357},
316  {0463, 0535, 0733, 0745},
317  {0117, 01365, 01633, 01653},
318  {02327, 02353, 02671, 03175},
319  {04767, 05723, 06265, 07455},
320  {011145, 012477, 015573, 016727},
321  {021113, 023175, 035527, 035537},
322 };
323 
324 
325 // MFD codes R=1/5
326 int Conv_Code_MFD_5[9][5] = {
327  {0, 0, 0, 0, 0},
328  {0, 0, 0, 0, 0},
329  {0, 0, 0, 0, 0},
330  {07, 07, 07, 05, 05},
331  {017, 017, 013, 015, 015},
332  {037, 027, 033, 025, 035},
333  {075, 071, 073, 065, 057},
334  {0175, 0131, 0135, 0135, 0147},
335  {0257, 0233, 0323, 0271, 0357},
336 };
337 
338 // MFD codes R=1/6
339 int Conv_Code_MFD_6[9][6] = {
340  {0, 0, 0, 0, 0, 0},
341  {0, 0, 0, 0, 0, 0},
342  {0, 0, 0, 0, 0, 0},
343  {07, 07, 07, 07, 05, 05},
344  {017, 017, 013, 013, 015, 015},
345  {037, 035, 027, 033, 025, 035},
346  {073, 075, 055, 065, 047, 057},
347  {0173, 0151, 0135, 0135, 0163, 0137},
348  {0253, 0375, 0331, 0235, 0313, 0357},
349 };
350 
351 // MFD codes R=1/7
352 int Conv_Code_MFD_7[9][7] = {
353  {0, 0, 0, 0, 0, 0, 0},
354  {0, 0, 0, 0, 0, 0, 0},
355  {0, 0, 0, 0, 0, 0, 0},
356  {07, 07, 07, 07, 05, 05, 05},
357  {017, 017, 013, 013, 013, 015, 015},
358  {035, 027, 025, 027, 033, 035, 037},
359  {053, 075, 065, 075, 047, 067, 057},
360  {0165, 0145, 0173, 0135, 0135, 0147, 0137},
361  {0275, 0253, 0375, 0331, 0235, 0313, 0357},
362 };
363 
364 // MFD codes R=1/8
365 int Conv_Code_MFD_8[9][8] = {
366  {0, 0, 0, 0, 0, 0, 0, 0},
367  {0, 0, 0, 0, 0, 0, 0, 0},
368  {0, 0, 0, 0, 0, 0, 0, 0},
369  {07, 07, 05, 05, 05, 07, 07, 07},
370  {017, 017, 013, 013, 013, 015, 015, 017},
371  {037, 033, 025, 025, 035, 033, 027, 037},
372  {057, 073, 051, 065, 075, 047, 067, 057},
373  {0153, 0111, 0165, 0173, 0135, 0135, 0147, 0137},
374  {0275, 0275, 0253, 0371, 0331, 0235, 0313, 0357},
375 };
376 
377 int maxK_Conv_Code_MFD[9] = {0, 0, 14, 14, 14, 8, 8, 8, 8};
378 
379 void get_MFD_gen_pol(int n, int K, ivec & gen)
380 {
381  gen.set_size(n);
382 
383  switch (n) {
384  case 2: // Rate 1/2
385  it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[2], "This convolutional code doesn't exist in the tables");
386  gen(0) = Conv_Code_MFD_2[K][0];
387  gen(1) = Conv_Code_MFD_2[K][1];
388  break;
389  case 3: // Rate 1/3
390  it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[3], "This convolutional code doesn't exist in the tables");
391  gen(0) = Conv_Code_MFD_3[K][0];
392  gen(1) = Conv_Code_MFD_3[K][1];
393  gen(2) = Conv_Code_MFD_3[K][2];
394  break;
395  case 4: // Rate 1/4
396  it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[4], "This convolutional code doesn't exist in the tables");
397  gen(0) = Conv_Code_MFD_4[K][0];
398  gen(1) = Conv_Code_MFD_4[K][1];
399  gen(2) = Conv_Code_MFD_4[K][2];
400  gen(3) = Conv_Code_MFD_4[K][3];
401  break;
402  case 5: // Rate 1/5
403  it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[5], "This convolutional code doesn't exist in the tables");
404  gen(0) = Conv_Code_MFD_5[K][0];
405  gen(1) = Conv_Code_MFD_5[K][1];
406  gen(2) = Conv_Code_MFD_5[K][2];
407  gen(3) = Conv_Code_MFD_5[K][3];
408  gen(4) = Conv_Code_MFD_5[K][4];
409  break;
410  case 6: // Rate 1/6
411  it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[6], "This convolutional code doesn't exist in the tables");
412  gen(0) = Conv_Code_MFD_6[K][0];
413  gen(1) = Conv_Code_MFD_6[K][1];
414  gen(2) = Conv_Code_MFD_6[K][2];
415  gen(3) = Conv_Code_MFD_6[K][3];
416  gen(4) = Conv_Code_MFD_6[K][4];
417  gen(5) = Conv_Code_MFD_6[K][5];
418  break;
419  case 7: // Rate 1/7
420  it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[7], "This convolutional code doesn't exist in the tables");
421  gen(0) = Conv_Code_MFD_7[K][0];
422  gen(1) = Conv_Code_MFD_7[K][1];
423  gen(2) = Conv_Code_MFD_7[K][2];
424  gen(3) = Conv_Code_MFD_7[K][3];
425  gen(4) = Conv_Code_MFD_7[K][4];
426  gen(5) = Conv_Code_MFD_7[K][5];
427  gen(6) = Conv_Code_MFD_7[K][6];
428  break;
429  case 8: // Rate 1/8
430  it_assert(K >= 3 && K <= maxK_Conv_Code_MFD[8], "This convolutional code doesn't exist in the tables");
431  gen(0) = Conv_Code_MFD_8[K][0];
432  gen(1) = Conv_Code_MFD_8[K][1];
433  gen(2) = Conv_Code_MFD_8[K][2];
434  gen(3) = Conv_Code_MFD_8[K][3];
435  gen(4) = Conv_Code_MFD_8[K][4];
436  gen(5) = Conv_Code_MFD_8[K][5];
437  gen(6) = Conv_Code_MFD_8[K][6];
438  gen(7) = Conv_Code_MFD_8[K][7];
439  break;
440  default: // wrong rate
441  it_assert(false, "This convolutional code doesn't exist in the tables");
442  }
443 }
444 
445 // ODS codes R=1/2
446 int Conv_Code_ODS_2[17][2] = {
447  {0, 0},
448  {0, 0},
449  {0, 0},
450  {05, 07},
451  {015, 017},
452  {023, 035},
453  {053, 075},
454  {0133, 0171},
455  {0247, 0371},
456  {0561, 0753},
457  {01151, 01753},
458  {03345, 03613},
459  {05261, 07173},
460  {012767, 016461},
461  {027251, 037363},
462  {063057, 044735},
463  {0126723, 0152711},
464 };
465 
466 // ODS codes R=1/3
467 int Conv_Code_ODS_3[14][3] = {
468  {0, 0, 0},
469  {0, 0, 0},
470  {0, 0, 0},
471  {05, 07, 07},
472  {013, 015, 017},
473  {025, 033, 037},
474  {047, 053, 075},
475  {0133, 0165, 0171},
476  {0225, 0331, 0367},
477  {0575, 0623, 0727},
478  {01233, 01375, 01671},
479  {02335, 02531, 03477},
480  {05745, 06471, 07553},
481  {013261, 015167, 017451},
482 };
483 
484 // ODS codes R=1/4
485 int Conv_Code_ODS_4[12][4] = {
486  {0, 0, 0, 0},
487  {0, 0, 0, 0},
488  {0, 0, 0, 0},
489  {05, 05, 07, 07},
490  {013, 015, 015, 017},
491  {025, 027, 033, 037},
492  {051, 055, 067, 077},
493  {0117, 0127, 0155, 0171},
494  {0231, 0273, 0327, 0375},
495  {0473, 0513, 0671, 0765},
496  {01173, 01325, 01467, 01751},
497  {02565, 02747, 03311, 03723},
498 };
499 
500 int maxK_Conv_Code_ODS[5] = {0, 0, 16, 13, 11};
501 
502 void get_ODS_gen_pol(int n, int K, ivec & gen)
503 {
504  gen.set_size(n);
505 
506  switch (n) {
507  case 2: // Rate 1/2
508  it_assert(K >= 3 && K <= maxK_Conv_Code_ODS[2], "This convolutional code doesn't exist in the tables");
509  gen(0) = Conv_Code_ODS_2[K][0];
510  gen(1) = Conv_Code_ODS_2[K][1];
511  break;
512  case 3: // Rate 1/3
513  it_assert(K >= 3 && K <= maxK_Conv_Code_ODS[3], "This convolutional code doesn't exist in the tables");
514  gen(0) = Conv_Code_ODS_3[K][0];
515  gen(1) = Conv_Code_ODS_3[K][1];
516  gen(2) = Conv_Code_ODS_3[K][2];
517  break;
518  case 4: // Rate 1/4
519  it_assert(K >= 3 && K <= maxK_Conv_Code_ODS[4], "This convolutional code doesn't exist in the tables");
520  gen(0) = Conv_Code_ODS_4[K][0];
521  gen(1) = Conv_Code_ODS_4[K][1];
522  gen(2) = Conv_Code_ODS_4[K][2];
523  gen(3) = Conv_Code_ODS_4[K][3];
524  break;
525  default: // wrong rate
526  it_assert(false, "This convolutional code doesn't exist in the tables");
527  }
528 }
529 
531 
532 // --------------- Public functions -------------------------
533 
535  int inverse_rate, int constraint_length)
536 {
537  ivec gen;
538 
539  if (type_of_code == MFD)
540  get_MFD_gen_pol(inverse_rate, constraint_length, gen);
541  else if (type_of_code == ODS)
542  get_ODS_gen_pol(inverse_rate, constraint_length, gen);
543  else
544  it_assert(false, "This convolutional code doesn't exist in the tables");
545 
546  set_generator_polynomials(gen, constraint_length);
547 }
548 
549 /*
550  Set generator polynomials. Given in Proakis integer form
551 */
553  int constraint_length)
554 {
555  it_error_if(constraint_length <= 0, "Convolutional_Code::set_generator_polynomials(): Constraint length out of range");
556  gen_pol = gen;
557  n = gen.size();
558  it_error_if(n <= 0, "Convolutional_Code::set_generator_polynomials(): Invalid code rate");
559 
560  // Generate table look-up of weight of integers of size K bits
561  if (constraint_length != K) {
562  K = constraint_length;
563  xor_int_table.set_size(pow2i(K), false);
564  for (int i = 0; i < pow2i(K); i++) {
565  xor_int_table(i) = (weight_int(K, i) & 1);
566  }
567  }
568 
569  trunc_length = 5 * K;
570  m = K - 1;
571  no_states = pow2i(m);
572  gen_pol_rev.set_size(n, false);
573  rate = 1.0 / n;
574 
575  for (int i = 0; i < n; i++) {
576  gen_pol_rev(i) = reverse_int(K, gen_pol(i));
577  }
578 
579  int zero_output, one_output;
580  output_reverse_int.set_size(no_states, 2, false);
581 
582  for (int i = 0; i < no_states; i++) {
583  output_reverse(i, zero_output, one_output);
584  output_reverse_int(i, 0) = zero_output;
585  output_reverse_int(i, 1) = one_output;
586  }
587 
588  // initialise memory structures
589  visited_state.set_size(no_states);
590  visited_state = false;
591  visited_state(start_state) = true; // mark start state
592 
593  sum_metric.set_size(no_states);
594  sum_metric.clear();
595 
596  trunc_ptr = 0;
597  trunc_state = 0;
598 
599 }
600 
601 // Reset encoder and decoder states
603 {
604  init_encoder();
605 
606  visited_state = false;
607  visited_state(start_state) = true; // mark start state
608 
609  sum_metric.clear();
610 
611  trunc_ptr = 0;
612  trunc_state = 0;
613 }
614 
615 
616 /*
617  Encode a binary vector of inputs using specified method
618 */
619 void Convolutional_Code::encode(const bvec &input, bvec &output)
620 {
621  switch (cc_method) {
622  case Trunc:
623  encode_trunc(input, output);
624  break;
625  case Tailbite:
626  encode_tailbite(input, output);
627  break;
628  case Tail:
629  default:
630  encode_tail(input, output);
631  break;
632  };
633 }
634 
635 /*
636  Encode a binary vector of inputs starting from state set by the
637  set_state function
638 */
639 void Convolutional_Code::encode_trunc(const bvec &input, bvec &output)
640 {
641  int temp;
642  output.set_size(input.size() * n, false);
643 
644  for (int i = 0; i < input.size(); i++) {
645  encoder_state |= static_cast<int>(input(i)) << m;
646  for (int j = 0; j < n; j++) {
647  temp = encoder_state & gen_pol(j);
648  output(i * n + j) = xor_int_table(temp);
649  }
650  encoder_state >>= 1;
651  }
652 }
653 
654 /*
655  Encode a binary vector of inputs starting from zero state and also adds
656  a tail of K-1 zeros to force the encoder into the zero state. Well
657  suited for packet transmission.
658 */
659 void Convolutional_Code::encode_tail(const bvec &input, bvec &output)
660 {
661  int temp;
662  output.set_size((input.size() + m) * n, false);
663 
664  // always start from state 0
665  encoder_state = 0;
666 
667  for (int i = 0; i < input.size(); i++) {
668  encoder_state |= static_cast<int>(input(i)) << m;
669  for (int j = 0; j < n; j++) {
670  temp = encoder_state & gen_pol(j);
671  output(i * n + j) = xor_int_table(temp);
672  }
673  encoder_state >>= 1;
674  }
675 
676  // add tail of m = K-1 zeros
677  for (int i = input.size(); i < input.size() + m; i++) {
678  for (int j = 0; j < n; j++) {
679  temp = encoder_state & gen_pol(j);
680  output(i * n + j) = xor_int_table(temp);
681  }
682  encoder_state >>= 1;
683  }
684 }
685 
686 /*
687  Encode a binary vector of inputs starting using tail biting
688 */
689 void Convolutional_Code::encode_tailbite(const bvec &input, bvec &output)
690 {
691  int temp;
692  output.set_size(input.size() * n, false);
693 
694  // Set the start state equal to the end state:
695  encoder_state = 0;
696  bvec last_bits = input.right(m);
697  for (int i = 0; i < m; i++) {
698  encoder_state |= static_cast<int>(last_bits(i)) << m;
699  encoder_state >>= 1;
700  }
701 
702  for (int i = 0; i < input.size(); i++) {
703  encoder_state |= static_cast<int>(input(i)) << m;
704  for (int j = 0; j < n; j++) {
705  temp = encoder_state & gen_pol(j);
706  output(i * n + j) = xor_int_table(temp);
707  }
708  encoder_state >>= 1;
709  }
710 }
711 
712 /*
713  Encode a binary bit starting from the internal encoder state.
714  To initialize the encoder state use set_start_state() and init_encoder()
715 */
716 void Convolutional_Code::encode_bit(const bin &input, bvec &output)
717 {
718  int temp;
719  output.set_size(n, false);
720 
721  encoder_state |= static_cast<int>(input) << m;
722  for (int j = 0; j < n; j++) {
723  temp = encoder_state & gen_pol(j);
724  output(j) = xor_int_table(temp);
725  }
726  encoder_state >>= 1;
727 }
728 
729 
730 // --------------- Hard-decision decoding is not implemented -----------------
731 
732 void Convolutional_Code::decode(const bvec &, bvec &)
733 {
734  it_error("Convolutional_Code::decode(): Hard-decision decoding not implemented");
735 }
736 
737 bvec Convolutional_Code::decode(const bvec &)
738 {
739  it_error("Convolutional_Code::decode(): Hard-decision decoding not implemented");
740  return bvec();
741 }
742 
743 
744 /*
745  Decode a block of encoded data using specified method
746 */
747 void Convolutional_Code::decode(const vec &received_signal, bvec &output)
748 {
749  switch (cc_method) {
750  case Trunc:
751  decode_trunc(received_signal, output);
752  break;
753  case Tailbite:
754  decode_tailbite(received_signal, output);
755  break;
756  case Tail:
757  default:
758  decode_tail(received_signal, output);
759  break;
760  };
761 }
762 
763 /*
764  Decode a block of encoded data where encode_tail has been used.
765  Thus is assumes a decoder start state of zero and that a tail of
766  K-1 zeros has been added. No memory truncation.
767 */
768 void Convolutional_Code::decode_tail(const vec &received_signal, bvec &output)
769 {
770  int block_length = received_signal.size() / n; // no input symbols
771  it_error_if(block_length - m <= 0,
772  "Convolutional_Code::decode_tail(): Input sequence to short");
773  int S0, S1;
774  vec temp_sum_metric(no_states), temp_rec(n), delta_metrics;
775  Array<bool> temp_visited_state(no_states);
776  double temp_metric_zero, temp_metric_one;
777 
778  path_memory.set_size(no_states, block_length, false);
779  output.set_size(block_length - m, false); // no tail in the output
780 
781  // clear visited states
782  visited_state = false;
783  temp_visited_state = visited_state;
784  visited_state(0) = true; // starts in the zero state
785 
786  // clear accumulated metrics
787  sum_metric.clear();
788 
789  // evolve up to m where not all states visited
790  for (int l = 0; l < m; l++) { // all transitions including the tail
791  temp_rec = received_signal.mid(l * n, n);
792 
793  // calculate all metrics for all codewords at the same time
794  calc_metric(temp_rec, delta_metrics);
795 
796  for (int s = 0; s < no_states; s++) { // all states
797  // S0 and S1 are the states that expanded end at state s
798  previous_state(s, S0, S1);
799  if (visited_state(S0)) { // expand trellis if state S0 is visited
800  temp_metric_zero = sum_metric(S0)
801  + delta_metrics(output_reverse_int(s, 0));
802  temp_visited_state(s) = true;
803  }
804  else {
805  temp_metric_zero = std::numeric_limits<double>::max();
806  }
807  if (visited_state(S1)) { // expand trellis if state S0 is visited
808  temp_metric_one = sum_metric(S1)
809  + delta_metrics(output_reverse_int(s, 1));
810  temp_visited_state(s) = true;
811  }
812  else {
813  temp_metric_one = std::numeric_limits<double>::max();
814  }
815  if (temp_metric_zero < temp_metric_one) { // path zero survives
816  temp_sum_metric(s) = temp_metric_zero;
817  path_memory(s, l) = 0;
818  }
819  else { // path one survives
820  temp_sum_metric(s) = temp_metric_one;
821  path_memory(s, l) = 1;
822  }
823  } // all states, s
824  sum_metric = temp_sum_metric;
825  visited_state = temp_visited_state;
826  } // all transitions, l
827 
828  // evolve from m and to the end of the block
829  for (int l = m; l < block_length; l++) { // all transitions except the tail
830  temp_rec = received_signal.mid(l * n, n);
831 
832  // calculate all metrics for all codewords at the same time
833  calc_metric(temp_rec, delta_metrics);
834 
835  for (int s = 0; s < no_states; s++) { // all states
836  // S0 and S1 are the states that expanded end at state s
837  previous_state(s, S0, S1);
838  temp_metric_zero = sum_metric(S0)
839  + delta_metrics(output_reverse_int(s, 0));
840  temp_metric_one = sum_metric(S1)
841  + delta_metrics(output_reverse_int(s, 1));
842  if (temp_metric_zero < temp_metric_one) { // path zero survives
843  temp_sum_metric(s) = temp_metric_zero;
844  path_memory(s, l) = 0;
845  }
846  else { // path one survives
847  temp_sum_metric(s) = temp_metric_one;
848  path_memory(s, l) = 1;
849  }
850  } // all states, s
851  sum_metric = temp_sum_metric;
852  } // all transitions, l
853 
854  // minimum metric is the zeroth state due to the tail
855  int min_metric_state = 0;
856  // trace back to remove tail of zeros
857  for (int l = block_length - 1; l > block_length - 1 - m; l--) {
858  // previous state calculation
859  min_metric_state = previous_state(min_metric_state,
860  path_memory(min_metric_state, l));
861  }
862 
863  // trace back to calculate output sequence
864  for (int l = block_length - 1 - m; l >= 0; l--) {
865  output(l) = get_input(min_metric_state);
866  // previous state calculation
867  min_metric_state = previous_state(min_metric_state,
868  path_memory(min_metric_state, l));
869  }
870 }
871 
872 
873 /*
874  Decode a block of encoded data where encode_tailbite has been used.
875 */
876 void Convolutional_Code::decode_tailbite(const vec &received_signal,
877  bvec &output)
878 {
879  int block_length = received_signal.size() / n; // no input symbols
880  it_error_if(block_length <= 0,
881  "Convolutional_Code::decode_tailbite(): Input sequence to short");
882  int S0, S1;
883  vec temp_sum_metric(no_states), temp_rec(n), delta_metrics;
884  Array<bool> temp_visited_state(no_states);
885  double temp_metric_zero, temp_metric_one;
886  double best_metric = std::numeric_limits<double>::max();
887  bvec best_output(block_length), temp_output(block_length);
888 
889  path_memory.set_size(no_states, block_length, false);
890  output.set_size(block_length, false);
891 
892 
893  // Try all start states ss
894  for (int ss = 0; ss < no_states; ss++) {
895  //Clear the visited_state vector:
896  visited_state = false;
897  temp_visited_state = visited_state;
898  visited_state(ss) = true; // starts in the state s
899 
900  // clear accumulated metrics
901  sum_metric.zeros();
902 
903  for (int l = 0; l < block_length; l++) { // all transitions
904  temp_rec = received_signal.mid(l * n, n);
905  // calculate all metrics for all codewords at the same time
906  calc_metric(temp_rec, delta_metrics);
907 
908  for (int s = 0; s < no_states; s++) { // all states
909  // S0 and S1 are the states that expanded end at state s
910  previous_state(s, S0, S1);
911  if (visited_state(S0)) { // expand trellis if state S0 is visited
912  temp_metric_zero = sum_metric(S0)
913  + delta_metrics(output_reverse_int(s, 0));
914  temp_visited_state(s) = true;
915  }
916  else {
917  temp_metric_zero = std::numeric_limits<double>::max();
918  }
919  if (visited_state(S1)) { // expand trellis if state S0 is visited
920  temp_metric_one = sum_metric(S1)
921  + delta_metrics(output_reverse_int(s, 1));
922  temp_visited_state(s) = true;
923  }
924  else {
925  temp_metric_one = std::numeric_limits<double>::max();
926  }
927  if (temp_metric_zero < temp_metric_one) { // path zero survives
928  temp_sum_metric(s) = temp_metric_zero;
929  path_memory(s, l) = 0;
930  }
931  else { // path one survives
932  temp_sum_metric(s) = temp_metric_one;
933  path_memory(s, l) = 1;
934  }
935  } // all states, s
936  sum_metric = temp_sum_metric;
937  visited_state = temp_visited_state;
938  } // all transitions, l
939 
940  // minimum metric is the ss state due to the tailbite
941  int min_metric_state = ss;
942 
943  // trace back to calculate output sequence
944  for (int l = block_length - 1; l >= 0; l--) {
945  temp_output(l) = get_input(min_metric_state);
946  // previous state calculation
947  min_metric_state = previous_state(min_metric_state,
948  path_memory(min_metric_state, l));
949  }
950  if (sum_metric(ss) < best_metric) {
951  best_metric = sum_metric(ss);
952  best_output = temp_output;
953  }
954  } // all start states ss
955  output = best_output;
956 }
957 
958 
959 /*
960  Viterbi decoding using truncation of memory (default = 5*K)
961 */
962 void Convolutional_Code::decode_trunc(const vec &received_signal,
963  bvec &output)
964 {
965  int block_length = received_signal.size() / n; // no input symbols
966  it_error_if(block_length <= 0,
967  "Convolutional_Code::decode_trunc(): Input sequence to short");
968  int S0, S1;
969  vec temp_sum_metric(no_states), temp_rec(n), delta_metrics;
970  Array<bool> temp_visited_state(no_states);
971  double temp_metric_zero, temp_metric_one;
972 
973  path_memory.set_size(no_states, trunc_length, false);
974  output.set_size(0);
975 
976  // copy visited states
977  temp_visited_state = visited_state;
978 
979  for (int i = 0; i < block_length; i++) {
980  // update path memory pointer
981  trunc_ptr = (trunc_ptr + 1) % trunc_length;
982 
983  temp_rec = received_signal.mid(i * n, n);
984  // calculate all metrics for all codewords at the same time
985  calc_metric(temp_rec, delta_metrics);
986 
987  for (int s = 0; s < no_states; s++) { // all states
988  // the states that expanded end at state s
989  previous_state(s, S0, S1);
990  if (visited_state(S0)) { // expand trellis if state S0 is visited
991  temp_metric_zero = sum_metric(S0)
992  + delta_metrics(output_reverse_int(s, 0));
993  temp_visited_state(s) = true;
994  }
995  else {
996  temp_metric_zero = std::numeric_limits<double>::max();
997  }
998  if (visited_state(S1)) { // expand trellis if state S0 is visited
999  temp_metric_one = sum_metric(S1)
1000  + delta_metrics(output_reverse_int(s, 1));
1001  temp_visited_state(s) = true;
1002  }
1003  else {
1004  temp_metric_one = std::numeric_limits<double>::max();
1005  }
1006  if (temp_metric_zero < temp_metric_one) { // path zero survives
1007  temp_sum_metric(s) = temp_metric_zero;
1008  path_memory(s, trunc_ptr) = 0;
1009  }
1010  else { // path one survives
1011  temp_sum_metric(s) = temp_metric_one;
1012  path_memory(s, trunc_ptr) = 1;
1013  }
1014  } // all states, s
1015  sum_metric = temp_sum_metric;
1016  visited_state = temp_visited_state;
1017 
1018  // find minimum metric
1019  int min_metric_state = min_index(sum_metric);
1020 
1021  // normalise accumulated metrics
1022  sum_metric -= sum_metric(min_metric_state);
1023 
1024  // check if we had enough metrics to generate output
1025  if (trunc_state >= trunc_length) {
1026  // trace back to calculate the output symbol
1027  for (int j = trunc_length; j > 0; j--) {
1028  // previous state calculation
1029  min_metric_state =
1030  previous_state(min_metric_state,
1031  path_memory(min_metric_state,
1032  (trunc_ptr + j) % trunc_length));
1033  }
1034  output.ins(output.size(), get_input(min_metric_state));
1035  }
1036  else { // if not increment trunc_state counter
1037  trunc_state++;
1038  }
1039  } // end for (int i = 0; i < block_length; l++)
1040 }
1041 
1042 
1043 /*
1044  Calculate the inverse sequence
1045 
1046  Assumes that encode_tail is used in the encoding process. Returns false
1047  if there is an error in the coded sequence (not a valid codeword). Do
1048  not check that the tail forces the encoder into the zeroth state.
1049 */
1050 bool Convolutional_Code::inverse_tail(const bvec coded_sequence, bvec &input)
1051 {
1052  int state = 0, zero_state, one_state, zero_temp, one_temp, i, j;
1053  bvec zero_output(n), one_output(n);
1054 
1055  int block_length = coded_sequence.size() / n - m; // no input symbols
1056  it_error_if(block_length <= 0, "The input sequence is to short");
1057  input.set_length(block_length, false); // not include the tail in the output
1058 
1059 
1060  for (i = 0; i < block_length; i++) {
1061  zero_state = state;
1062  one_state = state | (1 << m);
1063  for (j = 0; j < n; j++) {
1064  zero_temp = zero_state & gen_pol(j);
1065  one_temp = one_state & gen_pol(j);
1066  zero_output(j) = xor_int_table(zero_temp);
1067  one_output(j) = xor_int_table(one_temp);
1068  }
1069  if (coded_sequence.mid(i*n, n) == zero_output) {
1070  input(i) = bin(0);
1071  state = zero_state >> 1;
1072  }
1073  else if (coded_sequence.mid(i*n, n) == one_output) {
1074  input(i) = bin(1);
1075  state = one_state >> 1;
1076  }
1077  else
1078  return false;
1079  }
1080 
1081  return true;
1082 }
1083 
1084 /*
1085  Check if catastrophic. Returns true if catastrophic
1086 */
1088 {
1089  int start, S, W0, W1, S0, S1;
1090  bvec visited(1 << m);
1091 
1092  for (start = 1; start < no_states; start++) {
1093  visited.zeros();
1094  S = start;
1095  visited(S) = 1;
1096 
1097  node1:
1098  S0 = next_state(S, 0);
1099  S1 = next_state(S, 1);
1100  weight(S, W0, W1);
1101  if (S1 < start) goto node0;
1102  if (W1 > 0) goto node0;
1103 
1104  if (visited(S1) == bin(1))
1105  return true;
1106  S = S1; // goto node1
1107  visited(S) = 1;
1108  goto node1;
1109 
1110  node0:
1111  if (S0 < start) continue; // goto end;
1112  if (W0 > 0) continue; // goto end;
1113 
1114  if (visited(S0) == bin(1))
1115  return true;
1116  S = S0;
1117  visited(S) = 1;
1118  goto node1;
1119 
1120  //end:
1121  //void;
1122  }
1123 
1124  return false;
1125 }
1126 
1127 /*
1128  Calculate distance profile. If reverse = true calculate for the reverse code instead.
1129 */
1130 void Convolutional_Code::distance_profile(ivec &dist_prof, int dmax, bool reverse)
1131 {
1132  int max_stack_size = 50000;
1133  ivec S_stack(max_stack_size), W_stack(max_stack_size), t_stack(max_stack_size);
1134 
1135  int stack_pos = -1, t, S, W, W0, w0, w1;
1136 
1137  dist_prof.set_size(K, false);
1138  dist_prof.zeros();
1139  dist_prof += dmax; // just select a big number!
1140  if (reverse)
1141  W = weight_reverse(0, 1);
1142  else
1143  W = weight(0, 1);
1144 
1145  S = next_state(0, 1); // first state 0 and one as input, S = 1<<(m-1);
1146  t = 0;
1147  dist_prof(0) = W;
1148 
1149 node1:
1150  if (reverse)
1151  weight_reverse(S, w0, w1);
1152  else
1153  weight(S, w0, w1);
1154 
1155  if (t < m) {
1156  W0 = W + w0;
1157  if (W0 < dist_prof(m)) { // store node0 (S, W0, and t+1)
1158  stack_pos++;
1159  if (stack_pos >= max_stack_size) {
1160  max_stack_size = 2 * max_stack_size;
1161  S_stack.set_size(max_stack_size, true);
1162  W_stack.set_size(max_stack_size, true);
1163  t_stack.set_size(max_stack_size, true);
1164  }
1165 
1166  S_stack(stack_pos) = next_state(S, 0); //S>>1;
1167  W_stack(stack_pos) = W0;
1168  t_stack(stack_pos) = t + 1;
1169  }
1170  }
1171  else goto stack;
1172 
1173  W += w1;
1174  if (W > dist_prof(m))
1175  goto stack;
1176 
1177  t++;
1178  S = next_state(S, 1); //S = (S>>1)|(1<<(m-1));
1179 
1180  if (W < dist_prof(t))
1181  dist_prof(t) = W;
1182 
1183  if (t == m) goto stack;
1184  goto node1;
1185 
1186 
1187 stack:
1188  if (stack_pos >= 0) {
1189  // get root node from stack
1190  S = S_stack(stack_pos);
1191  W = W_stack(stack_pos);
1192  t = t_stack(stack_pos);
1193  stack_pos--;
1194 
1195  if (W < dist_prof(t))
1196  dist_prof(t) = W;
1197 
1198  if (t == m) goto stack;
1199 
1200  goto node1;
1201  }
1202 }
1203 
1204 /*
1205  Calculate spectrum
1206 
1207  Calculates both the weight spectrum (Ad) and the information weight spectrum (Cd) and
1208  returns it as ivec:s in the 0:th and 1:st component of spectrum, respectively. Suitable
1209  for calculating many terms in the spectra (uses an breadth first algorithm). It is assumed
1210  that the code is non-catastrophic or else it is a possibility for an eternal loop.
1211  dmax = an upper bound on the free distance
1212  no_terms = no_terms including the dmax term that should be calculated
1213 */
1215 {
1216  imat Ad_states(no_states, dmax + no_terms), Cd_states(no_states, dmax + no_terms);
1217  imat Ad_temp(no_states, dmax + no_terms), Cd_temp(no_states, dmax + no_terms);
1218  ivec mindist(no_states), mindist_temp(1 << m);
1219 
1220  spectrum.set_size(2);
1221  spectrum(0).set_size(dmax + no_terms, false);
1222  spectrum(1).set_size(dmax + no_terms, false);
1223  spectrum(0).zeros();
1224  spectrum(1).zeros();
1225  Ad_states.zeros();
1226  Cd_states.zeros();
1227  mindist.zeros();
1228  int wmax = dmax + no_terms;
1229  ivec visited_states(no_states), visited_states_temp(no_states);
1230  bool proceede;
1231  int d, w0, w1, s, s0, s1;
1232 
1233  visited_states.zeros(); // 0 = false
1234  s = next_state(0, 1); // Start in state from 0 with an one input.
1235  visited_states(s) = 1; // 1 = true. Start in state 2^(m-1).
1236  w1 = weight(0, 1);
1237  Ad_states(s, w1) = 1;
1238  Cd_states(s, w1) = 1;
1239  mindist(s) = w1;
1240  proceede = true;
1241 
1242  while (proceede) {
1243  proceede = false;
1244  Ad_temp.zeros();
1245  Cd_temp.zeros();
1246  mindist_temp.zeros();
1247  visited_states_temp.zeros(); //false
1248  for (s = 1; s < no_states; s++) {
1249  if ((mindist(s) > 0) && (mindist(s) < wmax)) {
1250  proceede = true;
1251  weight(s, w0, w1);
1252  s0 = next_state(s, 0);
1253  for (d = mindist(s); d < (wmax - w0); d++) {
1254  Ad_temp(s0, d + w0) += Ad_states(s, d);
1255  Cd_temp(s0, d + w0) += Cd_states(s, d);
1256  visited_states_temp(s0) = 1; //true
1257  }
1258 
1259  s1 = next_state(s, 1);
1260  for (d = mindist(s); d < (wmax - w1); d++) {
1261  Ad_temp(s1, d + w1) += Ad_states(s, d);
1262  Cd_temp(s1, d + w1) += Cd_states(s, d) + Ad_states(s, d);
1263  visited_states_temp(s1) = 1; //true
1264  }
1265  if (mindist_temp(s0) > 0) { mindist_temp(s0) = (mindist(s) + w0) < mindist_temp(s0) ? mindist(s) + w0 : mindist_temp(s0); }
1266  else { mindist_temp(s0) = mindist(s) + w0; }
1267  if (mindist_temp(s1) > 0) { mindist_temp(s1) = (mindist(s) + w1) < mindist_temp(s1) ? mindist(s) + w1 : mindist_temp(s1); }
1268  else { mindist_temp(s1) = mindist(s) + w1; }
1269 
1270  }
1271  }
1272  Ad_states = Ad_temp;
1273  Cd_states = Cd_temp;
1274  spectrum(0) += Ad_temp.get_row(0);
1275  spectrum(1) += Cd_temp.get_row(0);
1276  visited_states = visited_states_temp;
1277  mindist = elem_mult(mindist_temp, visited_states);
1278  }
1279 }
1280 
1281 /*
1282  Cederwall's fast algorithm
1283 
1284  See IT No. 6, pp. 1146-1159, Nov. 1989 for details.
1285  Calculates both the weight spectrum (Ad) and the information weight spectrum (Cd) and
1286  returns it as ivec:s in the 0:th and 1:st component of spectrum, respectively. The FAST algorithm
1287  is good for calculating only a few terms in the spectrum. If many terms are desired, use calc_spectrum instead.
1288  The algorithm returns -1 if the code tested is worse that the input dfree and Cdfree.
1289  It returns 0 if the code MAY be catastrophic (assuming that test_catastrophic is true),
1290  and returns 1 if everything went right.
1291 
1292  \arg \c dfree the free distance of the code (or an upper bound)
1293  \arg \c no_terms including the dfree term that should be calculated
1294  \ar \c Cdfree is the best value of information weight spectrum found so far
1295 */
1296 int Convolutional_Code::fast(Array<ivec> &spectrum, const int dfree, const int no_terms, const int Cdfree, const bool test_catastrophic)
1297 {
1298  int cat_treshold = 7 * K; // just a big number, but not to big!
1299  int i;
1300  ivec dist_prof(K), dist_prof_rev(K);
1301  distance_profile(dist_prof, dfree);
1302  distance_profile(dist_prof_rev, dfree, true); // for the reverse code
1303 
1304  int dist_prof_rev0 = dist_prof_rev(0);
1305 
1306  bool reverse = false; // false = use dist_prof
1307 
1308  // is the reverse distance profile better?
1309  for (i = 0; i < K; i++) {
1310  if (dist_prof_rev(i) > dist_prof(i)) {
1311  reverse = true;
1312  dist_prof_rev0 = dist_prof(0);
1313  dist_prof = dist_prof_rev;
1314  break;
1315  }
1316  else if (dist_prof_rev(i) < dist_prof(i)) {
1317  break;
1318  }
1319  }
1320 
1321  int d = dfree + no_terms - 1;
1322  int max_stack_size = 50000;
1323  ivec S_stack(max_stack_size), W_stack(max_stack_size), b_stack(max_stack_size), c_stack(max_stack_size);
1324  int stack_pos = -1;
1325 
1326  // F1.
1327  int mf = 1, b = 1;
1328  spectrum.set_size(2);
1329  spectrum(0).set_size(dfree + no_terms, false); // ad
1330  spectrum(1).set_size(dfree + no_terms, false); // cd
1331  spectrum(0).zeros();
1332  spectrum(1).zeros();
1333  int S, S0, S1, W0, W1, w0, w1, c = 0;
1334  S = next_state(0, 1); //first state zero and one as input
1335  int W = d - dist_prof_rev0;
1336 
1337 
1338 F2:
1339  S0 = next_state(S, 0);
1340  S1 = next_state(S, 1);
1341 
1342  if (reverse) {
1343  weight(S, w0, w1);
1344  }
1345  else {
1346  weight_reverse(S, w0, w1);
1347  }
1348  W0 = W - w0;
1349  W1 = W - w1;
1350  if (mf < m) goto F6;
1351 
1352  //F3:
1353  if (W0 >= 0) {
1354  spectrum(0)(d - W0)++;
1355  spectrum(1)(d - W0) += b;
1356  // The code is worse than the best found.
1357  if (((d - W0) < dfree) || (((d - W0) == dfree) && (spectrum(1)(d - W0) > Cdfree)))
1358  return -1;
1359  }
1360 
1361 
1362 F4:
1363  if ((W1 < dist_prof(m - 1)) || (W < dist_prof(m))) goto F5;
1364  // select node 1
1365  if (test_catastrophic && W == W1) {
1366  c++;
1367  if (c > cat_treshold)
1368  return 0;
1369  }
1370  else {
1371  c = 0;
1372  W = W1;
1373  }
1374 
1375  S = S1;
1376  mf = 1;
1377  b++;
1378  goto F2;
1379 
1380 F5:
1381  //if (stack_pos == -1) return n_values;
1382  if (stack_pos == -1) goto end;
1383  // get node from stack
1384  S = S_stack(stack_pos);
1385  W = W_stack(stack_pos);
1386  b = b_stack(stack_pos);
1387  c = c_stack(stack_pos);
1388  stack_pos--;
1389  mf = 1;
1390  goto F2;
1391 
1392 F6:
1393  if (W0 < dist_prof(m - mf - 1)) goto F4;
1394 
1395  //F7:
1396  if ((W1 >= dist_prof(m - 1)) && (W >= dist_prof(m))) {
1397  // save node 1
1398  if (test_catastrophic && stack_pos > 10000)
1399  return 0;
1400 
1401  stack_pos++;
1402  if (stack_pos >= max_stack_size) {
1403  max_stack_size = 2 * max_stack_size;
1404  S_stack.set_size(max_stack_size, true);
1405  W_stack.set_size(max_stack_size, true);
1406  b_stack.set_size(max_stack_size, true);
1407  c_stack.set_size(max_stack_size, true);
1408  }
1409  S_stack(stack_pos) = S1;
1410  W_stack(stack_pos) = W1;
1411  b_stack(stack_pos) = b + 1; // because gone to node 1
1412  c_stack(stack_pos) = c; // because gone to node 1
1413  }
1414  // select node 0
1415  S = S0;
1416  if (test_catastrophic && W == W0) {
1417  c++;
1418  if (c > cat_treshold)
1419  return 0;
1420  }
1421  else {
1422  c = 0;
1423  W = W0;
1424  }
1425 
1426 
1427  mf++;
1428  goto F2;
1429 
1430 end:
1431  return 1;
1432 }
1433 
1434 //----------- These functions should be moved into some other place -------
1435 
1439 int reverse_int(int length, int in)
1440 {
1441  int i, j, out = 0;
1442 
1443  for (i = 0; i < (length >> 1); i++) {
1444  out = out | ((in & (1 << i)) << (length - 1 - (i << 1)));
1445  }
1446  for (j = 0; j < (length - i); j++) {
1447  out = out | ((in & (1 << (j + i))) >> ((j << 1) - (length & 1) + 1));
1448  //out = out | ( (in & (1<<j+i)) >> ((j<<1)-(length&1)+1) ); old version with preecedence problems in MSC
1449 
1450  }
1451  return out;
1452 }
1453 
1457 int weight_int(int length, int in)
1458 {
1459  int i, w = 0;
1460  for (i = 0; i < length; i++) {
1461  w += (in & (1 << i)) >> i;
1462  }
1463  return w;
1464 }
1465 
1469 int compare_spectra(ivec v1, ivec v2)
1470 {
1471  it_assert_debug(v1.size() == v2.size(), "compare_spectra: wrong sizes");
1472 
1473  for (int i = 0; i < v1.size(); i++) {
1474  if (v1(i) < v2(i)) {
1475  return 1;
1476  }
1477  else if (v1(i) > v2(i)) {
1478  return 0;
1479  }
1480  }
1481  return -1;
1482 }
1483 
1489 int compare_spectra(ivec v1, ivec v2, vec weight_profile)
1490 {
1491  double t1 = 0, t2 = 0;
1492  for (int i = 0; i < v1.size(); i++) {
1493  t1 += (double)v1(i) * weight_profile(i);
1494  t2 += (double)v2(i) * weight_profile(i);
1495  }
1496  if (t1 < t2) return 1;
1497  else if (t1 > t2) return 0;
1498  else return -1;
1499 }
1500 
1501 } // namespace itpp
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:13 for IT++ by Doxygen 1.8.1.2