IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
rec_syst_conv_code.cpp
Go to the documentation of this file.
1 
30 
31 
32 namespace itpp
33 {
34 
36 double(*com_log)(double, double) = NULL;
37 
39 // This wrapper is because "com_log = std::max;" below caused an error
40 inline double max(double x, double y) { return std::max(x, y); }
42 
43 // ----------------- Protected functions -----------------------------
44 
45 int Rec_Syst_Conv_Code::calc_state_transition(const int instate, const int input, ivec &parity)
46 {
47  int i, j, in = 0, temp = (gen_pol_rev(0) & (instate << 1)), parity_temp, parity_bit;
48 
49  for (i = 0; i < K; i++) {
50  in = (temp & 1) ^ in;
51  temp = temp >> 1;
52  }
53  in = in ^ input;
54 
55  parity.set_size(n - 1, false);
56  for (j = 0; j < (n - 1); j++) {
57  parity_temp = ((instate << 1) + in) & gen_pol_rev(j + 1);
58  parity_bit = 0;
59  for (i = 0; i < K; i++) {
60  parity_bit = (parity_temp & 1) ^ parity_bit;
61  parity_temp = parity_temp >> 1;
62  }
63  parity(j) = parity_bit;
64  }
65  return in + ((instate << 1) & ((1 << m) - 1));
66 }
67 
68 // --------------- Public functions -------------------------
69 void Rec_Syst_Conv_Code::set_generator_polynomials(const ivec &gen, int constraint_length)
70 {
71  int j;
72  gen_pol = gen;
73  n = gen.size();
74  K = constraint_length;
75  m = K - 1;
76  rate = 1.0 / n;
77 
78  gen_pol_rev.set_size(n, false);
79  for (int i = 0; i < n; i++) {
80  gen_pol_rev(i) = reverse_int(K, gen_pol(i));
81  }
82 
83  Nstates = (1 << m);
84  state_trans.set_size(Nstates, 2, false);
85  rev_state_trans.set_size(Nstates, 2, false);
86  output_parity.set_size(Nstates, 2*(n - 1), false);
87  rev_output_parity.set_size(Nstates, 2*(n - 1), false);
88  int s0, s1, s_prim;
89  ivec p0, p1;
90  for (s_prim = 0; s_prim < Nstates; s_prim++) {
91  s0 = calc_state_transition(s_prim, 0, p0);
92  state_trans(s_prim, 0) = s0;
93  rev_state_trans(s0, 0) = s_prim;
94  for (j = 0; j < (n - 1); j++) {
95  output_parity(s_prim, 2*j + 0) = p0(j);
96  rev_output_parity(s0, 2*j + 0) = p0(j);
97  }
98 
99  s1 = calc_state_transition(s_prim, 1, p1);
100  state_trans(s_prim, 1) = s1;
101  rev_state_trans(s1, 1) = s_prim;
102  for (j = 0; j < (n - 1); j++) {
103  output_parity(s_prim, 2*j + 1) = p1(j);
104  rev_output_parity(s1, 2*j + 1) = p1(j);
105  }
106  }
107 
108  ln2 = std::log(2.0);
109 
110  //The default value of Lc is 1:
111  Lc = 1.0;
112 }
113 
115 {
116  Lc = 4.0 * std::sqrt(Ec) / N0;
117 }
118 
120 {
121  Lc = in_Lc;
122 }
123 
124 void Rec_Syst_Conv_Code::encode_tail(const bvec &input, bvec &tail, bmat &parity_bits)
125 {
126  int i, j, length = input.size(), target_state;
127  parity_bits.set_size(length + m, n - 1, false);
128  tail.set_size(m, false);
129 
130  encoder_state = 0;
131  for (i = 0; i < length; i++) {
132  for (j = 0; j < (n - 1); j++) {
133  parity_bits(i, j) = output_parity(encoder_state, 2 * j + int(input(i)));
134  }
135  encoder_state = state_trans(encoder_state, int(input(i)));
136  }
137 
138  // add tail of m=K-1 zeros
139  for (i = 0; i < m; i++) {
140  target_state = (encoder_state << 1) & ((1 << m) - 1);
141  if (state_trans(encoder_state, 0) == target_state) { tail(i) = bin(0); }
142  else { tail(i) = bin(1); }
143  for (j = 0; j < (n - 1); j++) {
144  parity_bits(length + i, j) = output_parity(encoder_state, 2 * j + int(tail(i)));
145  }
146  encoder_state = target_state;
147  }
148  terminated = true;
149 }
150 
151 void Rec_Syst_Conv_Code::encode(const bvec &input, bmat &parity_bits)
152 {
153  int i, j, length = input.size();
154  parity_bits.set_size(length, n - 1, false);
155 
156  encoder_state = 0;
157  for (i = 0; i < length; i++) {
158  for (j = 0; j < (n - 1); j++) {
159  parity_bits(i, j) = output_parity(encoder_state, 2 * j + int(input(i)));
160  }
161  encoder_state = state_trans(encoder_state, int(input(i)));
162  }
163  terminated = false;
164 }
165 
166 void Rec_Syst_Conv_Code::map_decode(const vec &rec_systematic, const mat &rec_parity, const vec &extrinsic_input,
167  vec &extrinsic_output, bool in_terminated)
168 {
169  double gamma_k_e, nom, den, temp0, temp1, exp_temp0, exp_temp1;
170  int j, s0, s1, k, kk, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
171  ivec p0, p1;
172 
173  alpha.set_size(Nstates, block_length + 1, false);
174  beta.set_size(Nstates, block_length + 1, false);
175  gamma.set_size(2*Nstates, block_length + 1, false);
176  denom.set_size(block_length + 1, false);
177  denom.clear();
178  extrinsic_output.set_size(block_length, false);
179 
180  if (in_terminated) { terminated = true; }
181 
182  //Calculate gamma
183  for (k = 1; k <= block_length; k++) {
184  kk = k - 1;
185  for (s_prim = 0; s_prim < Nstates; s_prim++) {
186  exp_temp0 = 0.0;
187  exp_temp1 = 0.0;
188  for (j = 0; j < (n - 1); j++) {
189  exp_temp0 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 0)); /* Is this OK? */
190  exp_temp1 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 1));
191  }
192  // gamma(2*s_prim+0,k) = std::exp( 0.5*(extrinsic_input(kk) + Lc*rec_systematic(kk))) * std::exp( exp_temp0 );
193  // gamma(2*s_prim+1,k) = std::exp(-0.5*(extrinsic_input(kk) + Lc*rec_systematic(kk))) * std::exp( exp_temp1 );
194  /* == Changed to trunc_exp() to address bug 1088420 which
195  described a numerical instability problem in map_decode()
196  at high SNR. This should be regarded as a temporary fix and
197  it is not necessarily a waterproof one: multiplication of
198  probabilities still can result in values out of
199  range. (Range checking for the multiplication operator was
200  not implemented as it was felt that this would sacrifice
201  too much runtime efficiency. Some margin was added to the
202  numerical hardlimits below to reflect this. The hardlimit
203  values below were taken as the minimum range that a
204  "double" should support reduced by a few orders of
205  magnitude to make sure multiplication of several values
206  does not exceed the limits.)
207 
208  It is suggested to use the QLLR based log-domain() decoders
209  instead of map_decode() as they are much faster and more
210  numerically stable.
211 
212  EGL 8/06. == */
213  gamma(2*s_prim + 0, k) = trunc_exp(0.5 * (extrinsic_input(kk) + Lc * rec_systematic(kk)) + exp_temp0);
214  gamma(2*s_prim + 1, k) = trunc_exp(-0.5 * (extrinsic_input(kk) + Lc * rec_systematic(kk)) + exp_temp1);
215  }
216  }
217 
218  //Initiate alpha
219  alpha.set_col(0, zeros(Nstates));
220  alpha(0, 0) = 1.0;
221 
222  //Calculate alpha and denom going forward through the trellis
223  for (k = 1; k <= block_length; k++) {
224  for (s = 0; s < Nstates; s++) {
225  s_prim0 = rev_state_trans(s, 0);
226  s_prim1 = rev_state_trans(s, 1);
227  temp0 = alpha(s_prim0, k - 1) * gamma(2 * s_prim0 + 0, k);
228  temp1 = alpha(s_prim1, k - 1) * gamma(2 * s_prim1 + 1, k);
229  alpha(s, k) = temp0 + temp1;
230  denom(k) += temp0 + temp1;
231  }
232  alpha.set_col(k, alpha.get_col(k) / denom(k));
233  }
234 
235  //Initiate beta
236  if (terminated) {
237  beta.set_col(block_length, zeros(Nstates));
238  beta(0, block_length) = 1.0;
239  }
240  else {
241  beta.set_col(block_length, alpha.get_col(block_length));
242  }
243 
244  //Calculate beta going backward in the trellis
245  for (k = block_length; k >= 2; k--) {
246  for (s_prim = 0; s_prim < Nstates; s_prim++) {
247  s0 = state_trans(s_prim, 0);
248  s1 = state_trans(s_prim, 1);
249  beta(s_prim, k - 1) = (beta(s0, k) * gamma(2 * s_prim + 0, k)) + (beta(s1, k) * gamma(2 * s_prim + 1, k));
250  }
251  beta.set_col(k - 1, beta.get_col(k - 1) / denom(k));
252  }
253 
254  //Calculate extrinsic output for each bit
255  for (k = 1; k <= block_length; k++) {
256  kk = k - 1;
257  nom = 0;
258  den = 0;
259  for (s_prim = 0; s_prim < Nstates; s_prim++) {
260  s0 = state_trans(s_prim, 0);
261  s1 = state_trans(s_prim, 1);
262  exp_temp0 = 0.0;
263  exp_temp1 = 0.0;
264  for (j = 0; j < (n - 1); j++) {
265  exp_temp0 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 0));
266  exp_temp1 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 1));
267  }
268  // gamma_k_e = std::exp( exp_temp0 );
269  gamma_k_e = trunc_exp(exp_temp0);
270  nom += alpha(s_prim, k - 1) * gamma_k_e * beta(s0, k);
271 
272  // gamma_k_e = std::exp( exp_temp1 );
273  gamma_k_e = trunc_exp(exp_temp1);
274  den += alpha(s_prim, k - 1) * gamma_k_e * beta(s1, k);
275  }
276  // extrinsic_output(kk) = std::log(nom/den);
277  extrinsic_output(kk) = trunc_log(nom / den);
278  }
279 
280 }
281 
282 void Rec_Syst_Conv_Code::log_decode(const vec &rec_systematic, const mat &rec_parity,
283  const vec &extrinsic_input, vec &extrinsic_output, bool in_terminated, std::string metric)
284 {
285  if (metric == "TABLE") {
286  /* Use the QLLR decoder. This can probably be done more
287  efficiently since it converts floating point vectors to QLLR.
288  However we have to live with this for the time being. */
289  QLLRvec rec_systematic_q = llrcalc.to_qllr(rec_systematic);
290  QLLRmat rec_parity_q = llrcalc.to_qllr(rec_parity);
291  QLLRvec extrinsic_input_q = llrcalc.to_qllr(extrinsic_input);
292  QLLRvec extrinsic_output_q(length(extrinsic_output));
293  log_decode(rec_systematic_q, rec_parity_q, extrinsic_input_q,
294  extrinsic_output_q, in_terminated);
295  extrinsic_output = llrcalc.to_double(extrinsic_output_q);
296  return;
297  }
298 
299  double nom, den, exp_temp0, exp_temp1, rp, temp0, temp1;
300  int i, j, s0, s1, k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
301  ivec p0, p1;
302 
303  //Set the internal metric:
304  if (metric == "LOGMAX") { com_log = max; }
305  else if (metric == "LOGMAP") { com_log = log_add; }
306  else {
307  it_error("Rec_Syst_Conv_Code::log_decode: Illegal metric parameter");
308  }
309 
310  alpha.set_size(Nstates, block_length + 1, false);
311  beta.set_size(Nstates, block_length + 1, false);
312  gamma.set_size(2*Nstates, block_length + 1, false);
313  extrinsic_output.set_size(block_length, false);
314  denom.set_size(block_length + 1, false);
315  for (k = 0; k <= block_length; k++) { denom(k) = -infinity; }
316 
317  if (in_terminated) { terminated = true; }
318 
319  //Check that Lc = 1.0
320  it_assert(Lc == 1.0,
321  "Rec_Syst_Conv_Code::log_decode: This function assumes that Lc = 1.0. Please use proper scaling of the input data");
322 
323  //Calculate gamma
324  for (k = 1; k <= block_length; k++) {
325  kk = k - 1;
326  for (s_prim = 0; s_prim < Nstates; s_prim++) {
327  exp_temp0 = 0.0;
328  exp_temp1 = 0.0;
329  for (j = 0; j < (n - 1); j++) {
330  rp = rec_parity(kk, j);
331  if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; }
332  else { exp_temp0 -= rp; }
333  if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; }
334  else { exp_temp1 -= rp; }
335  }
336  gamma(2*s_prim + 0, k) = 0.5 * ((extrinsic_input(kk) + rec_systematic(kk)) + exp_temp0);
337  gamma(2*s_prim + 1, k) = -0.5 * ((extrinsic_input(kk) + rec_systematic(kk)) - exp_temp1);
338  }
339  }
340 
341  //Initiate alpha
342  for (j = 1; j < Nstates; j++) { alpha(j, 0) = -infinity; }
343  alpha(0, 0) = 0.0;
344 
345  //Calculate alpha, going forward through the trellis
346  for (k = 1; k <= block_length; k++) {
347  for (s = 0; s < Nstates; s++) {
348  s_prim0 = rev_state_trans(s, 0);
349  s_prim1 = rev_state_trans(s, 1);
350  temp0 = alpha(s_prim0, k - 1) + gamma(2 * s_prim0 + 0, k);
351  temp1 = alpha(s_prim1, k - 1) + gamma(2 * s_prim1 + 1, k);
352  alpha(s, k) = com_log(temp0, temp1);
353  denom(k) = com_log(alpha(s, k), denom(k));
354  }
355  //Normalization of alpha
356  for (l = 0; l < Nstates; l++) { alpha(l, k) -= denom(k); }
357  }
358 
359  //Initiate beta
360  if (terminated) {
361  for (i = 1; i < Nstates; i++) { beta(i, block_length) = -infinity; }
362  beta(0, block_length) = 0.0;
363  }
364  else {
365  beta.set_col(block_length, alpha.get_col(block_length));
366  }
367 
368  //Calculate beta going backward in the trellis
369  for (k = block_length; k >= 1; k--) {
370  for (s_prim = 0; s_prim < Nstates; s_prim++) {
371  s0 = state_trans(s_prim, 0);
372  s1 = state_trans(s_prim, 1);
373  beta(s_prim, k - 1) = com_log(beta(s0, k) + gamma(2 * s_prim + 0, k) , beta(s1, k) + gamma(2 * s_prim + 1, k));
374  }
375  //Normalization of beta
376  for (l = 0; l < Nstates; l++) { beta(l, k - 1) -= denom(k); }
377  }
378 
379  //Calculate extrinsic output for each bit
380  for (k = 1; k <= block_length; k++) {
381  kk = k - 1;
382  nom = -infinity;
383  den = -infinity;
384  for (s_prim = 0; s_prim < Nstates; s_prim++) {
385  s0 = state_trans(s_prim, 0);
386  s1 = state_trans(s_prim, 1);
387  exp_temp0 = 0.0;
388  exp_temp1 = 0.0;
389  for (j = 0; j < (n - 1); j++) {
390  rp = rec_parity(kk, j);
391  if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; }
392  else { exp_temp0 -= rp; }
393  if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; }
394  else { exp_temp1 -= rp; }
395  }
396  nom = com_log(nom, alpha(s_prim, kk) + 0.5 * exp_temp0 + beta(s0, k));
397  den = com_log(den, alpha(s_prim, kk) + 0.5 * exp_temp1 + beta(s1, k));
398  }
399  extrinsic_output(kk) = nom - den;
400  }
401 
402 }
403 
404 void Rec_Syst_Conv_Code::log_decode_n2(const vec &rec_systematic, const vec &rec_parity,
405  const vec &extrinsic_input, vec &extrinsic_output, bool in_terminated, std::string metric)
406 {
407  if (metric == "TABLE") { // use the QLLR decoder; also see comment under log_decode()
408  QLLRvec rec_systematic_q = llrcalc.to_qllr(rec_systematic);
409  QLLRvec rec_parity_q = llrcalc.to_qllr(rec_parity);
410  QLLRvec extrinsic_input_q = llrcalc.to_qllr(extrinsic_input);
411  QLLRvec extrinsic_output_q(length(extrinsic_output));
412  log_decode_n2(rec_systematic_q, rec_parity_q, extrinsic_input_q,
413  extrinsic_output_q, in_terminated);
414  extrinsic_output = llrcalc.to_double(extrinsic_output_q);
415  return;
416  }
417 
418  // const double INF = 10e300; // replaced by DEFINE to be file-wide in scope
419  double nom, den, exp_temp0, exp_temp1, rp;
420  int k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
421  int ext_info_length = extrinsic_input.length();
422  ivec p0, p1;
423  double ex, norm;
424 
425  //Set the internal metric:
426  if (metric == "LOGMAX") { com_log = max; }
427  else if (metric == "LOGMAP") { com_log = log_add; }
428  else {
429  it_error("Rec_Syst_Conv_Code::log_decode_n2: Illegal metric parameter");
430  }
431 
432  alpha.set_size(Nstates, block_length + 1, false);
433  beta.set_size(Nstates, block_length + 1, false);
434  gamma.set_size(2*Nstates, block_length + 1, false);
435  extrinsic_output.set_size(ext_info_length, false);
436  //denom.set_size(block_length+1,false); for (k=0; k<=block_length; k++) { denom(k) = -infinity; }
437 
438  if (in_terminated) { terminated = true; }
439 
440  //Check that Lc = 1.0
441  it_assert(Lc == 1.0,
442  "Rec_Syst_Conv_Code::log_decode_n2: This function assumes that Lc = 1.0. Please use proper scaling of the input data");
443 
444  //Initiate alpha
445  for (s = 1; s < Nstates; s++) { alpha(s, 0) = -infinity; }
446  alpha(0, 0) = 0.0;
447 
448  //Calculate alpha and gamma going forward through the trellis
449  for (k = 1; k <= block_length; k++) {
450  kk = k - 1;
451  if (kk < ext_info_length) {
452  ex = 0.5 * (extrinsic_input(kk) + rec_systematic(kk));
453  }
454  else {
455  ex = 0.5 * rec_systematic(kk);
456  }
457  rp = 0.5 * rec_parity(kk);
458  for (s = 0; s < Nstates; s++) {
459  s_prim0 = rev_state_trans(s, 0);
460  s_prim1 = rev_state_trans(s, 1);
461  if (output_parity(s_prim0 , 0)) { exp_temp0 = -rp; }
462  else { exp_temp0 = rp; }
463  if (output_parity(s_prim1 , 1)) { exp_temp1 = -rp; }
464  else { exp_temp1 = rp; }
465  gamma(2*s_prim0 , k) = ex + exp_temp0;
466  gamma(2*s_prim1 + 1, k) = -ex + exp_temp1;
467  alpha(s, k) = com_log(alpha(s_prim0, kk) + gamma(2 * s_prim0 , k),
468  alpha(s_prim1, kk) + gamma(2 * s_prim1 + 1, k));
469  //denom(k) = com_log( alpha(s,k), denom(k) );
470  }
471  norm = alpha(0, k); //norm = denom(k);
472  for (l = 0; l < Nstates; l++) { alpha(l, k) -= norm; }
473  }
474 
475  //Initiate beta
476  if (terminated) {
477  for (s = 1; s < Nstates; s++) { beta(s, block_length) = -infinity; }
478  beta(0, block_length) = 0.0;
479  }
480  else {
481  beta.set_col(block_length, alpha.get_col(block_length));
482  }
483 
484  //Calculate beta going backward in the trellis
485  for (k = block_length; k >= 1; k--) {
486  kk = k - 1;
487  for (s_prim = 0; s_prim < Nstates; s_prim++) {
488  beta(s_prim, kk) = com_log(beta(state_trans(s_prim, 0), k) + gamma(2 * s_prim, k),
489  beta(state_trans(s_prim, 1), k) + gamma(2 * s_prim + 1, k));
490  }
491  norm = beta(0, k); //norm = denom(k);
492  for (l = 0; l < Nstates; l++) { beta(l, k) -= norm; }
493  }
494 
495  //Calculate extrinsic output for each bit
496  for (k = 1; k <= block_length; k++) {
497  kk = k - 1;
498  if (kk < ext_info_length) {
499  nom = -infinity;
500  den = -infinity;
501  rp = 0.5 * rec_parity(kk);
502  for (s_prim = 0; s_prim < Nstates; s_prim++) {
503  if (output_parity(s_prim , 0)) { exp_temp0 = -rp; }
504  else { exp_temp0 = rp; }
505  if (output_parity(s_prim , 1)) { exp_temp1 = -rp; }
506  else { exp_temp1 = rp; }
507  nom = com_log(nom, alpha(s_prim, kk) + exp_temp0 + beta(state_trans(s_prim, 0), k));
508  den = com_log(den, alpha(s_prim, kk) + exp_temp1 + beta(state_trans(s_prim, 1), k));
509  }
510  extrinsic_output(kk) = nom - den;
511  }
512  }
513 }
514 
515 // === Below new decoder functions by EGL, using QLLR arithmetics ===========
516 
517 void Rec_Syst_Conv_Code::log_decode(const QLLRvec &rec_systematic, const QLLRmat &rec_parity,
518  const QLLRvec &extrinsic_input,
519  QLLRvec &extrinsic_output, bool in_terminated)
520 {
521 
522  int nom, den, exp_temp0, exp_temp1, rp, temp0, temp1;
523  int i, j, s0, s1, k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
524  // ivec p0, p1;
525 
526  alpha_q.set_size(Nstates, block_length + 1, false);
527  beta_q.set_size(Nstates, block_length + 1, false);
528  gamma_q.set_size(2*Nstates, block_length + 1, false);
529  extrinsic_output.set_size(block_length, false);
530  denom_q.set_size(block_length + 1, false);
531  for (k = 0; k <= block_length; k++) { denom_q(k) = -QLLR_MAX; }
532 
533  if (in_terminated) { terminated = true; }
534 
535  //Check that Lc = 1.0
536  it_assert(Lc == 1.0,
537  "Rec_Syst_Conv_Code::log_decode: This function assumes that Lc = 1.0. Please use proper scaling of the input data");
538 
539  //Calculate gamma_q
540  for (k = 1; k <= block_length; k++) {
541  kk = k - 1;
542  for (s_prim = 0; s_prim < Nstates; s_prim++) {
543  exp_temp0 = 0;
544  exp_temp1 = 0;
545  for (j = 0; j < (n - 1); j++) {
546  rp = rec_parity(kk, j);
547  if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; }
548  else { exp_temp0 -= rp; }
549  if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; }
550  else { exp_temp1 -= rp; }
551  }
552  // right shift cannot be used due to implementation dependancy of how sign is handled?
553  gamma_q(2*s_prim + 0, k) = ((extrinsic_input(kk) + rec_systematic(kk)) + exp_temp0) / 2;
554  gamma_q(2*s_prim + 1, k) = - ((extrinsic_input(kk) + rec_systematic(kk)) - exp_temp1) / 2;
555  }
556  }
557 
558  //Initiate alpha_q
559  for (j = 1; j < Nstates; j++) { alpha_q(j, 0) = -QLLR_MAX; }
560  alpha_q(0, 0) = 0;
561 
562  //Calculate alpha_q, going forward through the trellis
563  for (k = 1; k <= block_length; k++) {
564  for (s = 0; s < Nstates; s++) {
565  s_prim0 = rev_state_trans(s, 0);
566  s_prim1 = rev_state_trans(s, 1);
567  temp0 = alpha_q(s_prim0, k - 1) + gamma_q(2 * s_prim0 + 0, k);
568  temp1 = alpha_q(s_prim1, k - 1) + gamma_q(2 * s_prim1 + 1, k);
569  // alpha_q(s,k) = com_log( temp0, temp1 );
570  // denom_q(k) = com_log( alpha_q(s,k), denom_q(k) );
571  alpha_q(s, k) = llrcalc.jaclog(temp0, temp1);
572  denom_q(k) = llrcalc.jaclog(alpha_q(s, k), denom_q(k));
573  }
574  //Normalization of alpha_q
575  for (l = 0; l < Nstates; l++) { alpha_q(l, k) -= denom_q(k); }
576  }
577 
578  //Initiate beta_q
579  if (terminated) {
580  for (i = 1; i < Nstates; i++) { beta_q(i, block_length) = -QLLR_MAX; }
581  beta_q(0, block_length) = 0;
582  }
583  else {
584  beta_q.set_col(block_length, alpha_q.get_col(block_length));
585  }
586 
587  //Calculate beta_q going backward in the trellis
588  for (k = block_length; k >= 1; k--) {
589  for (s_prim = 0; s_prim < Nstates; s_prim++) {
590  s0 = state_trans(s_prim, 0);
591  s1 = state_trans(s_prim, 1);
592  // beta_q(s_prim,k-1) = com_log( beta_q(s0,k) + gamma_q(2*s_prim+0,k) , beta_q(s1,k) + gamma_q(2*s_prim+1,k) );
593  beta_q(s_prim, k - 1) = llrcalc.jaclog(beta_q(s0, k) + gamma_q(2 * s_prim + 0, k) , beta_q(s1, k) + gamma_q(2 * s_prim + 1, k));
594  }
595  //Normalization of beta_q
596  for (l = 0; l < Nstates; l++) { beta_q(l, k - 1) -= denom_q(k); }
597  }
598 
599  //Calculate extrinsic output for each bit
600  for (k = 1; k <= block_length; k++) {
601  kk = k - 1;
602  nom = -QLLR_MAX;
603  den = -QLLR_MAX;
604  for (s_prim = 0; s_prim < Nstates; s_prim++) {
605  s0 = state_trans(s_prim, 0);
606  s1 = state_trans(s_prim, 1);
607  exp_temp0 = 0;
608  exp_temp1 = 0;
609  for (j = 0; j < (n - 1); j++) {
610  rp = rec_parity(kk, j);
611  if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; }
612  else { exp_temp0 -= rp; }
613  if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; }
614  else { exp_temp1 -= rp; }
615  }
616  // nom = com_log(nom, alpha_q(s_prim,kk) + 0.5*exp_temp0 + beta_q(s0,k) );
617  // den = com_log(den, alpha_q(s_prim,kk) + 0.5*exp_temp1 + beta_q(s1,k) );
618  nom = llrcalc.jaclog(nom, alpha_q(s_prim, kk) + exp_temp0 / 2 + beta_q(s0, k));
619  den = llrcalc.jaclog(den, alpha_q(s_prim, kk) + exp_temp1 / 2 + beta_q(s1, k));
620  }
621  extrinsic_output(kk) = nom - den;
622  }
623 
624 }
625 
626 
627 
628 void Rec_Syst_Conv_Code::log_decode_n2(const QLLRvec &rec_systematic,
629  const QLLRvec &rec_parity,
630  const QLLRvec &extrinsic_input,
631  QLLRvec &extrinsic_output,
632  bool in_terminated)
633 {
634  int nom, den, exp_temp0, exp_temp1, rp;
635  int k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
636  int ext_info_length = extrinsic_input.length();
637  ivec p0, p1;
638  int ex, norm;
639 
640 
641  alpha_q.set_size(Nstates, block_length + 1, false);
642  beta_q.set_size(Nstates, block_length + 1, false);
643  gamma_q.set_size(2*Nstates, block_length + 1, false);
644  extrinsic_output.set_size(ext_info_length, false);
645  //denom.set_size(block_length+1,false); for (k=0; k<=block_length; k++) { denom(k) = -infinity; }
646 
647  if (in_terminated) { terminated = true; }
648 
649  //Check that Lc = 1.0
650  it_assert(Lc == 1.0,
651  "Rec_Syst_Conv_Code::log_decode_n2: This function assumes that Lc = 1.0. Please use proper scaling of the input data");
652 
653  //Initiate alpha
654  for (s = 1; s < Nstates; s++) { alpha_q(s, 0) = -QLLR_MAX; }
655  alpha_q(0, 0) = 0;
656 
657  //Calculate alpha and gamma going forward through the trellis
658  for (k = 1; k <= block_length; k++) {
659  kk = k - 1;
660  if (kk < ext_info_length) {
661  ex = (extrinsic_input(kk) + rec_systematic(kk)) / 2;
662  }
663  else {
664  ex = rec_systematic(kk) / 2;
665  }
666  rp = rec_parity(kk) / 2;
667  for (s = 0; s < Nstates; s++) {
668  s_prim0 = rev_state_trans(s, 0);
669  s_prim1 = rev_state_trans(s, 1);
670  if (output_parity(s_prim0 , 0)) { exp_temp0 = -rp; }
671  else { exp_temp0 = rp; }
672  if (output_parity(s_prim1 , 1)) { exp_temp1 = -rp; }
673  else { exp_temp1 = rp; }
674  gamma_q(2*s_prim0 , k) = ex + exp_temp0;
675  gamma_q(2*s_prim1 + 1, k) = -ex + exp_temp1;
676  alpha_q(s, k) = llrcalc.jaclog(alpha_q(s_prim0, kk) + gamma_q(2 * s_prim0 , k),
677  alpha_q(s_prim1, kk) + gamma_q(2 * s_prim1 + 1, k));
678  //denom(k) = com_log( alpha(s,k), denom(k) );
679  }
680  norm = alpha_q(0, k); //norm = denom(k);
681  for (l = 0; l < Nstates; l++) { alpha_q(l, k) -= norm; }
682  }
683 
684  //Initiate beta
685  if (terminated) {
686  for (s = 1; s < Nstates; s++) { beta_q(s, block_length) = -QLLR_MAX; }
687  beta_q(0, block_length) = 0;
688  }
689  else {
690  beta_q.set_col(block_length, alpha_q.get_col(block_length));
691  }
692 
693  //Calculate beta going backward in the trellis
694  for (k = block_length; k >= 1; k--) {
695  kk = k - 1;
696  for (s_prim = 0; s_prim < Nstates; s_prim++) {
697  beta_q(s_prim, kk) = llrcalc.jaclog(beta_q(state_trans(s_prim, 0), k) + gamma_q(2 * s_prim, k),
698  beta_q(state_trans(s_prim, 1), k) + gamma_q(2 * s_prim + 1, k));
699  }
700  norm = beta_q(0, k); //norm = denom(k);
701  for (l = 0; l < Nstates; l++) { beta_q(l, k) -= norm; }
702  }
703 
704  //Calculate extrinsic output for each bit
705  for (k = 1; k <= block_length; k++) {
706  kk = k - 1;
707  if (kk < ext_info_length) {
708  nom = -QLLR_MAX;
709  den = -QLLR_MAX;
710  rp = rec_parity(kk) / 2;
711  for (s_prim = 0; s_prim < Nstates; s_prim++) {
712  if (output_parity(s_prim , 0)) { exp_temp0 = -rp; }
713  else { exp_temp0 = rp; }
714  if (output_parity(s_prim , 1)) { exp_temp1 = -rp; }
715  else { exp_temp1 = rp; }
716  nom = llrcalc.jaclog(nom, alpha_q(s_prim, kk) + exp_temp0 + beta_q(state_trans(s_prim, 0), k));
717  den = llrcalc.jaclog(den, alpha_q(s_prim, kk) + exp_temp1 + beta_q(state_trans(s_prim, 1), k));
718  }
719  extrinsic_output(kk) = nom - den;
720  }
721  }
722 }
723 
725 {
726  llrcalc = in_llrcalc;
727 }
728 
729 
730 } // namespace itpp
SourceForge Logo

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