IT++ Logo
siso_eq.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/comm/siso.h>
00030 #include <limits>
00031 #ifndef INFINITY
00032 #define INFINITY std::numeric_limits<double>::infinity()
00033 #endif
00034 
00035 namespace itpp
00036 {
00037 void SISO::gen_chtrellis(void)
00038 // generate trellis for precoded FIR channels with real coefficients
00039 {
00040     //get parameters
00041     int mem_len = impulse_response.cols()-1;//memory length of the channel
00042     int p_order = prec_gen.length()-1;//order of the precoder polynomial
00043 
00044     //other variables
00045     register int n,k,j;
00046     double inputs[] = {1.0,-1.0};//1->-1, 0->+1
00047     int index;
00048     double feedback[2];
00049 
00050     //create channel trellis
00051     int equiv_ch_mem_len = std::max(mem_len, p_order);
00052     chtrellis.stateNb = (1<<equiv_ch_mem_len);
00053     chtrellis.output = new double[2*chtrellis.stateNb];
00054     chtrellis.nextState = new int[2*chtrellis.stateNb];
00055     chtrellis.prevState = new int[2*chtrellis.stateNb];
00056     chtrellis.input = new int[2*chtrellis.stateNb];
00057 
00058     //initialize trellis
00059     itpp::ivec enc_mem(equiv_ch_mem_len);
00060 #pragma omp parallel for private(n,enc_mem,k,feedback,j)
00061     for (n=0; n<chtrellis.stateNb; n++) //initial state
00062     {
00063         enc_mem = 1-2*itpp::to_ivec(itpp::dec2bin(equiv_ch_mem_len, n));//1->-1, 0->+1
00064         //output
00065         for (k=0; k<2; k++)
00066         {
00067             feedback[k] = inputs[k];
00068             for (j=1; j<=p_order; j++)
00069                 if (prec_gen[j])
00070                     feedback[k] *= enc_mem[j-1];//xor truth table must remain the same
00071             chtrellis.output[n+k*chtrellis.stateNb] = feedback[k]*impulse_response(0,0);
00072             for (j=1; j<=mem_len; j++)
00073                 chtrellis.output[n+k*chtrellis.stateNb] += (enc_mem[j-1]*impulse_response(0,j));
00074         }
00075         //next state
00076         for (j=(equiv_ch_mem_len-1); j>0; j--)
00077             enc_mem[j] = enc_mem[j-1];
00078         for (k=0; k<2; k++)
00079         {
00080             enc_mem[0] = int(feedback[k]);
00081             chtrellis.nextState[n+k*chtrellis.stateNb] = itpp::bin2dec(itpp::to_bvec((1-enc_mem)/2), true);//-1->1, +1->0
00082         }
00083     }
00084 
00085 #pragma omp parallel for private(j,index,n,k)
00086     for (j=0; j<chtrellis.stateNb; j++)
00087     {
00088         index = 0;
00089         for (n=0; n<chtrellis.stateNb; n++)
00090         {
00091             for (k=0; k<2; k++)
00092             {
00093                 if (chtrellis.nextState[n+k*chtrellis.stateNb]==j)
00094                 {
00095                     chtrellis.prevState[j+index*chtrellis.stateNb] = n;
00096                     chtrellis.input[j+index*chtrellis.stateNb] = k;//this is an index to the channel input
00097                     index++;
00098                 }
00099             }
00100         }
00101     }
00102 }
00103 
00104 void SISO::equalizer_logMAP(itpp::vec &extrinsic_data, const itpp::vec &rec_sig, const itpp::vec &apriori_data)
00105 /*
00106   extrinsic_data - extrinsic information of channel input symbols
00107   rec_sig - received symbols
00108   apriori_data - a priori information of channel input symbols
00109  */
00110 {
00111     //get parameters
00112     int N = rec_sig.length();//length of the received frame
00113     //other parameters
00114     register int n,k,m;
00115     int pstates[2];
00116     int nstates[2];
00117     int inputs[2];
00118     double C[2];//log(gamma)
00119     double sum;
00120     double sumbis;
00121     double buffer;
00122 
00123     //initialize trellis
00124     gen_chtrellis();
00125     //initialize log(alpha) and log(beta)
00126     double* A = new double[chtrellis.stateNb*(N+1)];
00127     double* B = new double[chtrellis.stateNb*(N+1)];
00128     A[0] = 0;
00129     B[N*chtrellis.stateNb] = 0;
00130     sum = (tail?-INFINITY:0);
00131 #pragma omp parallel for private(n)
00132     for (n=1; n<chtrellis.stateNb; n++)
00133     {
00134         A[n] = -INFINITY;
00135         B[n+N*chtrellis.stateNb] = sum;//if tail==false the final state is not known
00136     }
00137 
00138 #pragma omp parallel sections private(n,sum,m,k,C)
00139     {
00140         //forward recursion
00141         for (n=1; n<=N; n++)
00142         {
00143             sum = 0;//normalisation factor
00144             for (m=0; m<chtrellis.stateNb; m++) //final state
00145             {
00146                 for (k=0; k<2; k++)
00147                 {
00148                     pstates[k] = chtrellis.prevState[m+chtrellis.stateNb*k];//determine previous states
00149                     inputs[k] = chtrellis.input[m+chtrellis.stateNb*k];//determine input
00150                     C[k] = (inputs[k])*apriori_data[n-1]-itpp::sqr(rec_sig[n-1]-chtrellis.output[pstates[k]+chtrellis.stateNb*inputs[k]])/(2*sigma2);//compute log of gamma
00151                 }
00152                 A[m+n*chtrellis.stateNb] = itpp::log_add(A[pstates[0]+(n-1)*chtrellis.stateNb]+C[0], A[pstates[1]+(n-1)*chtrellis.stateNb]+C[1]);
00153                 sum += std::exp(A[m+n*chtrellis.stateNb]);
00154             }
00155             //normalisation
00156             sum = std::log(sum);
00157             for (m=0; m<chtrellis.stateNb; m++)
00158             {
00159                 A[m+n*chtrellis.stateNb] -= sum;
00160             }
00161         }
00162 
00163         //backward recursion
00164 #pragma omp section
00165         for (n=N-1; n>=0; n--)
00166         {
00167             sum = 0;//normalisation factor
00168             for (m=0; m<chtrellis.stateNb; m++) //initial state
00169             {
00170                 for (k=0; k<2; k++)
00171                 {
00172                     nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states
00173                     C[k] = (k)*apriori_data[n]-itpp::sqr(rec_sig[n]-chtrellis.output[m+k*chtrellis.stateNb])/(2*sigma2);//compute log of gamma
00174                 }
00175                 B[m+n*chtrellis.stateNb] = itpp::log_add(B[nstates[0]+(n+1)*chtrellis.stateNb]+C[0], B[nstates[1]+(n+1)*chtrellis.stateNb]+C[1]);
00176                 sum += std::exp(B[m+n*chtrellis.stateNb]);
00177             }
00178             //normalisation
00179             sum = std::log(sum);
00180             for (m=0; m<chtrellis.stateNb; m++)
00181             {
00182                 B[m+n*chtrellis.stateNb] -= sum;
00183             }
00184         }
00185     }
00186 
00187     //compute extrinsic_data
00188     extrinsic_data.set_size(N);
00189 #pragma omp parallel for private(n,sum,sumbis,m,k,nstates,C,buffer)
00190     for (n=1; n<=N; n++)
00191     {
00192         sum = 0;//could be replaced by a vector
00193         sumbis = 0;
00194         for (m=0; m<chtrellis.stateNb; m++) //initial state
00195         {
00196             for (k=0; k<2; k++) //input index
00197             {
00198                 nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states
00199                 C[k] = (k)*apriori_data[n-1]-itpp::sqr(rec_sig[n-1]-chtrellis.output[m+k*chtrellis.stateNb])/(2*sigma2);//compute log of gamma
00200                 buffer = std::exp(A[m+(n-1)*chtrellis.stateNb]+C[k]+B[nstates[k]+n*chtrellis.stateNb]);
00201                 if (k)
00202                     sum += buffer;//1
00203                 else
00204                     sumbis += buffer;//0
00205             }
00206         }
00207         extrinsic_data[n-1] = std::log(sum/sumbis)-apriori_data[n-1];
00208     }
00209 
00210     //free memory
00211     delete[] chtrellis.output;
00212     delete[] chtrellis.nextState;
00213     delete[] chtrellis.prevState;
00214     delete[] chtrellis.input;
00215     delete[] A;
00216     delete[] B;
00217 }
00218 
00219 void SISO::equalizer_maxlogMAP(itpp::vec &extrinsic_data, const itpp::vec &rec_sig, const itpp::vec &apriori_data)
00220 /*
00221   extrinsic_data - extrinsic information for channel input symbols
00222   rec_sig - received symbols
00223   apriori_data - a priori information for channel input symbols
00224  */
00225 {
00226     //get parameters
00227     int N = rec_sig.length();//length of the received frame
00228     //other parameters
00229     register int n,k,m;
00230     int pstates[2];
00231     int nstates[2];
00232     int inputs[2];
00233     double C[2];//log(gamma)
00234     double sum;
00235     double sumbis;
00236     double buffer;
00237 
00238     //initialize trellis
00239     gen_chtrellis();
00240     //initialize log(alpha) and log(beta)
00241     double* A = new double[chtrellis.stateNb*(N+1)];
00242     double* B = new double[chtrellis.stateNb*(N+1)];
00243     A[0] = 0;
00244     for (n=1; n<chtrellis.stateNb; n++)
00245         A[n] = -INFINITY;
00246     B[N*chtrellis.stateNb] = 0;
00247     sum = (tail?-INFINITY:0);
00248     for (n=1; n<chtrellis.stateNb; n++)
00249         B[n+N*chtrellis.stateNb] = sum;//if tail==false the final state is not known
00250 
00251 #pragma omp parallel sections private(n,sum,m,k,C)
00252     {
00253         //forward recursion
00254         for (n=1; n<=N; n++)
00255         {
00256             sum = -INFINITY;//normalization factor
00257             for (m=0; m<chtrellis.stateNb; m++) //final state
00258             {
00259                 for (k=0; k<2; k++)
00260                 {
00261                     pstates[k] = chtrellis.prevState[m+chtrellis.stateNb*k];//determine previous states
00262                     inputs[k] = chtrellis.input[m+chtrellis.stateNb*k];//determine input
00263                     C[k] = (inputs[k])*apriori_data[n-1]-itpp::sqr(rec_sig[n-1]-chtrellis.output[pstates[k]+chtrellis.stateNb*inputs[k]])/(2*sigma2);//compute log of gamma
00264                 }
00265                 A[m+n*chtrellis.stateNb] = std::max(A[pstates[0]+(n-1)*chtrellis.stateNb]+C[0], A[pstates[1]+(n-1)*chtrellis.stateNb]+C[1]);
00266                 sum = std::max(sum, A[m+n*chtrellis.stateNb]);
00267             }
00268             for (m=0; m<chtrellis.stateNb; m++)
00269                 A[m+n*chtrellis.stateNb] -= sum;
00270         }
00271 
00272         //backward recursion
00273 #pragma omp section
00274         for (n=N-1; n>=0; n--)
00275         {
00276             sum = -INFINITY;//normalisation factor
00277             for (m=0; m<chtrellis.stateNb; m++) //initial state
00278             {
00279                 for (k=0; k<2; k++)
00280                 {
00281                     nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states
00282                     C[k] = (k)*apriori_data[n]-itpp::sqr(rec_sig[n]-chtrellis.output[m+k*chtrellis.stateNb])/(2*sigma2);//compute log of gamma
00283                 }
00284                 B[m+n*chtrellis.stateNb] = std::max(B[nstates[0]+(n+1)*chtrellis.stateNb]+C[0], B[nstates[1]+(n+1)*chtrellis.stateNb]+C[1]);
00285                 sum = std::max(sum, B[m+n*chtrellis.stateNb]);
00286             }
00287             for (m=0; m<chtrellis.stateNb; m++)
00288                 B[m+n*chtrellis.stateNb] -= sum;
00289         }
00290     }
00291 
00292     //compute extrinsic_data
00293     extrinsic_data.set_size(N);
00294     for (n=1; n<=N; n++)
00295     {
00296         sum = -INFINITY;
00297         sumbis = -INFINITY;
00298         for (m=0; m<chtrellis.stateNb; m++) //initial state
00299         {
00300             for (k=0; k<2; k++) //input index
00301             {
00302                 nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states
00303                 C[k] = (k)*apriori_data[n-1]-itpp::sqr(rec_sig[n-1]-chtrellis.output[m+k*chtrellis.stateNb])/(2*sigma2);//compute log of gamma
00304                 buffer = A[m+(n-1)*chtrellis.stateNb]+C[k]+B[nstates[k]+n*chtrellis.stateNb];
00305                 if (k)
00306                     sum = std::max(sum, buffer);//1
00307                 else
00308                     sumbis = std::max(sumbis, buffer);//0
00309             }
00310         }
00311         extrinsic_data[n-1] = (sum-sumbis)-apriori_data[n-1];
00312     }
00313 
00314     //free memory
00315     delete[] chtrellis.output;
00316     delete[] chtrellis.nextState;
00317     delete[] chtrellis.prevState;
00318     delete[] chtrellis.input;
00319     delete[] A;
00320     delete[] B;
00321 }
00322 }//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