IT++ Logo
siso_nsc.cpp
Go to the documentation of this file.
1
29#include <itpp/comm/siso.h>
30#include <itpp/base/itcompat.h>
31#include <limits>
32#ifndef INFINITY
33#define INFINITY std::numeric_limits<double>::infinity()
34#endif
35
36namespace itpp
37{
38void SISO::gen_nsctrellis(void)
39/*
40 generate trellis for non systematic non recursive convolutional codes
41 r - number of outputs of the CC
42 mem_len - memory length of the CC
43 */
44{
45 //get parameters
46 int r = gen.rows();
47 int mem_len = gen.cols()-1;
48 //other parameters
49 register int n,k,j,p;
50 itpp::bin inputs[] = {0,1};
51 int index;
52
53 nsctrellis.stateNb = (1<<mem_len);
54 nsctrellis.output = new double[nsctrellis.stateNb*2*r];
55 nsctrellis.nextState = new int[nsctrellis.stateNb*2];
56 nsctrellis.prevState = new int[nsctrellis.stateNb*2];
57 nsctrellis.input = new int[nsctrellis.stateNb*2];
58
59 itpp::bvec enc_mem(mem_len);
60 itpp::bin out;
61#pragma omp parallel for private(n,k,j,p,enc_mem,out)
62 for (n=0; n<nsctrellis.stateNb; n++) //initial state
63 {
64 enc_mem = itpp::dec2bin(mem_len, n);
65 //output
66 for (k=0; k<2; k++)
67 {
68 for (j=0; j<r; j++)
69 {
70 out = inputs[k]*gen(j,0);
71 for (p=1; p<=mem_len; p++)
72 {
73 out ^= (enc_mem[p-1]*gen(j,p));//0 or 1
74 }
75 nsctrellis.output[n+k*nsctrellis.stateNb+j*nsctrellis.stateNb*2] = double(out);
76 }
77 }
78 //next state
79 for (j=(mem_len-1); j>0; j--)
80 {
81 enc_mem[j] = enc_mem[j-1];
82 }
83 for (k=0; k<2; k++)
84 {
85 enc_mem[0] = inputs[k];
86 nsctrellis.nextState[n+k*nsctrellis.stateNb] = itpp::bin2dec(enc_mem, true);
87 }
88 }
89
90#pragma omp parallel for private(n,k,j,index)
91 for (j=0; j<nsctrellis.stateNb; j++)
92 {
93 index = 0;
94 for (n=0; n<nsctrellis.stateNb; n++)
95 {
96 for (k=0; k<2; k++)
97 {
98 if (nsctrellis.nextState[n+k*nsctrellis.stateNb]==j)
99 {
100 nsctrellis.prevState[j+index*nsctrellis.stateNb] = n;
101 nsctrellis.input[j+index*nsctrellis.stateNb] = k;//0 or 1
102 index++;
103 }
104 }
105 }
106 }
107}
108
109void SISO::nsc_logMAP(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data, const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data)
110/*
111 * generalized decoder for NSC codes (after the NSC code a scrambler of pattern phi is used) using log MAP alg.
112 * extrinsic_coded - extrinsic information of coded bits (output if the shaping filter)
113 * extrinsic_data - extrinsic information of data bits
114 * intrinsic_coded - intrinsic information of coded bits
115 * apriori_data - a priori information of data bits
116 */
117{
118 //get parameters
119 int N = apriori_data.length();
120 int Nc = scrambler_pattern.length();
121 int r = gen.rows();//number of outputs of the CC
122 //other parameters
123 register int n,k,m,mp,j,i;
124 int pstates[2];
125 int nstates[2];
126 int inputs[2];
127 double C[2];//log(gamma)
128 double sum;
129 double sumbis;
130 int index;
131
132 //initialize trellis
133 gen_nsctrellis();
134 //initialize log(alpha) and log(beta)
135 double* A = new double[nsctrellis.stateNb*(N+1)];
136 double* B = new double[nsctrellis.stateNb*(N+1)];
137 A[0] = 0;
138 B[N*nsctrellis.stateNb] = 0;
139 sum = (tail?-INFINITY:0);
140#pragma omp parallel for private(n)
141 for (n=1; n<nsctrellis.stateNb; n++)
142 {
143 A[n] = -INFINITY;
144 B[n+N*nsctrellis.stateNb] = sum;//if tail==false the final state is not known
145 }
146
147#pragma omp parallel sections private(n,sum,m,k,i,j,C)
148 {
149 //forward recursion
150 for (n=1; n<(N+1); n++)
151 {
152 sum = 0;//normalization factor
153 for (m=0; m<nsctrellis.stateNb; m++) //final state
154 {
155 for (k=0; k<2; k++)
156 {
157 pstates[k] = nsctrellis.prevState[m+nsctrellis.stateNb*k];//determine previous states
158 inputs[k] = nsctrellis.input[m+nsctrellis.stateNb*k];//determine input
159 C[k] = (inputs[k]?apriori_data[n-1]:0);//compute log of gamma
160 }
161 for (i=0; i<2; i++)//for each C[i]
162 {
163 for (k=0; k<r; k++)
164 {
165 for (j=0; j<Nc; j++)
166 {
167 C[i] += nsctrellis.output[pstates[i]+inputs[i]*nsctrellis.stateNb+k*nsctrellis.stateNb*2]*(1-2*double(scrambler_pattern[j]))*intrinsic_coded[j+k*Nc+(n-1)*Nc*r];
168 }
169 }
170 }
171 A[m+n*nsctrellis.stateNb] = itpp::log_add(A[pstates[0]+(n-1)*nsctrellis.stateNb]+C[0], A[pstates[1]+(n-1)*nsctrellis.stateNb]+C[1]);
172 sum += std::exp(A[m+n*nsctrellis.stateNb]);
173 }
174 //normalization
175 sum = std::log(sum);
176 if (std::isinf(sum))
177 {
178 continue;
179 }
180 for (m=0; m<nsctrellis.stateNb; m++)
181 {
182 A[m+n*nsctrellis.stateNb] -= sum;
183 }
184 }
185
186 //backward recursion
187#pragma omp section
188 for (n=N-1; n>=0; n--)
189 {
190 sum = 0;//normalisation factor
191 for (m=0; m<nsctrellis.stateNb; m++) //initial state
192 {
193 for (k=0; k<2; k++)
194 {
195 nstates[k] = nsctrellis.nextState[m+k*nsctrellis.stateNb];//determine next states
196 C[k] = (k?apriori_data[n]:0);//compute log of gamma
197 }
198 for (i=0; i<2; i++)
199 {
200 for (k=0; k<r; k++)
201 {
202 for (j=0; j<Nc; j++)
203 {
204 C[i] += nsctrellis.output[m+i*nsctrellis.stateNb+k*nsctrellis.stateNb*2]*(1-2*double(scrambler_pattern[j]))*intrinsic_coded[j+k*Nc+n*Nc*r];
205 }
206 }
207 }
208 B[m+n*nsctrellis.stateNb] = itpp::log_add(B[nstates[0]+(n+1)*nsctrellis.stateNb]+C[0], B[nstates[1]+(n+1)*nsctrellis.stateNb]+C[1]);
209 sum += std::exp(B[m+n*nsctrellis.stateNb]);
210 }
211 //normalisation
212 sum = std::log(sum);
213 if (std::isinf(sum))
214 {
215 continue;
216 }
217 for (m=0; m<nsctrellis.stateNb; m++)
218 {
219 B[m+n*nsctrellis.stateNb] -= sum;
220 }
221 }
222 }
223
224 //compute extrinsic_data
225 extrinsic_data.set_size(N);
226#pragma omp parallel for private(n,k,sum,sumbis)
227 for (n=1; n<(N+1); n++)
228 {
229 sum = 0;
230 sumbis = 0;
231 for (k=0; k<(nsctrellis.stateNb/2); k++)
232 {
233 sum += std::exp(A[k+(nsctrellis.stateNb/2)+n*nsctrellis.stateNb]+B[k+(nsctrellis.stateNb/2)+n*nsctrellis.stateNb]);//nominator
234 sumbis += std::exp(A[k+n*nsctrellis.stateNb]+B[k+n*nsctrellis.stateNb]);//denominator
235 }
236 extrinsic_data[n-1] = std::log(sum/sumbis)-apriori_data[n-1];
237 }
238
239 //compute extrinsic_coded
240 double *sum0 = new double[r];
241 double *sum1 = new double[r];
242 extrinsic_coded.set_size(N*Nc*r);
243 for (n=0; n<N; n++)
244 {
245 for (k=0; k<r; k++)
246 {
247 sum0[k] = 0;
248 sum1[k] = 0;
249 }
250 for (mp=0; mp<nsctrellis.stateNb; mp++) //initial state
251 {
252 for (i=0; i<2; i++)
253 {
254 m = nsctrellis.nextState[mp+i*nsctrellis.stateNb];//final state
255 //compute log of sigma
256 index = (m>=(nsctrellis.stateNb/2));//0 if input is 0, 1 if input is 1
257 C[0] = (index?apriori_data[n]:0);
258 for (k=0; k<r; k++)
259 {
260 for (j=0; j<Nc; j++)
261 {
262 C[0] += nsctrellis.output[mp+index*nsctrellis.stateNb+k*nsctrellis.stateNb*2]*(1-2*double(scrambler_pattern[j]))*intrinsic_coded[j+k*Nc+n*Nc*r];
263 }
264 }
265 C[1] = A[mp+n*nsctrellis.stateNb]+C[0]+B[m+(n+1)*nsctrellis.stateNb];//this is only a buffer
266 //compute sums
267 for (k=0; k<r; k++)
268 {
269 if (nsctrellis.output[mp+index*nsctrellis.stateNb+k*nsctrellis.stateNb*2])
270 {
271 sum1[k] += std::exp(C[1]);
272 }
273 else
274 {
275 sum0[k] += std::exp(C[1]);
276 }
277 }
278 }
279 }
280 for (k=0; k<r; k++)
281 {
282 for (j=0; j<Nc; j++)
283 {
284 index = j+k*Nc+n*r*Nc;
285 extrinsic_coded[index] = (1-2*double(scrambler_pattern[j]))*std::log(sum1[k]/sum0[k])-intrinsic_coded[index];
286 }
287 }
288 }
289
290 //free memory
291 delete[] nsctrellis.output;
292 delete[] nsctrellis.nextState;
293 delete[] nsctrellis.prevState;
294 delete[] nsctrellis.input;
295 delete[] A;
296 delete[] B;
297 delete[] sum0;
298 delete[] sum1;
299}
300
301void SISO::nsc_maxlogMAP(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data, const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data)
302/*
303 * generalized decoder for NSC codes (after the NSC code a scrambler of pattern phi is used) using max log MAP alg.
304 * extrinsic_coded - extrinsic information of coded bits (output if the shaping filter)
305 * extrinsic_data - extrinsic information of data bits
306 * intrinsic_coded - intrinsic information of coded bits
307 * apriori_data - a priori information of data bits
308 */
309{
310 //get parameters
311 int N = apriori_data.length();
312 int Nc = scrambler_pattern.length();
313 int r = gen.rows();//number of outputs of the CC
314 //other parameters
315 register int n,k,m,mp,j,i;
316 int pstates[2];
317 int nstates[2];
318 int inputs[2];
319 double C[2];//log(gamma)
320 double sum;
321 double sumbis;
322 int index;
323
324 //initialize trellis
325 gen_nsctrellis();
326 //initialize log(alpha) and log(beta)
327 double* A = new double[nsctrellis.stateNb*(N+1)];
328 double* B = new double[nsctrellis.stateNb*(N+1)];
329 A[0] = 0;
330 B[N*nsctrellis.stateNb] = 0;
331 sum = (tail?-INFINITY:0);
332#pragma omp parallel for private(n)
333 for (n=1; n<nsctrellis.stateNb; n++)
334 {
335 A[n] = -INFINITY;
336 B[n+N*nsctrellis.stateNb] = sum;//if tail==false the final state is not known
337 }
338
339 //forward recursion
340#pragma omp parallel sections private(n,sum,m,k,i,j,C)
341 {
342 for (n=1; n<(N+1); n++)
343 {
344 sum = -INFINITY;//normalisation factor
345 for (m=0; m<nsctrellis.stateNb; m++) //final state
346 {
347 for (k=0; k<2; k++)
348 {
349 pstates[k] = nsctrellis.prevState[m+nsctrellis.stateNb*k];//determine previous states
350 inputs[k] = nsctrellis.input[m+nsctrellis.stateNb*k];//determine input
351 C[k] = (inputs[k]?apriori_data[n-1]:0);//compute log of gamma
352 }
353 for (i=0; i<2; i++)//for each C[i]
354 {
355 for (k=0; k<r; k++)
356 {
357 for (j=0; j<Nc; j++)
358 {
359 C[i] += nsctrellis.output[pstates[i]+inputs[i]*nsctrellis.stateNb+k*nsctrellis.stateNb*2]*(1-2*double(scrambler_pattern[j]))*intrinsic_coded[j+k*Nc+(n-1)*Nc*r];
360 }
361 }
362 }
363 A[m+n*nsctrellis.stateNb] = std::max(A[pstates[0]+(n-1)*nsctrellis.stateNb]+C[0], A[pstates[1]+(n-1)*nsctrellis.stateNb]+C[1]);
364 sum = std::max(sum, A[m+n*nsctrellis.stateNb]);
365 }
366 //normalization
367 if (std::isinf(sum))
368 {
369 continue;
370 }
371 for (m=0; m<nsctrellis.stateNb; m++)
372 {
373 A[m+n*nsctrellis.stateNb] -= sum;
374 }
375 }
376
377 //backward recursion
378#pragma omp section
379 for (n=N-1; n>=0; n--)
380 {
381 sum = -INFINITY;//normalisation factor
382 for (m=0; m<nsctrellis.stateNb; m++) //initial state
383 {
384 for (k=0; k<2; k++)
385 {
386 nstates[k] = nsctrellis.nextState[m+k*nsctrellis.stateNb];//determine next states
387 C[k] = (k?apriori_data[n]:0);//compute log of gamma
388 }
389 for (i=0; i<2; i++)
390 {
391 for (k=0; k<r; k++)
392 {
393 for (j=0; j<Nc; j++)
394 {
395 C[i] += nsctrellis.output[m+i*nsctrellis.stateNb+k*nsctrellis.stateNb*2]*(1-2*double(scrambler_pattern[j]))*intrinsic_coded[j+k*Nc+n*Nc*r];
396 }
397 }
398 }
399 B[m+n*nsctrellis.stateNb] = std::max(B[nstates[0]+(n+1)*nsctrellis.stateNb]+C[0], B[nstates[1]+(n+1)*nsctrellis.stateNb]+C[1]);
400 sum = std::max(sum, B[m+n*nsctrellis.stateNb]);
401 }
402 //normalisation
403 if (std::isinf(sum))
404 {
405 continue;
406 }
407 for (m=0; m<nsctrellis.stateNb; m++)
408 {
409 B[m+n*nsctrellis.stateNb] -= sum;
410 }
411 }
412 }
413
414 //compute extrinsic_data
415 extrinsic_data.set_size(N);
416#pragma omp parallel for private(n,k,sum,sumbis)
417 for (n=1; n<(N+1); n++)
418 {
419 sum = -INFINITY;
420 sumbis = -INFINITY;
421 for (k=0; k<(nsctrellis.stateNb/2); k++)
422 {
423 sum = std::max(sum, A[k+(nsctrellis.stateNb/2)+n*nsctrellis.stateNb]+B[k+(nsctrellis.stateNb/2)+n*nsctrellis.stateNb]);//nominator
424 sumbis = std::max(sumbis, A[k+n*nsctrellis.stateNb]+B[k+n*nsctrellis.stateNb]);//denominator
425 }
426 extrinsic_data[n-1] = (sum-sumbis)-apriori_data[n-1];
427 }
428
429 //compute extrinsic_coded
430 double *sum0 = new double[r];
431 double *sum1 = new double[r];
432 extrinsic_coded.set_size(N*Nc*r);
433 for (n=0; n<N; n++)
434 {
435 for (k=0; k<r; k++)
436 {
437 sum0[k] = -INFINITY;
438 sum1[k] = -INFINITY;
439 }
440 for (mp=0; mp<nsctrellis.stateNb; mp++) //initial state
441 {
442 for (i=0; i<2; i++)
443 {
444 m = nsctrellis.nextState[mp+i*nsctrellis.stateNb];//final state
445 //compute log of sigma
446 index = (m>=(nsctrellis.stateNb/2));//0 if input is 0, 1 if input is 1
447 C[0] = (index?apriori_data[n]:0);
448 for (k=0; k<r; k++)
449 {
450 for (j=0; j<Nc; j++)
451 {
452 C[0] += nsctrellis.output[mp+index*nsctrellis.stateNb+k*nsctrellis.stateNb*2]*(1-2*double(scrambler_pattern[j]))*intrinsic_coded[j+k*Nc+n*Nc*r];
453 }
454 }
455 C[1] = A[mp+n*nsctrellis.stateNb]+C[0]+B[m+(n+1)*nsctrellis.stateNb];//this is only a buffer
456 //compute sums
457 for (k=0; k<r; k++)
458 {
459 if (nsctrellis.output[mp+index*nsctrellis.stateNb+k*nsctrellis.stateNb*2])
460 {
461 sum1[k] = std::max(sum1[k], C[1]);
462 }
463 else
464 {
465 sum0[k] = std::max(sum0[k], C[1]);
466 }
467 }
468 }
469 }
470 for (k=0; k<r; k++)
471 {
472 for (j=0; j<Nc; j++)
473 {
474 index = j+k*Nc+n*r*Nc;
475 extrinsic_coded[index] = (1-2*double(scrambler_pattern[j]))*(sum1[k]-sum0[k])-intrinsic_coded[index];
476 }
477 }
478 }
479
480 //free memory
481 delete[] nsctrellis.output;
482 delete[] nsctrellis.nextState;
483 delete[] nsctrellis.prevState;
484 delete[] nsctrellis.input;
485 delete[] A;
486 delete[] B;
487 delete[] sum0;
488 delete[] sum1;
489}
490}//end namespace tr
Binary arithmetic (boolean) class.
Definition binary.h:57
double log_add(double log_a, double log_b)
Safe substitute for log(exp(log_a) + exp(log_b))
Definition log_exp.cpp:40
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition matfunc.h:59
IT++ compatibility types and functions.
itpp namespace
Definition itmex.h:37
ITPP_EXPORT int bin2dec(const bvec &inbvec, bool msb_first=true)
Convert a bvec to decimal int with the first bit as MSB if msb_first == true.
ITPP_EXPORT bvec dec2bin(int length, int index)
Convert a decimal int index to bvec using length bits in the representation.
Definitions for Soft Input Soft Output (SISO) modules class.

Generated on Tue Aug 17 2021 10:59:15 for IT++ by Doxygen 1.12.0