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_rsctrellis(void) 00039 //generates 1/2 RSC trellis structure for binary symbols 00040 //the states are numbered from 0 00041 { 00042 int mem_len = gen.cols()-1; 00043 register int n,k,j; 00044 itpp::bin feedback,out; 00045 int buffer; 00046 00047 rsctrellis.numStates = (1<<mem_len); 00048 rsctrellis.prevStates = new int[2*rsctrellis.numStates]; 00049 rsctrellis.nextStates = new int[2*rsctrellis.numStates]; 00050 rsctrellis.PARout = new double[2*rsctrellis.numStates]; 00051 rsctrellis.fm = new itpp::bin[rsctrellis.numStates]; 00052 00053 itpp::bvec cases(mem_len); 00054 for (n=0; n<2; n++) 00055 { 00056 #pragma omp parallel for private(k, j, cases, feedback, out, buffer) 00057 for (k=0; k<rsctrellis.numStates; k++) 00058 { 00059 cases = itpp::dec2bin(mem_len, k); 00060 //feedback 00061 feedback = (itpp::bin)n; 00062 for (j=1; j<(mem_len+1); j++) 00063 { 00064 feedback ^= (gen(0,j)*cases[j-1]); 00065 } 00066 //out 00067 out = feedback*gen(1,0); 00068 for (j=1; j<(mem_len+1); j++) 00069 { 00070 out ^= (gen(1,j)*cases[j-1]); 00071 } 00072 rsctrellis.PARout[k+n*rsctrellis.numStates] = (out?1.0:0.0);//parity bit 00073 rsctrellis.fm[k] = itpp::bin(n)^out; 00074 //shift 00075 for (j=mem_len-1; j>0; j--) 00076 { 00077 cases[j] = cases[j-1]; 00078 } 00079 cases[0] = feedback; 00080 //next and previous state 00081 buffer = itpp::bin2dec(cases, true); 00082 rsctrellis.nextStates[k+n*rsctrellis.numStates] = buffer;//next state 00083 rsctrellis.prevStates[buffer+n*rsctrellis.numStates] = k;//previous state 00084 } 00085 } 00086 } 00087 00088 void SISO::rsc_logMAP(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data, 00089 const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data) 00090 /* 00091 logMAP (SISO) decoder for RSC of rate 1/2 00092 extrinsic_coded - extrinsic information of coded bits 00093 extrinsic_data - extrinsic information of data (informational) bits 00094 intrinsic_coded - intrinsic information of coded (systematic and parity) bits 00095 apriori_data - a priori information of data (informational) bits 00096 Reference: Steven S. Pietrobon and Adrian S. Barbulescu, "A simplification of 00097 the modified Bahl decoding algorithm for systematic convolutional codes", Proc. ISITA, 1994 00098 */ 00099 { 00100 //get parameters 00101 int N = apriori_data.length(); 00102 //other parameters 00103 register int n,k; 00104 double buffer; 00105 int index; 00106 double A_min, A_max; 00107 double sum0, sum1; 00108 00109 //trellis generation 00110 gen_rsctrellis(); 00111 00112 //parameter initialization 00113 double* Lc1I = new double[N]; 00114 double* Lc2I = new double[N]; 00115 #pragma omp parallel for private(n) 00116 for (n=0; n<N; n++) 00117 { 00118 Lc1I[n] = intrinsic_coded[2*n]; 00119 Lc2I[n] = intrinsic_coded[2*n+1]; 00120 } 00121 double* A0 = new double[rsctrellis.numStates*N]; 00122 double* A1 = new double[rsctrellis.numStates*N]; 00123 double* A_mid = new double[N]; 00124 double* B0 = new double[rsctrellis.numStates*N]; 00125 double* B1 = new double[rsctrellis.numStates*N]; 00126 buffer = (tail?-INFINITY:0);//log(buffer) 00127 #pragma omp parallel for private(n,k) 00128 for (n=0; n<N; n++) 00129 { 00130 for (k=0; k<rsctrellis.numStates; k++) 00131 { 00132 A0[k+n*rsctrellis.numStates] = -INFINITY; 00133 A1[k+n*rsctrellis.numStates] = -INFINITY; 00134 B0[k+n*rsctrellis.numStates] = buffer; 00135 B1[k+n*rsctrellis.numStates] = buffer; 00136 } 00137 A_mid[n] = 0; 00138 } 00139 00140 //A 00141 A0[0] = Lc2I[0]*rsctrellis.PARout[0];//i=0 00142 A1[0] = Lc1I[0]+apriori_data[0]+Lc2I[0]*rsctrellis.PARout[rsctrellis.numStates];//i=1 00143 for (n=1; n<N; n++) 00144 { 00145 A_min = INFINITY; 00146 A_max = 0; 00147 for (k=0; k<rsctrellis.numStates; k++) 00148 { 00149 buffer = itpp::log_add(A0[rsctrellis.prevStates[k]+(n-1)*rsctrellis.numStates], 00150 A1[rsctrellis.prevStates[k+rsctrellis.numStates]+(n-1)*rsctrellis.numStates]);//log(alpha0+alpha1) 00151 A0[k+rsctrellis.numStates*n] = Lc2I[n]*rsctrellis.PARout[k]+buffer; 00152 A1[k+rsctrellis.numStates*n] = Lc1I[n]+apriori_data[n]+ 00153 Lc2I[n]*rsctrellis.PARout[k+rsctrellis.numStates]+buffer; 00154 //find min A(:,n) 00155 A_min = std::min(A_min, A0[k+rsctrellis.numStates*n]); 00156 //find max A(:,n) 00157 A_max = std::max(A_max, A0[k+rsctrellis.numStates*n]); 00158 } 00159 //normalization 00160 A_mid[n] = (A_min+A_max)/2; 00161 if (std::isinf(A_mid[n])) 00162 { 00163 continue; 00164 } 00165 for (k=0; k<rsctrellis.numStates; k++) 00166 { 00167 A0[k+rsctrellis.numStates*n] -= A_mid[n]; 00168 A1[k+rsctrellis.numStates*n] -= A_mid[n]; 00169 } 00170 } 00171 00172 //B 00173 B0[rsctrellis.prevStates[0]+(N-1)*rsctrellis.numStates] = 0; 00174 B1[rsctrellis.prevStates[rsctrellis.numStates]+(N-1)*rsctrellis.numStates] = 0; 00175 for (n=N-2; n>=0; n--) 00176 { 00177 for (k=0; k<rsctrellis.numStates; k++) 00178 { 00179 index = rsctrellis.nextStates[k]; 00180 B0[k+rsctrellis.numStates*n] = itpp::log_add(B0[index+(n+1)*rsctrellis.numStates]+ 00181 Lc2I[n+1]*rsctrellis.PARout[index], 00182 B1[index+(n+1)*rsctrellis.numStates]+ 00183 Lc1I[n+1]+apriori_data[n+1]+Lc2I[n+1]*rsctrellis.PARout[index+rsctrellis.numStates]); 00184 index = rsctrellis.nextStates[k+rsctrellis.numStates]; 00185 B1[k+rsctrellis.numStates*n] = itpp::log_add(B0[index+(n+1)*rsctrellis.numStates]+ 00186 Lc2I[n+1]*rsctrellis.PARout[index], 00187 B1[index+(n+1)*rsctrellis.numStates]+ 00188 Lc1I[n+1]+apriori_data[n+1]+Lc2I[n+1]*rsctrellis.PARout[index+rsctrellis.numStates]); 00189 00190 } 00191 if (std::isinf(A_mid[n+1])) 00192 { 00193 continue; 00194 } 00195 for (k=0; k<rsctrellis.numStates; k++) 00196 { 00197 B0[k+rsctrellis.numStates*n] -= A_mid[n+1]; 00198 B1[k+rsctrellis.numStates*n] -= A_mid[n+1]; 00199 } 00200 } 00201 00202 //updated LLR for information bits 00203 extrinsic_data.set_size(N); 00204 extrinsic_coded.set_size(2*N); 00205 #pragma omp parallel for private(n, k, sum0, sum1) 00206 for (n=0; n<N; n++) 00207 { 00208 sum0 = 0; 00209 sum1 = 0; 00210 for (k=0; k<rsctrellis.numStates; k++) 00211 { 00212 sum1 += std::exp(A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]); 00213 sum0 += std::exp(A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]); 00214 } 00215 extrinsic_data[n] = std::log(sum1/sum0)-apriori_data[n];//updated information must be independent of input LLR 00216 extrinsic_coded[2*n] = std::log(sum1/sum0)-Lc1I[n];//this information is used in serial concatenations 00217 } 00218 00219 //updated LLR for coded (parity) bits 00220 #pragma omp parallel for private(n, k, sum0, sum1) 00221 for (n=0; n<N; n++) 00222 { 00223 sum0 = 0; 00224 sum1 = 0; 00225 for (k=0; k<rsctrellis.numStates; k++) 00226 { 00227 if (rsctrellis.fm[k]) 00228 { 00229 sum1 += std::exp(A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]); 00230 sum0 += std::exp(A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]); 00231 } 00232 else 00233 { 00234 sum0 += std::exp(A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]); 00235 sum1 += std::exp(A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]); 00236 } 00237 } 00238 extrinsic_coded[2*n+1] = std::log(sum0/sum1)-Lc2I[n];//updated information must be independent of input LLR 00239 } 00240 00241 //destroy trellis 00242 delete[] rsctrellis.prevStates; 00243 delete[] rsctrellis.nextStates; 00244 delete[] rsctrellis.PARout; 00245 delete[] rsctrellis.fm; 00246 //destroy MAP parameters 00247 delete[] Lc1I; 00248 delete[] Lc2I; 00249 delete[] A0; 00250 delete[] A1; 00251 delete[] A_mid; 00252 delete[] B0; 00253 delete[] B1; 00254 } 00255 00256 void SISO::rsc_maxlogMAP(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data, 00257 const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data) 00258 /* 00259 maxlogMAP (SISO) decoder for RSC of rate 1/2 00260 extrinsic_coded - extrinsic information of coded bits 00261 extrinsic_data - extrinsic information of data (informational) bits 00262 intrinsic_coded - intrinsic information of coded (systematic and parity) bits 00263 apriori_data - a priori information of data (informational) bits 00264 Reference: Steven S. Pietrobon and Adrian S. Barbulescu, "A simplification of 00265 the modified Bahl decoding algorithm for systematic convolutional codes", Proc. ISITA, 1994 00266 */ 00267 { 00268 //get parameters 00269 int N = apriori_data.length(); 00270 //other parameters 00271 register int n,k; 00272 double buffer; 00273 int index; 00274 double A_min, A_max; 00275 double sum0, sum1; 00276 00277 //trellis generation 00278 gen_rsctrellis(); 00279 00280 //parameter initialization 00281 double* Lc1I = new double[N]; 00282 double* Lc2I = new double[N]; 00283 #pragma omp parallel for private(n) 00284 for (n=0; n<N; n++) 00285 { 00286 Lc1I[n] = intrinsic_coded[2*n]; 00287 Lc2I[n] = intrinsic_coded[2*n+1]; 00288 } 00289 double* A0 = new double[rsctrellis.numStates*N]; 00290 double* A1 = new double[rsctrellis.numStates*N]; 00291 double* A_mid = new double[N]; 00292 double* B0 = new double[rsctrellis.numStates*N]; 00293 double* B1 = new double[rsctrellis.numStates*N]; 00294 buffer = (tail?-INFINITY:0);//log(buffer) 00295 #pragma omp parallel for private(n,k) 00296 for (n=0; n<N; n++) 00297 { 00298 for (k=0; k<rsctrellis.numStates; k++) 00299 { 00300 A0[k+n*rsctrellis.numStates] = -INFINITY; 00301 A1[k+n*rsctrellis.numStates] = -INFINITY; 00302 B0[k+n*rsctrellis.numStates] = buffer; 00303 B1[k+n*rsctrellis.numStates] = buffer; 00304 } 00305 A_mid[n] = 0; 00306 } 00307 00308 //A 00309 A0[0] = Lc2I[0]*rsctrellis.PARout[0];//i=0 00310 A1[0] = Lc1I[0]+apriori_data[0]+Lc2I[0]*rsctrellis.PARout[rsctrellis.numStates];//i=1 00311 for (n=1; n<N; n++) 00312 { 00313 A_min = INFINITY; 00314 A_max = 0; 00315 for (k=0; k<rsctrellis.numStates; k++) 00316 { 00317 buffer = std::max(A0[rsctrellis.prevStates[k]+(n-1)*rsctrellis.numStates], 00318 A1[rsctrellis.prevStates[k+rsctrellis.numStates]+(n-1)*rsctrellis.numStates]);//log(alpha0+alpha1) 00319 A0[k+rsctrellis.numStates*n] = Lc2I[n]*rsctrellis.PARout[k]+buffer; 00320 A1[k+rsctrellis.numStates*n] = Lc1I[n]+apriori_data[n]+ 00321 Lc2I[n]*rsctrellis.PARout[k+rsctrellis.numStates]+buffer; 00322 //find min A(:,n) 00323 A_min = std::min(A_min, A0[k+rsctrellis.numStates*n]); 00324 //find max A(:,n) 00325 A_max = std::max(A_max, A0[k+rsctrellis.numStates*n]); 00326 } 00327 //normalization 00328 A_mid[n] = (A_min+A_max)/2; 00329 if (std::isinf(A_mid[n])) 00330 continue; 00331 for (k=0; k<rsctrellis.numStates; k++) 00332 { 00333 A0[k+rsctrellis.numStates*n] -= A_mid[n]; 00334 A1[k+rsctrellis.numStates*n] -= A_mid[n]; 00335 } 00336 } 00337 //B 00338 B0[rsctrellis.prevStates[0]+(N-1)*rsctrellis.numStates] = 0; 00339 B1[rsctrellis.prevStates[rsctrellis.numStates]+(N-1)*rsctrellis.numStates] = 0; 00340 for (n=N-2; n>=0; n--) 00341 { 00342 for (k=0; k<rsctrellis.numStates; k++) 00343 { 00344 index = rsctrellis.nextStates[k]; 00345 B0[k+rsctrellis.numStates*n] = std::max(B0[index+(n+1)*rsctrellis.numStates]+ 00346 Lc2I[n+1]*rsctrellis.PARout[index], 00347 B1[index+(n+1)*rsctrellis.numStates]+ 00348 Lc1I[n+1]+apriori_data[n+1]+Lc2I[n+1]*rsctrellis.PARout[index+rsctrellis.numStates]); 00349 index = rsctrellis.nextStates[k+rsctrellis.numStates]; 00350 B1[k+rsctrellis.numStates*n] = std::max(B0[index+(n+1)*rsctrellis.numStates]+ 00351 Lc2I[n+1]*rsctrellis.PARout[index], 00352 B1[index+(n+1)*rsctrellis.numStates]+ 00353 Lc1I[n+1]+apriori_data[n+1]+Lc2I[n+1]*rsctrellis.PARout[index+rsctrellis.numStates]); 00354 00355 } 00356 if (std::isinf(A_mid[n+1])) 00357 continue; 00358 for (k=0; k<rsctrellis.numStates; k++) 00359 { 00360 B0[k+rsctrellis.numStates*n] -= A_mid[n+1]; 00361 B1[k+rsctrellis.numStates*n] -= A_mid[n+1]; 00362 } 00363 } 00364 00365 //updated LLR for information bits 00366 extrinsic_data.set_size(N); 00367 extrinsic_coded.set_size(2*N); 00368 #pragma omp parallel for private(n, k, sum0, sum1) 00369 for (n=0; n<N; n++) 00370 { 00371 sum0 = -INFINITY; 00372 sum1 = -INFINITY; 00373 for (k=0; k<rsctrellis.numStates; k++) 00374 { 00375 sum1 = std::max(sum1, A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]); 00376 sum0 = std::max(sum0, A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]); 00377 } 00378 extrinsic_data[n] = (sum1-sum0)-apriori_data[n];//updated information must be independent of input LLR 00379 extrinsic_coded[2*n] = (sum1-sum0)-Lc1I[n]; 00380 } 00381 00382 //updated LLR for coded (parity) bits 00383 #pragma omp parallel for private(n, k, sum0, sum1) 00384 for (n=0; n<N; n++) 00385 { 00386 sum0 = -INFINITY; 00387 sum1 = -INFINITY; 00388 for (k=0; k<rsctrellis.numStates; k++) 00389 { 00390 if (rsctrellis.fm[k]) 00391 { 00392 sum1 = std::max(sum1, A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]); 00393 sum0 = std::max(sum0, A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]); 00394 } 00395 else 00396 { 00397 sum0 = std::max(sum0, A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]); 00398 sum1 = std::max(sum1, A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]); 00399 } 00400 } 00401 extrinsic_coded[2*n+1] = (sum0-sum1)-Lc2I[n];//updated information must be independent of input LLR 00402 } 00403 //destroy trellis 00404 delete[] rsctrellis.prevStates; 00405 delete[] rsctrellis.nextStates; 00406 delete[] rsctrellis.PARout; 00407 delete[] rsctrellis.fm; 00408 //destroy MAP parameters 00409 delete[] Lc1I; 00410 delete[] Lc2I; 00411 delete[] A0; 00412 delete[] A1; 00413 delete[] A_mid; 00414 delete[] B0; 00415 delete[] B1; 00416 } 00417 00418 void SISO::rsc_sova(itpp::vec &extrinsic_data, const itpp::vec &intrinsic_coded, 00419 const itpp::vec &apriori_data, const int &win_len) 00420 /* Soft Output Viterbi Algorithm (SOVA) optimized for 1/N encoders 00421 * Output: extrinsic_data - extrinsic information of data bits 00422 * Inputs: intrinsic_coded - intrinsic information of data bits 00423 * apriori_data - a priori information of data bits 00424 * win_len - window length used to represent the code trellis 00425 * 00426 * The original version has been written by Adrian Bohdanowicz (2003). 00427 * It is assumed that the BPSK mapping is: 0 -> +1, 1 -> -1. 00428 * Changes have been made to adapt the code for RSC codes of rate 1/2 00429 * and for soft input informations. 00430 * Improved SOVA has been implemented using a scaling factor and threshold for 00431 * the reliability information (see Wang [2003]). Even so, PCCC performance 00432 * are close to the original SOVA. 00433 */ 00434 { 00435 //number of code outputs 00436 int nb_outputs = gen.rows(); 00437 00438 //setup internal variables based on RSC trellis 00439 register int i,j,s; 00440 gen_rsctrellis();//trellis generation for 1/2 RSC codes 00441 int nb_states = rsctrellis.numStates; 00442 itpp::Array<itpp::mat> bin_out(2);//contains code output for each initial state and code input 00443 itpp::imat next_state(nb_states,2);//next state in the trellis 00444 for (i=0; i<2; i++) 00445 { 00446 bin_out(i).set_size(nb_states, nb_outputs); 00447 for (j=0; j<nb_states; j++) 00448 { 00449 bin_out(i)(j,0) = double(i);//systematic bit 00450 bin_out(i)(j,1) = rsctrellis.PARout[j+i*nb_states];//parity bit 00451 next_state(j,i) = rsctrellis.nextStates[j+i*nb_states]; 00452 } 00453 } 00454 itpp::vec bin_inp("0 1");//binary code inputs 00455 00456 int len = apriori_data.length();//number of decoding steps (total) 00457 00458 //allocate memory for the trellis window 00459 itpp::mat metr(nb_states,win_len+1);//path metric buffer 00460 metr.zeros(); 00461 metr += -INFINITY; 00462 metr(0,0) = 0;//initial state => (0,0) 00463 itpp::mat surv(nb_states,win_len+1);//survivor state buffer 00464 surv.zeros(); 00465 itpp::mat inpt(nb_states,win_len+1);//survivor input buffer (dec. output) 00466 inpt.zeros(); 00467 itpp::mat diff(nb_states,win_len+1);//path metric difference 00468 diff.zeros(); 00469 itpp::mat comp(nb_states,win_len+1);//competitor state buffer 00470 comp.zeros(); 00471 itpp::mat inpc(nb_states,win_len+1);//competitor input buffer 00472 inpc.zeros(); 00473 //soft output (sign with reliability) 00474 itpp::vec sft(len); 00475 sft.zeros(); 00476 sft += INFINITY; 00477 00478 //decode all the bits 00479 int Cur,Nxt,nxt,sur,b,tmp,idx; 00480 itpp::vec buf(nb_outputs); 00481 double llb,mtr,dif,cmp,inc,srv,inp; 00482 itpp::vec bin(nb_outputs); 00483 itpp::ivec cyclic_buffer(win_len); 00484 extrinsic_data.set_size(len); 00485 for (i = 0; i < len; i++) 00486 { 00487 //indices + precalculations 00488 Cur = i%(win_len+1);//curr trellis (cycl. buf) position 00489 Nxt = (i+1)%(win_len+1);//next trellis (cycl. buf) position 00490 buf = intrinsic_coded(i*nb_outputs,(i+1)*nb_outputs-1);//intrinsic_info portion to be processed 00491 llb = apriori_data(i);//SOVA: apriori_info portion to be processed 00492 metr.set_col(Nxt, -INFINITY*itpp::ones(nb_states)); 00493 00494 //forward recursion 00495 for (s = 0; s<nb_states; s++) 00496 { 00497 for (j = 0; j<2; j++) 00498 { 00499 nxt = next_state(s,j);//state after transition 00500 bin = bin_out(j).get_row(s);//transition output (encoder) 00501 mtr = bin*buf+metr(s,Cur);//transition metric 00502 mtr += bin_inp(j)*llb;//SOVA 00503 00504 if (metr(nxt,Nxt) < mtr) 00505 { 00506 diff(nxt,Nxt) = mtr-metr(nxt,Nxt);//SOVA 00507 comp(nxt,Nxt) = surv(nxt,Nxt);//SOVA 00508 inpc(nxt,Nxt) = inpt(nxt,Nxt);//SOVA 00509 00510 metr(nxt,Nxt) = mtr;//store the metric 00511 surv(nxt,Nxt) = s;//store the survival state 00512 inpt(nxt,Nxt) = j;//store the survival input 00513 } 00514 else 00515 { 00516 dif = metr(nxt,Nxt)-mtr; 00517 if (dif <= diff(nxt,Nxt)) 00518 { 00519 diff(nxt,Nxt) = dif;//SOVA 00520 comp(nxt,Nxt) = s;//SOVA 00521 inpc(nxt,Nxt) = j;//SOVA 00522 } 00523 } 00524 } 00525 } 00526 00527 //trace backwards 00528 if (i < (win_len-1)) 00529 { 00530 continue; 00531 }//proceed if the buffer has been filled 00532 mtr = itpp::max(metr.get_col(Nxt), sur);//find the initial state (max metric) 00533 b = i;//temporary bit index 00534 for (j=0; j<win_len; j++)//indices in a 'cyclic buffer' operation 00535 { 00536 cyclic_buffer(j) = (Nxt-j)%(win_len+1); 00537 cyclic_buffer(j) = (cyclic_buffer(j)<0)?(cyclic_buffer(j)+win_len+1):cyclic_buffer(j); 00538 } 00539 00540 for (j=0; j<win_len; j++)//for all the bits in the buffer 00541 { 00542 inp = inpt(sur,cyclic_buffer(j));//current bit-decoder output (encoder input) 00543 extrinsic_data(b) = inp;//store the hard output 00544 00545 tmp = cyclic_buffer(j); 00546 cmp = comp(sur, tmp);//SOVA: competitor state (previous) 00547 inc = inpc(sur, tmp);//SOVA: competitor bit output 00548 dif = diff(sur, tmp);//SOVA: corresp. path metric difference 00549 srv = surv(sur, tmp);//SOVA: temporary survivor path state 00550 00551 for (s=j+1; s<=win_len; s++)//check all buffer bits srv and cmp paths 00552 { 00553 if (inp != inc) 00554 { 00555 idx = b-((s-1)-j);//calculate index: [enc.k*(b-(k-1)+j-1)+1:enc.k*(b-(k-1)+j)] 00556 sft(idx) = std::min(sft(idx), dif);//update LLRs for bits that are different 00557 } 00558 if (srv == cmp) 00559 { 00560 break; 00561 }//stop if surv and comp merge (no need to continue) 00562 if (s == win_len) 00563 { 00564 break; 00565 }//stop if the end (otherwise: error) 00566 tmp = cyclic_buffer(s); 00567 inp = inpt(int(srv), tmp);//previous surv bit 00568 inc = inpt(int(cmp), tmp);//previous comp bit 00569 srv = surv(int(srv), tmp);//previous surv state 00570 cmp = surv(int(cmp), tmp);//previous comp state 00571 } 00572 sur = int(surv(sur, cyclic_buffer(j)));//state for the previous surv bit 00573 b--;//update bit index 00574 } 00575 } 00576 00577 // provide soft output with +/- sign: 00578 sft = threshold(sft, SOVA_threshold); 00579 extrinsic_data = 00580 itpp::elem_mult((2.0*extrinsic_data-1.0), SOVA_scaling_factor*sft)-apriori_data; 00581 00582 //free trellis memory 00583 delete[] rsctrellis.prevStates; 00584 delete[] rsctrellis.nextStates; 00585 delete[] rsctrellis.PARout; 00586 delete[] rsctrellis.fm; 00587 } 00588 00589 void SISO::rsc_viterbi(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data, 00590 const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data, const int &win_len) 00591 /* Soft Input Soft Output module based on Viterbi algorithm 00592 * Output: extrinsic_data - extrinsic information of data bits 00593 * Inputs: intrinsic_coded - intrinsic information of data bits 00594 * apriori_data - a priori information of data bits 00595 * win_len - window length used to represent the code trellis 00596 * 00597 * The implemented algorithm follows M. Kerner and O. Amrani, ''Iterative Decoding 00598 * Using Optimum Soft Input - Hard Output Module`` (2009), in: 00599 * IEEE Transactions on Communications, 57:7(1881-1885) 00600 */ 00601 { 00602 //number of code outputs 00603 int nb_outputs = gen.rows(); 00604 00605 //setup internal variables based on RSC trellis 00606 register int i,j,s; 00607 gen_rsctrellis();//trellis generation for 1/2 RSC codes 00608 int nb_states = rsctrellis.numStates; 00609 itpp::Array<itpp::mat> bin_out(2);//contains code output for each initial state and code input 00610 itpp::imat next_state(nb_states,2);//next state in the trellis 00611 for (i=0; i<2; i++) 00612 { 00613 bin_out(i).set_size(nb_states, nb_outputs); 00614 for (j=0; j<nb_states; j++) 00615 { 00616 bin_out(i)(j,0) = double(i);//systematic bit 00617 bin_out(i)(j,1) = rsctrellis.PARout[j+i*nb_states];//parity bit 00618 next_state(j,i) = rsctrellis.nextStates[j+i*nb_states]; 00619 } 00620 } 00621 itpp::vec bin_inp("0 1");//binary code inputs 00622 00623 int len = apriori_data.length();//number of decoding steps (total) 00624 00625 //allocate memory for the trellis window 00626 itpp::mat metr(nb_states,win_len+1);//path metric buffer 00627 metr.zeros(); 00628 metr += -INFINITY; 00629 metr(0,0) = 0;//initial state => (0,0) 00630 itpp::mat surv(nb_states,win_len+1);//survivor state buffer 00631 surv.zeros(); 00632 itpp::mat inpt(nb_states,win_len+1);//survivor input bits buffer (dec. output) 00633 inpt.zeros(); 00634 itpp::mat parity(nb_states,win_len+1);//survivor parity bits buffer (dec. output) 00635 parity.zeros(); 00636 00637 //decode all the bits 00638 int Cur,Nxt,nxt,sur,b; 00639 itpp::vec buf(nb_outputs); 00640 double llb,mtr; 00641 itpp::vec bin(nb_outputs); 00642 itpp::ivec cyclic_buffer(win_len); 00643 extrinsic_coded.set_size(2*len);//initialize memory for output 00644 extrinsic_data.set_size(len); 00645 for (i = 0; i < len; i++) 00646 { 00647 //indices + precalculations 00648 Cur = i%(win_len+1);//curr trellis (cycl. buf) position 00649 Nxt = (i+1)%(win_len+1);//next trellis (cycl. buf) position 00650 buf = intrinsic_coded(i*nb_outputs,(i+1)*nb_outputs-1);//intrinsic_info portion to be processed 00651 llb = apriori_data(i);//SOVA: apriori_info portion to be processed 00652 metr.set_col(Nxt, -INFINITY*itpp::ones(nb_states)); 00653 00654 //forward recursion 00655 for (s = 0; s<nb_states; s++) 00656 { 00657 for (j = 0; j<2; j++) 00658 { 00659 nxt = next_state(s,j);//state after transition 00660 bin = bin_out(j).get_row(s);//transition output (encoder) 00661 mtr = bin*buf+metr(s,Cur);//transition metric 00662 mtr += bin_inp(j)*llb;//add a priori info contribution 00663 00664 if (metr(nxt,Nxt) < mtr) 00665 { 00666 metr(nxt,Nxt) = mtr;//store the metric 00667 surv(nxt,Nxt) = s;//store the survival state 00668 inpt(nxt,Nxt) = bin(0);//j;//store the survival input 00669 parity(nxt,Nxt) = bin(1);//store survival parity bit 00670 } 00671 } 00672 } 00673 00674 //trace backwards 00675 if (i < (win_len-1)) 00676 { 00677 continue; 00678 }//proceed if the buffer has been filled 00679 mtr = itpp::max(metr.get_col(Nxt), sur);//find the initial state (max metric) 00680 b = i;//temporary bit index 00681 for (j=0; j<win_len; j++)//indices in a 'cyclic buffer' operation 00682 { 00683 cyclic_buffer(j) = (Nxt-j)%(win_len+1); 00684 cyclic_buffer(j) = (cyclic_buffer(j)<0)?(cyclic_buffer(j)+win_len+1):cyclic_buffer(j); 00685 } 00686 00687 for (j=0; j<win_len; j++)//for all the bits in the buffer 00688 { 00689 extrinsic_data(b) = inpt(sur,cyclic_buffer(j));//current bit-decoder output (encoder input) 00690 extrinsic_coded(2*b) = extrinsic_data(b); 00691 extrinsic_coded(2*b+1) = parity(sur,cyclic_buffer(j)); 00692 sur = int(surv(sur, cyclic_buffer(j)));//state for the previous surv bit 00693 b--;//update bit index 00694 } 00695 } 00696 00697 //free trellis memory 00698 delete[] rsctrellis.prevStates; 00699 delete[] rsctrellis.nextStates; 00700 delete[] rsctrellis.PARout; 00701 delete[] rsctrellis.fm; 00702 00703 //output only the hard output if the flag is true 00704 if (Viterbi_hard_output_flag==true) 00705 { 00706 return; 00707 } 00708 00709 // provide soft output (extrinsic information) for data bits 00710 //compute flipped and non-flipped bits positions 00711 itpp::bvec tmp(len); 00712 for (i=0; i<len; i++) 00713 { 00714 tmp(i) = ((apriori_data(i)+intrinsic_coded(2*i))>=0)^(extrinsic_data(i)==1.0); 00715 } 00716 itpp::ivec *ptr_idx_matching = new itpp::ivec; 00717 itpp::ivec *ptr_idx_nonmatching = new itpp::ivec; 00718 *ptr_idx_matching = itpp::find(tmp==itpp::bin(0)); 00719 *ptr_idx_nonmatching = itpp::find(tmp==itpp::bin(1)); 00720 //Estimated Bit Error Rate 00721 int idx_nonmatching_len = ptr_idx_nonmatching->length(); 00722 double error_rate = double(idx_nonmatching_len)/double(len); 00723 //Logarithm of Likelihood Ratio 00724 double LLR; 00725 if (error_rate==0.0) 00726 { 00727 LLR = std::log(double(len)); 00728 } else if (error_rate==1.0) 00729 { 00730 LLR = -std::log(double(len)); 00731 } else 00732 { 00733 LLR = std::log((1.0-error_rate)/error_rate); 00734 } 00735 for (i=0; i<ptr_idx_matching->length(); i++) 00736 { 00737 extrinsic_data(ptr_idx_matching->get(i)) = Viterbi_scaling_factor[0]* 00738 (2.0*extrinsic_data(ptr_idx_matching->get(i))-1.0)*LLR; 00739 } 00740 for (i=0; i<idx_nonmatching_len; i++) 00741 { 00742 extrinsic_data(ptr_idx_nonmatching->get(i)) = Viterbi_scaling_factor[1]* 00743 (2.0*extrinsic_data(ptr_idx_nonmatching->get(i))-1.0)*LLR; 00744 } 00745 //extrinsic information 00746 extrinsic_data -= apriori_data; 00747 00748 //provide soft output for coded bits 00749 tmp.set_length(2*len); 00750 for (i=0; i<(2*len); i++) 00751 { 00752 tmp(i) = (intrinsic_coded(i)>=0)^(extrinsic_coded(i)==1.0); 00753 } 00754 delete ptr_idx_matching;//free old memory 00755 delete ptr_idx_nonmatching; 00756 ptr_idx_matching = new itpp::ivec;//allocate memory for new vectors 00757 ptr_idx_nonmatching = new itpp::ivec; 00758 *ptr_idx_matching = itpp::find(tmp==itpp::bin(0)); 00759 *ptr_idx_nonmatching = itpp::find(tmp==itpp::bin(1)); 00760 //Estimated Bit Error Rate 00761 idx_nonmatching_len = ptr_idx_nonmatching->length(); 00762 error_rate = double(idx_nonmatching_len)/double(2*len); 00763 //Logarithm of Likelihood Ratio 00764 if (error_rate==0.0) 00765 { 00766 LLR = std::log(double(2*len)); 00767 } else if (error_rate==1.0) 00768 { 00769 LLR = -std::log(double(2*len)); 00770 } else 00771 { 00772 LLR = std::log((1.0-error_rate)/error_rate); 00773 } 00774 for (i=0; i<ptr_idx_matching->length(); i++) 00775 { 00776 extrinsic_coded(ptr_idx_matching->get(i)) = Viterbi_scaling_factor[0]* 00777 (2.0*extrinsic_coded(ptr_idx_matching->get(i))-1.0)*LLR; 00778 } 00779 for (i=0; i<idx_nonmatching_len; i++) 00780 { 00781 extrinsic_coded(ptr_idx_nonmatching->get(i)) = Viterbi_scaling_factor[1]* 00782 (2.0*extrinsic_coded(ptr_idx_nonmatching->get(i))-1.0)*LLR; 00783 } 00784 delete ptr_idx_matching; 00785 delete ptr_idx_nonmatching; 00786 //extrinsic information 00787 extrinsic_coded -= intrinsic_coded; 00788 } 00789 00790 }
Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4