IT++ Logo
siso_rsc.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:32 for IT++ by Doxygen 1.7.4