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

Generated on Sat Jul 6 2013 10:54:22 for IT++ by Doxygen 1.8.2