00001 00029 #include <itpp/comm/siso.h> 00030 #include <itpp/base/itcompat.h> 00031 #include <limits> 00032 #ifndef INFINITY 00033 #define INFINITY std::numeric_limits<double>::infinity() 00034 #endif 00035 00036 namespace itpp 00037 { 00038 void SISO::gen_nsctrellis(void) 00039 /* 00040 generate trellis for non systematic non recursive convolutional codes 00041 r - number of outputs of the CC 00042 mem_len - memory length of the CC 00043 */ 00044 { 00045 //get parameters 00046 int r = gen.rows(); 00047 int mem_len = gen.cols()-1; 00048 //other parameters 00049 register int n,k,j,p; 00050 itpp::bin inputs[] = {0,1}; 00051 int index; 00052 00053 nsctrellis.stateNb = (1<<mem_len); 00054 nsctrellis.output = new double[nsctrellis.stateNb*2*r]; 00055 nsctrellis.nextState = new int[nsctrellis.stateNb*2]; 00056 nsctrellis.prevState = new int[nsctrellis.stateNb*2]; 00057 nsctrellis.input = new int[nsctrellis.stateNb*2]; 00058 00059 itpp::bvec enc_mem(mem_len); 00060 itpp::bin out; 00061 #pragma omp parallel for private(n,k,j,p,enc_mem,out) 00062 for (n=0; n<nsctrellis.stateNb; n++) //initial state 00063 { 00064 enc_mem = itpp::dec2bin(mem_len, n); 00065 //output 00066 for (k=0; k<2; k++) 00067 { 00068 for (j=0; j<r; j++) 00069 { 00070 out = inputs[k]*gen(j,0); 00071 for (p=1; p<=mem_len; p++) 00072 { 00073 out ^= (enc_mem[p-1]*gen(j,p));//0 or 1 00074 } 00075 nsctrellis.output[n+k*nsctrellis.stateNb+j*nsctrellis.stateNb*2] = double(out); 00076 } 00077 } 00078 //next state 00079 for (j=(mem_len-1); j>0; j--) 00080 { 00081 enc_mem[j] = enc_mem[j-1]; 00082 } 00083 for (k=0; k<2; k++) 00084 { 00085 enc_mem[0] = inputs[k]; 00086 nsctrellis.nextState[n+k*nsctrellis.stateNb] = itpp::bin2dec(enc_mem, true); 00087 } 00088 } 00089 00090 #pragma omp parallel for private(n,k,j,index) 00091 for (j=0; j<nsctrellis.stateNb; j++) 00092 { 00093 index = 0; 00094 for (n=0; n<nsctrellis.stateNb; n++) 00095 { 00096 for (k=0; k<2; k++) 00097 { 00098 if (nsctrellis.nextState[n+k*nsctrellis.stateNb]==j) 00099 { 00100 nsctrellis.prevState[j+index*nsctrellis.stateNb] = n; 00101 nsctrellis.input[j+index*nsctrellis.stateNb] = k;//0 or 1 00102 index++; 00103 } 00104 } 00105 } 00106 } 00107 } 00108 00109 void SISO::nsc_logMAP(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data, const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data) 00110 /* 00111 * generalized decoder for NSC codes (after the NSC code a scrambler of pattern phi is used) using log MAP alg. 00112 * extrinsic_coded - extrinsic information of coded bits (output if the shaping filter) 00113 * extrinsic_data - extrinsic information of data bits 00114 * intrinsic_coded - intrinsic information of coded bits 00115 * apriori_data - a priori information of data bits 00116 */ 00117 { 00118 //get parameters 00119 int N = apriori_data.length(); 00120 int Nc = scrambler_pattern.length(); 00121 int r = gen.rows();//number of outputs of the CC 00122 //other parameters 00123 register int n,k,m,mp,j,i; 00124 int pstates[2]; 00125 int nstates[2]; 00126 int inputs[2]; 00127 double C[2];//log(gamma) 00128 double sum; 00129 double sumbis; 00130 int index; 00131 00132 //initialize trellis 00133 gen_nsctrellis(); 00134 //initialize log(alpha) and log(beta) 00135 double* A = new double[nsctrellis.stateNb*(N+1)]; 00136 double* B = new double[nsctrellis.stateNb*(N+1)]; 00137 A[0] = 0; 00138 B[N*nsctrellis.stateNb] = 0; 00139 sum = (tail?-INFINITY:0); 00140 #pragma omp parallel for private(n) 00141 for (n=1; n<nsctrellis.stateNb; n++) 00142 { 00143 A[n] = -INFINITY; 00144 B[n+N*nsctrellis.stateNb] = sum;//if tail==false the final state is not known 00145 } 00146 00147 #pragma omp parallel sections private(n,sum,m,k,i,j,C) 00148 { 00149 //forward recursion 00150 for (n=1; n<(N+1); n++) 00151 { 00152 sum = 0;//normalization factor 00153 for (m=0; m<nsctrellis.stateNb; m++) //final state 00154 { 00155 for (k=0; k<2; k++) 00156 { 00157 pstates[k] = nsctrellis.prevState[m+nsctrellis.stateNb*k];//determine previous states 00158 inputs[k] = nsctrellis.input[m+nsctrellis.stateNb*k];//determine input 00159 C[k] = (inputs[k]?apriori_data[n-1]:0);//compute log of gamma 00160 } 00161 for (i=0; i<2; i++)//for each C[i] 00162 { 00163 for (k=0; k<r; k++) 00164 { 00165 for (j=0; j<Nc; j++) 00166 { 00167 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]; 00168 } 00169 } 00170 } 00171 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]); 00172 sum += std::exp(A[m+n*nsctrellis.stateNb]); 00173 } 00174 //normalization 00175 sum = std::log(sum); 00176 if (std::isinf(sum)) 00177 { 00178 continue; 00179 } 00180 for (m=0; m<nsctrellis.stateNb; m++) 00181 { 00182 A[m+n*nsctrellis.stateNb] -= sum; 00183 } 00184 } 00185 00186 //backward recursion 00187 #pragma omp section 00188 for (n=N-1; n>=0; n--) 00189 { 00190 sum = 0;//normalisation factor 00191 for (m=0; m<nsctrellis.stateNb; m++) //initial state 00192 { 00193 for (k=0; k<2; k++) 00194 { 00195 nstates[k] = nsctrellis.nextState[m+k*nsctrellis.stateNb];//determine next states 00196 C[k] = (k?apriori_data[n]:0);//compute log of gamma 00197 } 00198 for (i=0; i<2; i++) 00199 { 00200 for (k=0; k<r; k++) 00201 { 00202 for (j=0; j<Nc; j++) 00203 { 00204 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]; 00205 } 00206 } 00207 } 00208 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]); 00209 sum += std::exp(B[m+n*nsctrellis.stateNb]); 00210 } 00211 //normalisation 00212 sum = std::log(sum); 00213 if (std::isinf(sum)) 00214 { 00215 continue; 00216 } 00217 for (m=0; m<nsctrellis.stateNb; m++) 00218 { 00219 B[m+n*nsctrellis.stateNb] -= sum; 00220 } 00221 } 00222 } 00223 00224 //compute extrinsic_data 00225 extrinsic_data.set_size(N); 00226 #pragma omp parallel for private(n,k,sum,sumbis) 00227 for (n=1; n<(N+1); n++) 00228 { 00229 sum = 0; 00230 sumbis = 0; 00231 for (k=0; k<(nsctrellis.stateNb/2); k++) 00232 { 00233 sum += std::exp(A[k+(nsctrellis.stateNb/2)+n*nsctrellis.stateNb]+B[k+(nsctrellis.stateNb/2)+n*nsctrellis.stateNb]);//nominator 00234 sumbis += std::exp(A[k+n*nsctrellis.stateNb]+B[k+n*nsctrellis.stateNb]);//denominator 00235 } 00236 extrinsic_data[n-1] = std::log(sum/sumbis)-apriori_data[n-1]; 00237 } 00238 00239 //compute extrinsic_coded 00240 double *sum0 = new double[r]; 00241 double *sum1 = new double[r]; 00242 extrinsic_coded.set_size(N*Nc*r); 00243 for (n=0; n<N; n++) 00244 { 00245 for (k=0; k<r; k++) 00246 { 00247 sum0[k] = 0; 00248 sum1[k] = 0; 00249 } 00250 for (mp=0; mp<nsctrellis.stateNb; mp++) //initial state 00251 { 00252 for (i=0; i<2; i++) 00253 { 00254 m = nsctrellis.nextState[mp+i*nsctrellis.stateNb];//final state 00255 //compute log of sigma 00256 index = (m>=(nsctrellis.stateNb/2));//0 if input is 0, 1 if input is 1 00257 C[0] = (index?apriori_data[n]:0); 00258 for (k=0; k<r; k++) 00259 { 00260 for (j=0; j<Nc; j++) 00261 { 00262 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]; 00263 } 00264 } 00265 C[1] = A[mp+n*nsctrellis.stateNb]+C[0]+B[m+(n+1)*nsctrellis.stateNb];//this is only a buffer 00266 //compute sums 00267 for (k=0; k<r; k++) 00268 { 00269 if (nsctrellis.output[mp+index*nsctrellis.stateNb+k*nsctrellis.stateNb*2]) 00270 { 00271 sum1[k] += std::exp(C[1]); 00272 } 00273 else 00274 { 00275 sum0[k] += std::exp(C[1]); 00276 } 00277 } 00278 } 00279 } 00280 for (k=0; k<r; k++) 00281 { 00282 for (j=0; j<Nc; j++) 00283 { 00284 index = j+k*Nc+n*r*Nc; 00285 extrinsic_coded[index] = (1-2*double(scrambler_pattern[j]))*std::log(sum1[k]/sum0[k])-intrinsic_coded[index]; 00286 } 00287 } 00288 } 00289 00290 //free memory 00291 delete[] nsctrellis.output; 00292 delete[] nsctrellis.nextState; 00293 delete[] nsctrellis.prevState; 00294 delete[] nsctrellis.input; 00295 delete[] A; 00296 delete[] B; 00297 delete[] sum0; 00298 delete[] sum1; 00299 } 00300 00301 void SISO::nsc_maxlogMAP(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data, const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data) 00302 /* 00303 * generalized decoder for NSC codes (after the NSC code a scrambler of pattern phi is used) using max log MAP alg. 00304 * extrinsic_coded - extrinsic information of coded bits (output if the shaping filter) 00305 * extrinsic_data - extrinsic information of data bits 00306 * intrinsic_coded - intrinsic information of coded bits 00307 * apriori_data - a priori information of data bits 00308 */ 00309 { 00310 //get parameters 00311 int N = apriori_data.length(); 00312 int Nc = scrambler_pattern.length(); 00313 int r = gen.rows();//number of outputs of the CC 00314 //other parameters 00315 register int n,k,m,mp,j,i; 00316 int pstates[2]; 00317 int nstates[2]; 00318 int inputs[2]; 00319 double C[2];//log(gamma) 00320 double sum; 00321 double sumbis; 00322 int index; 00323 00324 //initialize trellis 00325 gen_nsctrellis(); 00326 //initialize log(alpha) and log(beta) 00327 double* A = new double[nsctrellis.stateNb*(N+1)]; 00328 double* B = new double[nsctrellis.stateNb*(N+1)]; 00329 A[0] = 0; 00330 B[N*nsctrellis.stateNb] = 0; 00331 sum = (tail?-INFINITY:0); 00332 #pragma omp parallel for private(n) 00333 for (n=1; n<nsctrellis.stateNb; n++) 00334 { 00335 A[n] = -INFINITY; 00336 B[n+N*nsctrellis.stateNb] = sum;//if tail==false the final state is not known 00337 } 00338 00339 //forward recursion 00340 #pragma omp parallel sections private(n,sum,m,k,i,j,C) 00341 { 00342 for (n=1; n<(N+1); n++) 00343 { 00344 sum = -INFINITY;//normalisation factor 00345 for (m=0; m<nsctrellis.stateNb; m++) //final state 00346 { 00347 for (k=0; k<2; k++) 00348 { 00349 pstates[k] = nsctrellis.prevState[m+nsctrellis.stateNb*k];//determine previous states 00350 inputs[k] = nsctrellis.input[m+nsctrellis.stateNb*k];//determine input 00351 C[k] = (inputs[k]?apriori_data[n-1]:0);//compute log of gamma 00352 } 00353 for (i=0; i<2; i++)//for each C[i] 00354 { 00355 for (k=0; k<r; k++) 00356 { 00357 for (j=0; j<Nc; j++) 00358 { 00359 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]; 00360 } 00361 } 00362 } 00363 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]); 00364 sum = std::max(sum, A[m+n*nsctrellis.stateNb]); 00365 } 00366 //normalization 00367 if (std::isinf(sum)) 00368 { 00369 continue; 00370 } 00371 for (m=0; m<nsctrellis.stateNb; m++) 00372 { 00373 A[m+n*nsctrellis.stateNb] -= sum; 00374 } 00375 } 00376 00377 //backward recursion 00378 #pragma omp section 00379 for (n=N-1; n>=0; n--) 00380 { 00381 sum = -INFINITY;//normalisation factor 00382 for (m=0; m<nsctrellis.stateNb; m++) //initial state 00383 { 00384 for (k=0; k<2; k++) 00385 { 00386 nstates[k] = nsctrellis.nextState[m+k*nsctrellis.stateNb];//determine next states 00387 C[k] = (k?apriori_data[n]:0);//compute log of gamma 00388 } 00389 for (i=0; i<2; i++) 00390 { 00391 for (k=0; k<r; k++) 00392 { 00393 for (j=0; j<Nc; j++) 00394 { 00395 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]; 00396 } 00397 } 00398 } 00399 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]); 00400 sum = std::max(sum, B[m+n*nsctrellis.stateNb]); 00401 } 00402 //normalisation 00403 if (std::isinf(sum)) 00404 { 00405 continue; 00406 } 00407 for (m=0; m<nsctrellis.stateNb; m++) 00408 { 00409 B[m+n*nsctrellis.stateNb] -= sum; 00410 } 00411 } 00412 } 00413 00414 //compute extrinsic_data 00415 extrinsic_data.set_size(N); 00416 #pragma omp parallel for private(n,k,sum,sumbis) 00417 for (n=1; n<(N+1); n++) 00418 { 00419 sum = -INFINITY; 00420 sumbis = -INFINITY; 00421 for (k=0; k<(nsctrellis.stateNb/2); k++) 00422 { 00423 sum = std::max(sum, A[k+(nsctrellis.stateNb/2)+n*nsctrellis.stateNb]+B[k+(nsctrellis.stateNb/2)+n*nsctrellis.stateNb]);//nominator 00424 sumbis = std::max(sumbis, A[k+n*nsctrellis.stateNb]+B[k+n*nsctrellis.stateNb]);//denominator 00425 } 00426 extrinsic_data[n-1] = (sum-sumbis)-apriori_data[n-1]; 00427 } 00428 00429 //compute extrinsic_coded 00430 double *sum0 = new double[r]; 00431 double *sum1 = new double[r]; 00432 extrinsic_coded.set_size(N*Nc*r); 00433 for (n=0; n<N; n++) 00434 { 00435 for (k=0; k<r; k++) 00436 { 00437 sum0[k] = -INFINITY; 00438 sum1[k] = -INFINITY; 00439 } 00440 for (mp=0; mp<nsctrellis.stateNb; mp++) //initial state 00441 { 00442 for (i=0; i<2; i++) 00443 { 00444 m = nsctrellis.nextState[mp+i*nsctrellis.stateNb];//final state 00445 //compute log of sigma 00446 index = (m>=(nsctrellis.stateNb/2));//0 if input is 0, 1 if input is 1 00447 C[0] = (index?apriori_data[n]:0); 00448 for (k=0; k<r; k++) 00449 { 00450 for (j=0; j<Nc; j++) 00451 { 00452 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]; 00453 } 00454 } 00455 C[1] = A[mp+n*nsctrellis.stateNb]+C[0]+B[m+(n+1)*nsctrellis.stateNb];//this is only a buffer 00456 //compute sums 00457 for (k=0; k<r; k++) 00458 { 00459 if (nsctrellis.output[mp+index*nsctrellis.stateNb+k*nsctrellis.stateNb*2]) 00460 { 00461 sum1[k] = std::max(sum1[k], C[1]); 00462 } 00463 else 00464 { 00465 sum0[k] = std::max(sum0[k], C[1]); 00466 } 00467 } 00468 } 00469 } 00470 for (k=0; k<r; k++) 00471 { 00472 for (j=0; j<Nc; j++) 00473 { 00474 index = j+k*Nc+n*r*Nc; 00475 extrinsic_coded[index] = (1-2*double(scrambler_pattern[j]))*(sum1[k]-sum0[k])-intrinsic_coded[index]; 00476 } 00477 } 00478 } 00479 00480 //free memory 00481 delete[] nsctrellis.output; 00482 delete[] nsctrellis.nextState; 00483 delete[] nsctrellis.prevState; 00484 delete[] nsctrellis.input; 00485 delete[] A; 00486 delete[] B; 00487 delete[] sum0; 00488 delete[] sum1; 00489 } 00490 }//end namespace tr
Generated on Sat Jul 9 2011 15:21:32 for IT++ by Doxygen 1.7.4