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
Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4