IT++ Logo
siso_nsc.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_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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4