IT++ Logo
siso_dem.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::find_half_const(int &select_half, itpp::vec &re_part, itpp::bmat &re_bin_part, itpp::vec &im_part, itpp::bmat &im_bin_part)
00038 /* finds real in imaginary parts of the constellation and its corresponding bits
00039  * this approach is used for equivalent channel according to Hassibi's model
00040  * the constellation must be quadratic and the number of bits per symbol must be a multiple of two
00041  */
00042 {
00043     //values needed for initializations
00044     int const_size = itpp::pow2i(nb_bits_symb);//constellation size
00045     int half_nb_bits_symb = nb_bits_symb/2;
00046     int half_len = itpp::pow2i(half_nb_bits_symb);//number of values of real(imaginary) part
00047     //initialize output variables
00048     select_half = 0;
00049     re_part.set_size(half_len);
00050     re_bin_part.set_size(half_len, half_nb_bits_symb);
00051     re_part.zeros();
00052     re_part(0) = constellation(0).real();
00053     im_part.set_size(half_len);
00054     im_bin_part.set_size(half_len, half_nb_bits_symb);
00055     im_part.zeros();
00056     im_part(0) = constellation(0).imag();
00057     //select half for real (imaginary) to binary correspondence
00058     if (nb_bits_symb%2)
00059     {
00060         print_err_msg("SISO::find_half_const: number of bits per symbol must be a multiple of two");
00061         return;
00062     }
00063     const double min_diff = 1e-3;
00064     itpp::ivec idx = itpp::find(itpp::abs(itpp::real(constellation)-re_part(0))<min_diff);
00065     if (idx.length()!=half_len)
00066     {
00067         print_err_msg("SISO::find_half_const: the constellation must be quadratic");
00068         return;
00069     }
00070     itpp::bvec temp(nb_bits_symb);
00071     register int n;
00072     for (n=0; n<2; n++)
00073     {
00074         temp = bin_constellation.get_row(idx(n));
00075         re_bin_part.set_row(n,temp(0,half_nb_bits_symb-1));
00076     }
00077     select_half = (re_bin_part.get_row(0)==re_bin_part.get_row(1))?0:1;
00078     //algorithm
00079     double buffer;
00080     temp = bin_constellation.get_row(0);
00081     re_bin_part.set_row(0,temp(select_half*half_nb_bits_symb,(1+select_half)*half_nb_bits_symb-1));
00082     im_bin_part.set_row(0,temp((1-select_half)*half_nb_bits_symb,(2-select_half)*half_nb_bits_symb-1));
00083     int re_idx = 0;
00084     int im_idx = 0;
00085     for (n=1; n<const_size; n++)
00086     {
00087         temp = bin_constellation.get_row(n);
00088         buffer = constellation(n).real();
00089         if (itpp::prod(re_part-buffer)>min_diff)
00090         {
00091             re_idx++;
00092             re_part(re_idx) = buffer;
00093             re_bin_part.set_row(re_idx, temp(select_half*half_nb_bits_symb,(1+select_half)*half_nb_bits_symb-1));
00094         }
00095         buffer = constellation(n).imag();
00096         if (itpp::prod(im_part-buffer)>min_diff)
00097         {
00098             im_idx++;
00099             im_part(im_idx) = buffer;
00100             im_bin_part.set_row(im_idx, temp((1-select_half)*half_nb_bits_symb,(2-select_half)*half_nb_bits_symb-1));
00101         }
00102     }
00103 }
00104 
00105 void SISO::EquivRecSig(itpp::vec &x_eq, const itpp::cmat &rec_sig)
00106 //finds equivalent received signal with real coefficients
00107 //the equivalent received signal follows the model of Hassibi's paper
00108 //ouput:
00109 //x_eq - equivalent received signal with real coefficients
00110 //inputs:
00111 //rec_sig - received signal
00112 {
00113     for (int k=0; k<nb_rec_ant; k++)
00114     {
00115         x_eq.set_subvector(k*2*block_duration, itpp::real(rec_sig.get_col(k)));
00116         x_eq.set_subvector(k*2*block_duration+block_duration, itpp::imag(rec_sig.get_col(k)));
00117     }
00118 }
00119 
00120 void SISO::EquivCh(itpp::mat &H_eq, const itpp::cvec &H)
00121 //finds equivalent channel with real coefficients following the model of Hassibi's paper
00122 //output:
00123 //H_eq - equivalent channel
00124 //input:
00125 //H - channel matrix
00126 {
00127     itpp::mat Aq(2*block_duration,2*nb_em_ant);
00128     itpp::mat Bq(2*block_duration,2*nb_em_ant);
00129     itpp::cmat temp(block_duration,nb_em_ant);
00130     itpp::vec h(2*nb_em_ant);
00131     itpp::mat AhBh(2*block_duration,2);
00132     register int n,k;
00133     for (k=0; k<symbols_block; k++)
00134     {
00135         temp = ST_gen1.get(k*block_duration,k*block_duration+block_duration-1,0,nb_em_ant-1);
00136         Aq.set_submatrix(0, 0, itpp::real(temp));
00137         Aq.set_submatrix(0, nb_em_ant, -itpp::imag(temp));
00138         Aq.set_submatrix(block_duration, 0, itpp::imag(temp));
00139         Aq.set_submatrix(block_duration, nb_em_ant, itpp::real(temp));
00140         temp = ST_gen2.get(k*block_duration,k*block_duration+block_duration-1,0,nb_em_ant-1);
00141         Bq.set_submatrix(0, 0, -itpp::imag(temp));
00142         Bq.set_submatrix(0, nb_em_ant, -itpp::real(temp));
00143         Bq.set_submatrix(block_duration, 0, itpp::real(temp));
00144         Bq.set_submatrix(block_duration, nb_em_ant, -itpp::imag(temp));
00145         for (n=0; n<nb_rec_ant; n++)
00146         {
00147             h.set_subvector(0, real(H.mid(n*nb_em_ant,nb_em_ant)));
00148             h.set_subvector(nb_em_ant, imag(H.mid(n*nb_em_ant,nb_em_ant)));
00149             AhBh.set_col(0, Aq*h);
00150             AhBh.set_col(1, Bq*h);
00151             H_eq.set_submatrix(2*block_duration*n, 2*k, AhBh);
00152         }
00153     }
00154 }
00155 
00156 void SISO::Hassibi_maxlogMAP(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
00157 //maxlogMAP algorithm for ST block codes using Hassibi's model
00158 {
00159     //general parameters
00160     int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks of ST matrices/int period
00161     double N0 = 2*sigma2;//noise DSP
00162     int nb_all_symb = itpp::pow2i(nb_bits_symb*symbols_block);//nb. of all possible input symbols as a binary vector
00163     double nom,denom;//nominator and denominator of extrinsic information
00164     itpp::bvec bin_frame(nb_bits_symb*symbols_block);//binary frame at channel input
00165     itpp::bmat mat_bin_frame(nb_bits_symb, symbols_block);
00166     itpp::vec symb_frame_eq(2*symbols_block);//frame of symbols at equivalent channel input
00167     double temp;
00168     itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);//equivalent channel matrix
00169     itpp::vec x_eq(2*block_duration*nb_rec_ant);//equivalent received signal
00170     register int ns,q,nb,n,k;
00171     int index;
00172     extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
00173     //main loop
00174     for (ns=0; ns<nb_subblocks; ns++)//for each subblock
00175     {
00176         //find equivalent channel matrix
00177         EquivCh(H_eq, c_impulse_response.get_col(ns));
00178         //find equivalent received signal
00179         EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
00180         //compute the LLR of each bit in a frame of symbols_block symbols
00181         for (q=0; q<symbols_block; q++)//for each symbol in a subblock
00182         {
00183             for (nb=0; nb<nb_bits_symb; nb++)//for a given bit try all possible sollutions for the input symbol vector
00184             {
00185                 nom = -INFINITY;
00186                 denom = -INFINITY;
00187                 for (n=0; n<nb_all_symb; n++)//all possible symbols
00188                 {
00189                     bin_frame = itpp::dec2bin(nb_bits_symb*symbols_block, n);
00190                     mat_bin_frame = itpp::reshape(bin_frame, nb_bits_symb, symbols_block);
00191                     for (k=0; k<symbols_block; k++)
00192                     {
00193                         symb_frame_eq(2*k) = constellation(itpp::bin2dec(mat_bin_frame.get_col(k))).real();
00194                         symb_frame_eq(1+2*k) = constellation(itpp::bin2dec(mat_bin_frame.get_col(k))).imag();
00195                     }
00196                     temp = -itpp::sum_sqr(x_eq-H_eq*symb_frame_eq)/N0+\
00197                            itpp::to_vec(bin_frame)*apriori_data.mid(ns*nb_bits_symb*symbols_block,nb_bits_symb*symbols_block);
00198                     if (bin_frame(nb+q*nb_bits_symb))
00199                         nom = std::max(nom, temp);
00200                     else
00201                         denom = std::max(denom, temp);
00202                 }
00203                 index = nb+q*nb_bits_symb+ns*nb_bits_symb*symbols_block;
00204                 extrinsic_data(index) = (nom-denom)-apriori_data(index);
00205             }//bits/symbol
00206         }//symbols/subblock
00207     }//subblocks
00208 }
00209 
00210 void SISO::GA(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
00211 // Gaussian Approximation algorithm for ST codes using Hassibi's model
00212 {
00213     //general parameters
00214     int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks
00215     int half_nb_bits_symb = nb_bits_symb/2;
00216     int half_len = itpp::pow2i(half_nb_bits_symb);//number of values of real(imaginary) part
00217 
00218     //correspondence between real and imaginary part of symbols and their binary representations
00219     int select_half;
00220     itpp::vec re_part;
00221     itpp::bmat re_bin_part;
00222     itpp::vec im_part;
00223     itpp::bmat im_bin_part;
00224     find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
00225 
00226     //equivalent channel
00227     itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
00228     itpp::vec E_re_s(symbols_block);
00229     itpp::vec E_im_s(symbols_block);
00230     itpp::vec Var_re_s(symbols_block);
00231     itpp::vec Var_im_s(symbols_block);
00232     itpp::vec Ey(2*block_duration*nb_rec_ant);
00233     itpp::mat Cy(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
00234     itpp::mat Cy_inv(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
00235     itpp::vec x_eq(2*block_duration*nb_rec_ant);
00236     itpp::vec EZeta(2*block_duration*nb_rec_ant);
00237     itpp::mat CZeta_inv(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
00238     double nom,denom;
00239     double temp;
00240     register int ns,q,k,p,cs;
00241     int index;
00242     extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
00243     for (ns=0; ns<nb_subblocks; ns++)//subblock by subblock
00244     {
00245         //mean and variance of real and imaginary parts of emitted symbols
00246         E_re_s.zeros();
00247         E_im_s.zeros();
00248         Var_re_s.zeros();
00249         Var_im_s.zeros();
00250         for (q=0; q<symbols_block; q++)
00251         {
00252             index = q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
00253             for (k=0; k<half_len; k++)
00254             {
00255                 E_re_s(q) += re_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)),\
00256                                                    apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))),\
00257                                                    1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))));
00258                 E_im_s(q) += im_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(im_bin_part.get_row(k)),\
00259                                                    apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))),\
00260                                                    1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))));
00261                 Var_re_s(q) += itpp::sqr(re_part(k))*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)),\
00262                                apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))),\
00263                                1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))));
00264                 Var_im_s(q) += itpp::sqr(im_part(k))*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(im_bin_part.get_row(k)),\
00265                                apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))),\
00266                                1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))));
00267             }
00268             Var_re_s(q) -= itpp::sqr(E_re_s(q));
00269             Var_im_s(q) -= itpp::sqr(E_im_s(q));
00270         }
00271 
00272         //find equivalent channel
00273         EquivCh(H_eq, c_impulse_response.get_col(ns));
00274 
00275         //compute E[y] and Cov[y]
00276         Ey.zeros();
00277         Cy = sigma2*itpp::eye(2*block_duration*nb_rec_ant);
00278         for (q=0; q<symbols_block; q++)
00279         {
00280             //real & imaginary
00281             Ey += (H_eq.get_col(2*q)*E_re_s(q)+H_eq.get_col(1+2*q)*E_im_s(q));
00282             Cy += (itpp::outer_product(H_eq.get_col(2*q), H_eq.get_col(2*q)*Var_re_s(q))+\
00283                    itpp::outer_product(H_eq.get_col(1+2*q), H_eq.get_col(1+2*q)*Var_im_s(q)));
00284         }
00285 
00286         //inverse of Cov[y]
00287         Cy_inv = itpp::inv(Cy);
00288 
00289         //find equivalent received signal
00290         EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
00291 
00292         //compute extrinsic information of coded bits
00293         for (q=0; q<symbols_block; q++)
00294         {
00295             //real part
00296             EZeta = Ey-H_eq.get_col(2*q)*E_re_s(q);
00297             CZeta_inv = Cy_inv+itpp::outer_product(Cy_inv*\
00298                                                    ((Var_re_s(q)/(1-(((H_eq.get_col(2*q)).transpose()*Cy_inv)*(H_eq.get_col(2*q)*Var_re_s(q)))(0)))*\
00299                                                     H_eq.get_col(2*q)), Cy_inv.transpose()*H_eq.get_col(2*q));
00300             index = select_half*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
00301             for (p=0; p<half_nb_bits_symb; p++)
00302             {
00303                 nom = -INFINITY;
00304                 denom = -INFINITY;
00305                 for (cs=0; cs<half_len; cs++)
00306                 {
00307                     temp = -0.5*((x_eq-H_eq.get_col(2*q)*re_part(cs)-EZeta).transpose()*CZeta_inv*(x_eq-H_eq.get_col(2*q)*re_part(cs)-EZeta))(0)+\
00308                            itpp::to_vec(re_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
00309                     if (re_bin_part(cs,p))
00310                         nom = std::max(nom, temp);
00311                     else
00312                         denom = std::max(denom, temp);
00313                 }
00314                 extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
00315             }
00316             //imaginary part
00317             EZeta = Ey-H_eq.get_col(1+2*q)*E_im_s(q);
00318             CZeta_inv = Cy_inv+itpp::outer_product(Cy_inv*\
00319                                                    ((Var_im_s(q)/(1-(((H_eq.get_col(1+2*q)).transpose()*Cy_inv)*(H_eq.get_col(1+2*q)*Var_im_s(q)))(0)))*\
00320                                                     H_eq.get_col(1+2*q)), Cy_inv.transpose()*H_eq.get_col(1+2*q));
00321             index = (1-select_half)*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
00322             for (p=0; p<half_nb_bits_symb; p++)
00323             {
00324                 nom = -INFINITY;
00325                 denom = -INFINITY;
00326                 for (cs=0; cs<half_len; cs++)
00327                 {
00328                     temp = -0.5*((x_eq-H_eq.get_col(1+2*q)*im_part(cs)-EZeta).transpose()*CZeta_inv*(x_eq-H_eq.get_col(1+2*q)*im_part(cs)-EZeta))(0)+\
00329                            itpp::to_vec(im_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
00330                     if (im_bin_part(cs,p))
00331                         nom = std::max(nom, temp);
00332                     else
00333                         denom = std::max(denom, temp);
00334                 }
00335                 extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
00336             }
00337         }
00338     }//subblock by subblock
00339 }
00340 
00341 void SISO::sGA(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
00342 //simplified Gaussian Approximation algorithm for ST codes using Hassibi's model
00343 {
00344     //general parameters
00345     int nb_subblocks = (int)(rec_sig.rows()/block_duration);//number of subblocks
00346     int half_nb_bits_symb = (int)(nb_bits_symb/2);
00347     int half_len = itpp::pow2i(half_nb_bits_symb);//number of values of real(imaginary) part
00348 
00349     //correspondence between real and imaginary part of symbols and their binary representations
00350     int select_half;
00351     itpp::vec re_part;
00352     itpp::bmat re_bin_part;
00353     itpp::vec im_part;
00354     itpp::bmat im_bin_part;
00355     find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
00356 
00357     //equivalent channel
00358     itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
00359 
00360     itpp::vec E_re_s(symbols_block);
00361     itpp::vec E_im_s(symbols_block);
00362     itpp::vec Var_re_s(symbols_block);
00363     itpp::vec Var_im_s(symbols_block);
00364     itpp::vec Ey(2*block_duration*nb_rec_ant);
00365     itpp::mat Cy(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
00366     itpp::vec x_eq(2*block_duration*nb_rec_ant);
00367     itpp::vec EZeta(2*block_duration*nb_rec_ant);
00368     itpp::vec CZeta(2*block_duration*nb_rec_ant);
00369     double nom,denom;
00370     double temp;
00371     register int ns,q,k,p,cs;
00372     int index;
00373     extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
00374     for (ns=0; ns<nb_subblocks; ns++)//subblock by subblock
00375     {
00376         //mean and variance of real and imaginary parts of emitted symbols
00377         E_re_s.zeros();
00378         E_im_s.zeros();
00379         Var_re_s.zeros();
00380         Var_im_s.zeros();
00381         for (q=0; q<symbols_block; q++)
00382         {
00383             index = q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
00384             for (k=0; k<half_len; k++)
00385             {
00386                 E_re_s(q) += re_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)),\
00387                                                    apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))),\
00388                                                    1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))));
00389                 E_im_s(q) += im_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(im_bin_part.get_row(k)),\
00390                                                    apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))),\
00391                                                    1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))));
00392                 Var_re_s(q) += itpp::sqr(re_part(k))*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)),\
00393                                apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))),\
00394                                1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))));
00395                 Var_im_s(q) += itpp::sqr(im_part(k))*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(im_bin_part.get_row(k)),\
00396                                apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))),\
00397                                1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))));
00398             }
00399             Var_re_s(q) -= itpp::sqr(E_re_s(q));
00400             Var_im_s(q) -= itpp::sqr(E_im_s(q));
00401         }
00402 
00403         //find equivalent channel
00404         EquivCh(H_eq, c_impulse_response.get_col(ns));
00405 
00406         //compute E[y] and Cov[y]
00407         Ey.zeros();
00408         Cy = sigma2*itpp::eye(2*block_duration*nb_rec_ant);
00409         for (q=0; q<symbols_block; q++)
00410         {
00411             //real & imaginary
00412             Ey += (H_eq.get_col(2*q)*E_re_s(q)+H_eq.get_col(1+2*q)*E_im_s(q));
00413             Cy += (itpp::outer_product(H_eq.get_col(2*q), H_eq.get_col(2*q)*Var_re_s(q))+\
00414                    itpp::outer_product(H_eq.get_col(1+2*q), H_eq.get_col(1+2*q)*Var_im_s(q)));
00415         }
00416 
00417         //find equivalent received signal
00418         EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
00419 
00420         //compute extrinsic INFINITYormation of coded bits
00421         for (q=0; q<symbols_block; q++)
00422         {
00423             //real part
00424             EZeta = Ey-H_eq.get_col(2*q)*E_re_s(q);
00425             CZeta = diag(Cy-itpp::outer_product(H_eq.get_col(2*q), H_eq.get_col(2*q)*Var_re_s(q)));
00426             index = select_half*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
00427             for (p=0; p<half_nb_bits_symb; p++)
00428             {
00429                 nom = -INFINITY;
00430                 denom = -INFINITY;
00431                 for (cs=0; cs<half_len; cs++)
00432                 {
00433                     temp = -0.5*itpp::sum(itpp::elem_div(sqr(x_eq-H_eq.get_col(2*q)*re_part(cs)-EZeta), CZeta))+\
00434                            itpp::to_vec(re_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
00435                     if (re_bin_part(cs,p))
00436                         nom = std::max(nom, temp);
00437                     else
00438                         denom = std::max(denom, temp);
00439                 }
00440                 extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
00441             }
00442             //imaginary part
00443             EZeta = Ey-H_eq.get_col(1+2*q)*E_im_s(q);
00444             CZeta = itpp::diag(Cy-itpp::outer_product(H_eq.get_col(1+2*q), H_eq.get_col(1+2*q)*Var_im_s(q)));
00445             index = (1-select_half)*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
00446             for (p=0; p<half_nb_bits_symb; p++)
00447             {
00448                 nom = -INFINITY;
00449                 denom = -INFINITY;
00450                 for (cs=0; cs<half_len; cs++)
00451                 {
00452                     temp = -0.5*itpp::sum(itpp::elem_div(sqr(x_eq-H_eq.get_col(1+2*q)*im_part(cs)-EZeta), CZeta))+\
00453                            itpp::to_vec(im_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
00454                     if (im_bin_part(cs,p))
00455                         nom = std::max(nom, temp);
00456                     else
00457                         denom = std::max(denom, temp);
00458                 }
00459                 extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
00460             }
00461         }
00462     }//subblock by subblock
00463 }
00464 
00465 void SISO::mmsePIC(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
00466 //MMSE Parallel Interference Canceller
00467 {
00468     //general parameters
00469     int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks
00470     int half_nb_bits_symb = nb_bits_symb/2;
00471     int half_const_len = itpp::pow2i(half_nb_bits_symb);
00472     int nb_bits_subblock = nb_bits_symb*symbols_block;//number of coded bits in an ST block
00473     itpp::vec Es(2*symbols_block);
00474     itpp::vec Vs(2*symbols_block);
00475     itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
00476     itpp::mat K(2*nb_rec_ant*block_duration,2*nb_rec_ant*block_duration);
00477     itpp::mat K_inv(2*nb_rec_ant*block_duration,2*nb_rec_ant*block_duration);
00478     itpp::vec x_eq(2*nb_rec_ant*block_duration);
00479     itpp::vec interf(2*symbols_block);
00480     itpp::vec temp(2*nb_rec_ant*block_duration);
00481     itpp::vec w(2*nb_rec_ant*block_duration);//filter impulse response
00482     double s_tilde;
00483     double mu_res;
00484     double sigma2_res;
00485     double nom,denom;
00486     double tmp;
00487     register int ns,q,k,s;
00488     int index;
00489 
00490     //correspondence between real and imaginary part of symbols and their binary representations
00491     int select_half;
00492     itpp::vec re_part;
00493     itpp::bmat re_bin_part;
00494     itpp::vec im_part;
00495     itpp::bmat im_bin_part;
00496     find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
00497     double part_var = 1/(double)(2*nb_em_ant);//real and imaginary part variance
00498 
00499     extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
00500     for (ns=0; ns<nb_subblocks; ns++)//compute block by block
00501     {
00502         //mean and variance of real and imaginary parts of emitted symbols
00503         Es.zeros();
00504         Vs.zeros();
00505         for (q=0; q<symbols_block; q++)
00506         {
00507             index = q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
00508             for (k=0; k<half_const_len; k++)
00509             {
00510                 Es(2*q) += re_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)), \
00511                                                  apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb))), \
00512                                                  (1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb)))));
00513                 Es(1+2*q) += im_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(im_bin_part.get_row(k)), \
00514                                                    apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb))), \
00515                                                    (1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb)))));
00516                 Vs(2*q) += itpp::sqr(re_part(k))*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)), \
00517                            apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb))), \
00518                            (1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb)))));
00519                 Vs(1+2*q) += itpp::sqr(im_part(k))*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(im_bin_part.get_row(k)), \
00520                              apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb))), \
00521                              (1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb)))));
00522             }
00523             Vs(2*q) -= itpp::sqr(Es(2*q));
00524             Vs(1+2*q) -= itpp::sqr(Es(1+2*q));
00525         }
00526 
00527         //find equivalent channel matrix
00528         EquivCh(H_eq, c_impulse_response.get_col(ns));
00529         //compute invariant inverse
00530         K = H_eq*diag(Vs)*H_eq.transpose()+sigma2*itpp::eye(2*block_duration*nb_rec_ant);
00531         K_inv = itpp::inv(K);
00532         //find equivalent received signal
00533         EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
00534         for (q=0; q<symbols_block; q++)//symbols/block
00535         {
00536             //compute the extrinsic information of coded bits
00537             //real part
00538             //IC + filtering (real and imaginary parts of one symbol)
00539             interf = Es;
00540             interf(2*q) = 0;//this is the symbol to recover
00541             temp = H_eq.get_col(2*q);
00542             w = (part_var*temp.transpose())*(K_inv-itpp::outer_product(K_inv*(((part_var-Vs(2*q))/ \
00543                                              (1+((temp.transpose()*K_inv)*(temp*(part_var-Vs(2*q))))(0)))*temp), K_inv.transpose()*temp));
00544             s_tilde = w*(x_eq-H_eq*interf);
00545             mu_res = w*temp;//mean of the filtered signal
00546             index = select_half*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
00547             for (k=0; k<half_nb_bits_symb; k++)
00548             {
00549                 nom = -INFINITY;
00550                 denom = -INFINITY;
00551                 for (s=0; s<half_const_len; s++)
00552                 {
00553                     sigma2_res = ((w.transpose()*(K+itpp::outer_product(temp, temp*(itpp::sqr(re_part(s))-Vs(2*q)))))*w)(0)-itpp::sqr(re_part(s)*mu_res);//variance of the filtered signal
00554                     tmp = -itpp::sqr(s_tilde-mu_res*re_part(s))/(2*sigma2_res)+ \
00555                           itpp::to_vec(re_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
00556                     if (re_bin_part(s,k))
00557                         nom = std::max(nom, tmp);
00558                     else
00559                         denom = std::max(denom, tmp);
00560                 }
00561                 extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
00562             }
00563             //end real part
00564             //imaginary part
00565             //IC + filtering (real and imaginary parts of one symbol)
00566             interf = Es;
00567             interf(2*q+1) = 0;//this is the symbol to recover
00568             temp = H_eq.get_col(2*q+1);
00569             w = (part_var*temp.transpose())*(K_inv-itpp::outer_product(K_inv*(((part_var-Vs(1+2*q))/ \
00570                                              (1+((temp.transpose()*K_inv)*(temp*(part_var-Vs(1+2*q))))(0)))*temp), K_inv.transpose()*temp));
00571             s_tilde = w*(x_eq-H_eq*interf);
00572             mu_res = w*temp;//mean of the filtered signal
00573             index = (1-select_half)*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
00574             for (k=0; k<half_nb_bits_symb; k++)
00575             {
00576                 nom = -INFINITY;
00577                 denom = -INFINITY;
00578                 for (s=0; s<half_const_len; s++)
00579                 {
00580                     sigma2_res = ((w.transpose()*(K+itpp::outer_product(temp, temp*(itpp::sqr(im_part(s))-Vs(1+2*q)))))*w)(0)-itpp::sqr(im_part(s)*mu_res);//variance of the filtered signal
00581                     tmp = -itpp::sqr(s_tilde-mu_res*im_part(s))/(2*sigma2_res)+ \
00582                           itpp::to_vec(im_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
00583                     if (im_bin_part(s,k))
00584                         nom = std::max(nom, tmp);
00585                     else
00586                         denom = std::max(denom, tmp);
00587                 }
00588                 extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
00589             }
00590             //end imaginary part
00591         }//symbols/block
00592     }//block by block
00593 }
00594 
00595 void SISO::zfPIC(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
00596 //ZF Parallel Interference Canceller
00597 {
00598     //general parameters
00599     int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks
00600     int half_nb_bits_symb = nb_bits_symb/2;
00601     int half_const_len = itpp::pow2i(half_nb_bits_symb);
00602     int nb_bits_subblock = nb_bits_symb*symbols_block;//number of coded bits in an ST block
00603     itpp::vec Es(2*symbols_block);
00604     itpp::vec Vs(2*symbols_block);
00605     itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
00606     itpp::mat K(2*nb_rec_ant*block_duration,2*nb_rec_ant*block_duration);
00607     itpp::vec x_eq(2*nb_rec_ant*block_duration);
00608     itpp::vec interf(2*symbols_block);
00609     itpp::vec temp(2*nb_rec_ant*block_duration);
00610     itpp::vec w(2*nb_rec_ant*block_duration);//filter impulse response
00611     double s_tilde;
00612     double mu_res;
00613     double sigma2_res;
00614     double nom,denom;
00615     double tmp;
00616     register int ns,q,k,s;
00617     int index;
00618 
00619     //correspondence between real and imaginary part of symbols and their binary representations
00620     int select_half;
00621     itpp::vec re_part;
00622     itpp::bmat re_bin_part;
00623     itpp::vec im_part;
00624     itpp::bmat im_bin_part;
00625     find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
00626 
00627     extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
00628     for (ns=0; ns<nb_subblocks; ns++)//compute block by block
00629     {
00630         //mean and variance of real and imaginary parts of emitted symbols
00631         Es.zeros();
00632         Vs.zeros();
00633         for (q=0; q<symbols_block; q++)
00634         {
00635             index = q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
00636             for (k=0; k<half_const_len; k++)
00637             {
00638                 Es(2*q) += re_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)), \
00639                                                  apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb))), \
00640                                                  (1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb)))));
00641                 Es(1+2*q) += im_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(im_bin_part.get_row(k)), \
00642                                                    apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb))), \
00643                                                    (1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb)))));
00644                 Vs(2*q) += itpp::sqr(re_part(k))*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)), \
00645                            apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb))), \
00646                            (1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb)))));
00647                 Vs(1+2*q) += itpp::sqr(im_part(k))*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(to_vec(im_bin_part.get_row(k)), \
00648                              apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb))), \
00649                              (1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb)))));
00650             }
00651             Vs(2*q) -= itpp::sqr(Es(2*q));
00652             Vs(1+2*q) -= itpp::sqr(Es(1+2*q));
00653         }
00654 
00655         //find equivalent channel matrix
00656         EquivCh(H_eq, c_impulse_response.get_col(ns));
00657         //compute invariant inverse
00658         K = H_eq*itpp::diag(Vs)*H_eq.transpose()+sigma2*itpp::eye(2*block_duration*nb_rec_ant);
00659         //find equivalent received signal
00660         EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
00661         for (q=0; q<symbols_block; q++)//symbols/block
00662         {
00663             //compute the extrinsic information of coded bits
00664             //real part
00665             //IC + filtering (real and imaginary parts of one symbol)
00666             interf = Es;
00667             interf(2*q) = 0;//this is the symbol to recover
00668             temp = H_eq.get_col(2*q);
00669             w = temp/(temp*temp);//filter impulse response
00670             s_tilde = w*(x_eq-H_eq*interf);
00671             mu_res = w*temp;//mean of the filtered signal
00672             index = select_half*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
00673             for (k=0; k<half_nb_bits_symb; k++)
00674             {
00675                 nom = -INFINITY;
00676                 denom = -INFINITY;
00677                 for (s=0; s<half_const_len; s++)
00678                 {
00679                     sigma2_res = ((w.transpose()*(K+itpp::outer_product(temp, temp*(itpp::sqr(re_part(s))-Vs(2*q)))))*w)(0)-itpp::sqr(re_part(s)*mu_res);
00680                     tmp = -itpp::sqr(s_tilde-mu_res*re_part(s))/(2*sigma2_res)+ \
00681                           itpp::to_vec(re_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
00682                     if (re_bin_part(s,k))
00683                         nom = std::max(nom, tmp);
00684                     else
00685                         denom = std::max(denom, tmp);
00686                 }
00687                 extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
00688             }
00689             //end real part
00690             //imaginary part
00691             //IC + filtering (real and imaginary parts of one symbol)
00692             interf = Es;
00693             interf(2*q+1) = 0;//this is the symbol to recover
00694             temp = H_eq.get_col(2*q+1);
00695             w = temp/(temp*temp);//filter impulse response
00696             s_tilde = w*(x_eq-H_eq*interf);
00697             mu_res = w*temp;//mean of the filtered signal
00698             index = (1-select_half)*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
00699             for (k=0; k<half_nb_bits_symb; k++)
00700             {
00701                 nom = -INFINITY;
00702                 denom = -INFINITY;
00703                 for (s=0; s<half_const_len; s++)
00704                 {
00705                     sigma2_res = ((w.transpose()*(K+itpp::outer_product(temp, temp*(itpp::sqr(im_part(s))-Vs(1+2*q)))))*w)(0)-itpp::sqr(im_part(s)*mu_res);
00706                     tmp = -itpp::sqr(s_tilde-mu_res*im_part(s))/(2*sigma2_res)+ \
00707                           itpp::to_vec(im_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
00708                     if (im_bin_part(s,k))
00709                         nom = std::max(nom, tmp);
00710                     else
00711                         denom = std::max(denom, tmp);
00712                 }
00713                 extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
00714             }
00715             //end imaginary part
00716         }//symbols/block
00717     }//block by block
00718 }
00719 
00720 void SISO::Alamouti_maxlogMAP(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
00721 //maxlogMAP algorithm for Alamouti ST code
00722 {
00723     //matched filter
00724     int int_len = apriori_data.length();//interleaver length
00725     int nb_symb = (int)(int_len/nb_bits_symb);//number of symbols/block
00726     itpp::cvec comb_sig(nb_symb);
00727     comb_sig.zeros();
00728     itpp::cmat conj_H = itpp::conj(c_impulse_response);
00729     itpp::cmat conj_X = itpp::conj(rec_sig);
00730     register int nr,n,cs;
00731     for (nr=0; nr<nb_rec_ant; nr++)
00732     {
00733         for (n=0; n<(nb_symb/2); n++)
00734         {
00735             comb_sig(2*n) += (conj_H(2*nr,n)*rec_sig(2*n,nr)+c_impulse_response(1+2*nr,n)*conj_X(1+2*n,nr));
00736             comb_sig(1+2*n) += (conj_H(1+2*nr,n)*rec_sig(2*n,nr)-c_impulse_response(2*nr,n)*conj_X(1+2*n,nr));
00737         }
00738     }
00739 
00740     //extrinsic information of coded bits
00741     int const_size = itpp::pow2i(nb_bits_symb);//constellation size
00742     double buffer;
00743     double nom,denom;
00744     double temp;
00745     int index;
00746     extrinsic_data.set_size(nb_bits_symb*nb_symb);
00747     for (n=0; n<nb_symb; n++)
00748     {
00749         buffer = itpp::sum_sqr(itpp::abs(c_impulse_response.get_col(n/2)));
00750         for (nr=0; nr<nb_bits_symb; nr++)
00751         {
00752             nom = -INFINITY;
00753             denom = -INFINITY;
00754             for (cs=0; cs<const_size; cs++)
00755             {
00756                 temp = -itpp::sqr(comb_sig(n)-buffer*constellation(cs))/(2*buffer*sigma2)+ \
00757                        itpp::to_vec(bin_constellation.get_row(cs))*apriori_data.mid(n*nb_bits_symb, nb_bits_symb);
00758                 if (bin_constellation(cs,nr))
00759                     nom = std::max(nom, temp);
00760                 else
00761                     denom = std::max(denom, temp);
00762             }
00763             index = n*nb_bits_symb+nr;
00764             extrinsic_data(index) = (nom-denom)-apriori_data(index);//extrinsic information
00765         }
00766     }
00767 }
00768 
00769 void SISO::demodulator_logMAP(itpp::vec &extrinsic_data, const itpp::cvec &rec_sig, const itpp::vec &apriori_data)
00771 {
00772     int nb_symb = rec_sig.length();
00773     int const_size = itpp::pow2i(nb_bits_symb);
00774     double nom,denom,temp;
00775     register int k,i,cs;
00776     int index;
00777     extrinsic_data.set_size(nb_bits_symb*nb_symb);
00778     for (k=0; k<nb_symb; k++)
00779     {
00780         for (i=0; i<nb_bits_symb; i++)
00781         {
00782             nom = 0;
00783             denom = 0;
00784             for (cs=0; cs<const_size; cs++)
00785             {
00786                 temp = -itpp::sqr(rec_sig(k)-c_impulse_response(0,k)*constellation(cs))/(2*sigma2)+\
00787                        itpp::to_vec(bin_constellation.get_row(cs))*apriori_data.mid(k*nb_bits_symb, nb_bits_symb);
00788                 if (bin_constellation(cs,i))
00789                     nom += std::exp(temp);
00790                 else
00791                     denom += std::exp(temp);
00792             }
00793             index = k*nb_bits_symb+i;
00794             extrinsic_data(index) = std::log(nom/denom)-apriori_data(index);//extrinsic information
00795         }
00796     }
00797 }
00798 
00799 void SISO::demodulator_maxlogMAP(itpp::vec &extrinsic_data, const itpp::cvec &rec_sig, const itpp::vec &apriori_data)
00801 {
00802     int nb_symb = rec_sig.length();
00803     int const_size = itpp::pow2i(nb_bits_symb);
00804     double nom,denom,temp;
00805     register int k,i,cs;
00806     int index;
00807     extrinsic_data.set_size(nb_bits_symb*nb_symb);
00808     for (k=0; k<nb_symb; k++)
00809     {
00810         for (i=0; i<nb_bits_symb; i++)
00811         {
00812             nom = -INFINITY;
00813             denom = -INFINITY;
00814             for (cs=0; cs<const_size; cs++)
00815             {
00816                 temp = -itpp::sqr(rec_sig(k)-c_impulse_response(0,k)*constellation(cs))/(2*sigma2)+\
00817                        itpp::to_vec(bin_constellation.get_row(cs))*apriori_data.mid(k*nb_bits_symb, nb_bits_symb);
00818                 if (bin_constellation(cs,i))
00819                     nom = std::max(nom, temp);
00820                 else
00821                     denom = std::max(denom, temp);
00822             }
00823             index = k*nb_bits_symb+i;
00824             extrinsic_data(index) = (nom-denom)-apriori_data(index);//extrinsic information
00825         }
00826     }
00827 }
00828 }//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