IT++ 4.3.1
rec_syst_conv_code.cpp
Go to the documentation of this file.
1
28
30
31
32namespace itpp
33{
34
36double(*com_log)(double, double) = NULL;
37
39// This wrapper is because "com_log = std::max;" below caused an error
40inline double max(double x, double y) { return std::max(x, y); }
42
43// ----------------- Protected functions -----------------------------
44
45int 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 -------------------------
69void 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
124void 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
151void 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
166void 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 mat alpha(Nstates, block_length + 1);
174 mat beta(Nstates, block_length + 1);
175 mat gamma(2*Nstates, block_length + 1);
176 vec denom(block_length + 1);
177 denom.clear();
178
179 extrinsic_output.set_size(block_length, false);
180
181 if (in_terminated) { terminated = true; }
182
183 //Calculate gamma
184 for (k = 1; k <= block_length; k++) {
185 kk = k - 1;
186 for (s_prim = 0; s_prim < Nstates; s_prim++) {
187 exp_temp0 = 0.0;
188 exp_temp1 = 0.0;
189 for (j = 0; j < (n - 1); j++) {
190 exp_temp0 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 0)); /* Is this OK? */
191 exp_temp1 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 1));
192 }
193 // gamma(2*s_prim+0,k) = std::exp( 0.5*(extrinsic_input(kk) + Lc*rec_systematic(kk))) * std::exp( exp_temp0 );
194 // gamma(2*s_prim+1,k) = std::exp(-0.5*(extrinsic_input(kk) + Lc*rec_systematic(kk))) * std::exp( exp_temp1 );
195 /* == Changed to trunc_exp() to address bug 1088420 which
196 described a numerical instability problem in map_decode()
197 at high SNR. This should be regarded as a temporary fix and
198 it is not necessarily a waterproof one: multiplication of
199 probabilities still can result in values out of
200 range. (Range checking for the multiplication operator was
201 not implemented as it was felt that this would sacrifice
202 too much runtime efficiency. Some margin was added to the
203 numerical hardlimits below to reflect this. The hardlimit
204 values below were taken as the minimum range that a
205 "double" should support reduced by a few orders of
206 magnitude to make sure multiplication of several values
207 does not exceed the limits.)
208
209 It is suggested to use the QLLR based log-domain() decoders
210 instead of map_decode() as they are much faster and more
211 numerically stable.
212
213 EGL 8/06. == */
214 gamma(2*s_prim + 0, k) = trunc_exp(0.5 * (extrinsic_input(kk) + Lc * rec_systematic(kk)) + exp_temp0);
215 gamma(2*s_prim + 1, k) = trunc_exp(-0.5 * (extrinsic_input(kk) + Lc * rec_systematic(kk)) + exp_temp1);
216 }
217 }
218
219 //Initiate alpha
220 alpha.set_col(0, zeros(Nstates));
221 alpha(0, 0) = 1.0;
222
223 //Calculate alpha and denom going forward through the trellis
224 for (k = 1; k <= block_length; k++) {
225 for (s = 0; s < Nstates; s++) {
226 s_prim0 = rev_state_trans(s, 0);
227 s_prim1 = rev_state_trans(s, 1);
228 temp0 = alpha(s_prim0, k - 1) * gamma(2 * s_prim0 + 0, k);
229 temp1 = alpha(s_prim1, k - 1) * gamma(2 * s_prim1 + 1, k);
230 alpha(s, k) = temp0 + temp1;
231 denom(k) += temp0 + temp1;
232 }
233 alpha.set_col(k, alpha.get_col(k) / denom(k));
234 }
235
236 //Initiate beta
237 if (terminated) {
238 beta.set_col(block_length, zeros(Nstates));
239 beta(0, block_length) = 1.0;
240 }
241 else {
242 beta.set_col(block_length, alpha.get_col(block_length));
243 }
244
245 //Calculate beta going backward in the trellis
246 for (k = block_length; k >= 2; k--) {
247 for (s_prim = 0; s_prim < Nstates; s_prim++) {
248 s0 = state_trans(s_prim, 0);
249 s1 = state_trans(s_prim, 1);
250 beta(s_prim, k - 1) = (beta(s0, k) * gamma(2 * s_prim + 0, k)) + (beta(s1, k) * gamma(2 * s_prim + 1, k));
251 }
252 beta.set_col(k - 1, beta.get_col(k - 1) / denom(k));
253 }
254
255 //Calculate extrinsic output for each bit
256 for (k = 1; k <= block_length; k++) {
257 kk = k - 1;
258 nom = 0;
259 den = 0;
260 for (s_prim = 0; s_prim < Nstates; s_prim++) {
261 s0 = state_trans(s_prim, 0);
262 s1 = state_trans(s_prim, 1);
263 exp_temp0 = 0.0;
264 exp_temp1 = 0.0;
265 for (j = 0; j < (n - 1); j++) {
266 exp_temp0 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 0));
267 exp_temp1 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 1));
268 }
269 // gamma_k_e = std::exp( exp_temp0 );
270 gamma_k_e = trunc_exp(exp_temp0);
271 nom += alpha(s_prim, k - 1) * gamma_k_e * beta(s0, k);
272
273 // gamma_k_e = std::exp( exp_temp1 );
274 gamma_k_e = trunc_exp(exp_temp1);
275 den += alpha(s_prim, k - 1) * gamma_k_e * beta(s1, k);
276 }
277 // extrinsic_output(kk) = std::log(nom/den);
278 extrinsic_output(kk) = trunc_log(nom / den);
279 }
280
281}
282
283void Rec_Syst_Conv_Code::log_decode(const vec &rec_systematic, const mat &rec_parity,
284 const vec &extrinsic_input, vec &extrinsic_output, bool in_terminated, std::string metric)
285{
286 if (metric == "TABLE") {
287 /* Use the QLLR decoder. This can probably be done more
288 efficiently since it converts floating point vectors to QLLR.
289 However we have to live with this for the time being. */
290 QLLRvec rec_systematic_q = llrcalc.to_qllr(rec_systematic);
291 QLLRmat rec_parity_q = llrcalc.to_qllr(rec_parity);
292 QLLRvec extrinsic_input_q = llrcalc.to_qllr(extrinsic_input);
293 QLLRvec extrinsic_output_q(length(extrinsic_output));
294 log_decode(rec_systematic_q, rec_parity_q, extrinsic_input_q,
295 extrinsic_output_q, in_terminated);
296 extrinsic_output = llrcalc.to_double(extrinsic_output_q);
297 return;
298 }
299
300 double nom, den, exp_temp0, exp_temp1, rp, temp0, temp1;
301 int i, j, s0, s1, k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
302 ivec p0, p1;
303
304 //Set the internal metric:
305 if (metric == "LOGMAX") { com_log = max; }
306 else if (metric == "LOGMAP") { com_log = log_add; }
307 else {
308 it_error("Rec_Syst_Conv_Code::log_decode: Illegal metric parameter");
309 }
310
311 mat alpha(Nstates, block_length + 1);
312 mat beta(Nstates, block_length + 1);
313 mat gamma(2*Nstates, block_length + 1);
314 extrinsic_output.set_size(block_length, false);
315 vec denom(block_length + 1);
316 for (k = 0; k <= block_length; k++) { denom(k) = -infinity; }
317
318 if (in_terminated) { terminated = true; }
319
320 //Check that Lc = 1.0
321 it_assert(Lc == 1.0,
322 "Rec_Syst_Conv_Code::log_decode: This function assumes that Lc = 1.0. Please use proper scaling of the input data");
323
324 //Calculate gamma
325 for (k = 1; k <= block_length; k++) {
326 kk = k - 1;
327 for (s_prim = 0; s_prim < Nstates; s_prim++) {
328 exp_temp0 = 0.0;
329 exp_temp1 = 0.0;
330 for (j = 0; j < (n - 1); j++) {
331 rp = rec_parity(kk, j);
332 if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; }
333 else { exp_temp0 -= rp; }
334 if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; }
335 else { exp_temp1 -= rp; }
336 }
337 gamma(2*s_prim + 0, k) = 0.5 * ((extrinsic_input(kk) + rec_systematic(kk)) + exp_temp0);
338 gamma(2*s_prim + 1, k) = -0.5 * ((extrinsic_input(kk) + rec_systematic(kk)) - exp_temp1);
339 }
340 }
341
342 //Initiate alpha
343 for (j = 1; j < Nstates; j++) { alpha(j, 0) = -infinity; }
344 alpha(0, 0) = 0.0;
345
346 //Calculate alpha, going forward through the trellis
347 for (k = 1; k <= block_length; k++) {
348 for (s = 0; s < Nstates; s++) {
349 s_prim0 = rev_state_trans(s, 0);
350 s_prim1 = rev_state_trans(s, 1);
351 temp0 = alpha(s_prim0, k - 1) + gamma(2 * s_prim0 + 0, k);
352 temp1 = alpha(s_prim1, k - 1) + gamma(2 * s_prim1 + 1, k);
353 alpha(s, k) = com_log(temp0, temp1);
354 denom(k) = com_log(alpha(s, k), denom(k));
355 }
356 //Normalization of alpha
357 for (l = 0; l < Nstates; l++) { alpha(l, k) -= denom(k); }
358 }
359
360 //Initiate beta
361 if (terminated) {
362 for (i = 1; i < Nstates; i++) { beta(i, block_length) = -infinity; }
363 beta(0, block_length) = 0.0;
364 }
365 else {
366 beta.set_col(block_length, alpha.get_col(block_length));
367 }
368
369 //Calculate beta going backward in the trellis
370 for (k = block_length; k >= 1; k--) {
371 for (s_prim = 0; s_prim < Nstates; s_prim++) {
372 s0 = state_trans(s_prim, 0);
373 s1 = state_trans(s_prim, 1);
374 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));
375 }
376 //Normalization of beta
377 for (l = 0; l < Nstates; l++) { beta(l, k - 1) -= denom(k); }
378 }
379
380 //Calculate extrinsic output for each bit
381 for (k = 1; k <= block_length; k++) {
382 kk = k - 1;
383 nom = -infinity;
384 den = -infinity;
385 for (s_prim = 0; s_prim < Nstates; s_prim++) {
386 s0 = state_trans(s_prim, 0);
387 s1 = state_trans(s_prim, 1);
388 exp_temp0 = 0.0;
389 exp_temp1 = 0.0;
390 for (j = 0; j < (n - 1); j++) {
391 rp = rec_parity(kk, j);
392 if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; }
393 else { exp_temp0 -= rp; }
394 if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; }
395 else { exp_temp1 -= rp; }
396 }
397 nom = com_log(nom, alpha(s_prim, kk) + 0.5 * exp_temp0 + beta(s0, k));
398 den = com_log(den, alpha(s_prim, kk) + 0.5 * exp_temp1 + beta(s1, k));
399 }
400 extrinsic_output(kk) = nom - den;
401 }
402}
403
404void 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 mat alpha(Nstates, block_length + 1);
433 mat beta(Nstates, block_length + 1);
434 mat gamma(2*Nstates, block_length + 1);
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
517void 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 QLLRmat alpha_q(Nstates, block_length + 1);
527 QLLRmat beta_q(Nstates, block_length + 1);
528 QLLRmat gamma_q(2*Nstates, block_length + 1);
529 extrinsic_output.set_size(block_length, false);
530 QLLRvec denom_q(block_length + 1);
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
628void 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 QLLRmat alpha_q(Nstates, block_length + 1);
642 QLLRmat beta_q(Nstates, block_length + 1);
643 QLLRmat gamma_q(2*Nstates, block_length + 1);
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
Log-likelihood algebra calculation unit.
Definition llr.h:125
virtual void map_decode(const vec &rec_systematic, const mat &rec_parity, const vec &extrinsic_input, vec &extrinsic_output, bool set_terminated=false)
Maximum A Posteriori (MAP) Probability symbol-by-symbol decoder.
void set_generator_polynomials(const ivec &gen, int constraint_length)
Set generator polynomials.
void encode(const bvec &input, bmat &parity_bits)
Encode a binary vector of inputs starting from zero state without adding of a tail.
void encode_tail(const bvec &input, bvec &tail, bmat &parity_bits)
Encode a binary vector of inputs and also adds a tail of K-1 zeros to force the encoder into the zero...
virtual void log_decode_n2(const vec &rec_systematic, const vec &rec_parity, const vec &extrinsic_input, vec &extrinsic_output, bool set_terminated=false, std::string metric="LOGMAX")
Special Log-MAP/Log-MAX decoder implementation for n = 2.
void set_llrcalc(LLR_calc_unit in_llrcalc)
Set the lookup table for algebra with quantized LLR values (see LLR_calc_unit class)
void set_awgn_channel_parameters(double Ec, double N0)
Sets the channel parameters needed for MAP-decoding.
void set_scaling_factor(double in_Lc)
Set scaling factor for the decoder.
virtual void log_decode(const vec &rec_systematic, const mat &rec_parity, const vec &extrinsic_input, vec &extrinsic_output, bool set_terminated=false, std::string metric="LOGMAX")
Log domain implementation of the Maximum A Posteriori (MAP) Probability symbol-by-symbol decoder.
Binary arithmetic (boolean) class.
Definition binary.h:57
#define it_error(s)
Abort unconditionally.
Definition itassert.h:126
#define it_assert(t, s)
Abort if t is not true.
Definition itassert.h:94
double trunc_log(double x)
Truncated natural logarithm function.
Definition log_exp.h:115
double log_add(double log_a, double log_b)
Safe substitute for log(exp(log_a) + exp(log_b))
Definition log_exp.cpp:40
double trunc_exp(double x)
Truncated exponential function.
Definition log_exp.h:137
int length(const Vec< T > &v)
Length of vector.
Definition matfunc.h:51
double gamma(double x)
Deprecated gamma function - please use tgamma() instead.
Definition elem_math.cpp:79
T max(const Vec< T > &v)
Maximum value of vector.
Definition min_max.h:45
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
double norm(const cvec &v)
Calculate the 2-norm: norm(v)=sqrt(sum(abs(v).^2))
Definition misc_stat.cpp:77
Mat< bin > bmat
bin matrix
Definition mat.h:508
itpp namespace
Definition itmex.h:37
int reverse_int(int length, int in)
double(* com_log)(double, double)
Pointer to logarithmic branch metric function.
Definitions of a Recursive Systematic Convolutional codec class.