IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
siso_dem.cpp
Go to the documentation of this file.
1 
29 #include <itpp/comm/siso.h>
30 #include <limits>
31 #ifndef INFINITY
32 #define INFINITY std::numeric_limits<double>::infinity()
33 #endif
34 
35 namespace itpp
36 {
37 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)
38 /* finds real in imaginary parts of the constellation and its corresponding bits
39  * this approach is used for equivalent channel according to Hassibi's model
40  * the constellation must be quadratic and the number of bits per symbol must be a multiple of two
41  */
42 {
43  //values needed for initializations
44  int const_size = itpp::pow2i(nb_bits_symb);//constellation size
45  int half_nb_bits_symb = nb_bits_symb/2;
46  int half_len = itpp::pow2i(half_nb_bits_symb);//number of values of real(imaginary) part
47  //initialize output variables
48  select_half = 0;
49  re_part.set_size(half_len);
50  re_bin_part.set_size(half_len, half_nb_bits_symb);
51  re_part.zeros();
52  re_part(0) = constellation(0).real();
53  im_part.set_size(half_len);
54  im_bin_part.set_size(half_len, half_nb_bits_symb);
55  im_part.zeros();
56  im_part(0) = constellation(0).imag();
57  //select half for real (imaginary) to binary correspondence
58  if (nb_bits_symb%2)
59  {
60  print_err_msg("SISO::find_half_const: number of bits per symbol must be a multiple of two");
61  return;
62  }
63  const double min_diff = 1e-3;
64  itpp::ivec idx = itpp::find(itpp::abs(itpp::real(constellation)-re_part(0))<min_diff);
65  if (idx.length()!=half_len)
66  {
67  print_err_msg("SISO::find_half_const: the constellation must be quadratic");
68  return;
69  }
70  itpp::bvec temp(nb_bits_symb);
71  register int n;
72  for (n=0; n<2; n++)
73  {
74  temp = bin_constellation.get_row(idx(n));
75  re_bin_part.set_row(n,temp(0,half_nb_bits_symb-1));
76  }
77  select_half = (re_bin_part.get_row(0)==re_bin_part.get_row(1))?0:1;
78  //algorithm
79  double buffer;
80  temp = bin_constellation.get_row(0);
81  re_bin_part.set_row(0,temp(select_half*half_nb_bits_symb,(1+select_half)*half_nb_bits_symb-1));
82  im_bin_part.set_row(0,temp((1-select_half)*half_nb_bits_symb,(2-select_half)*half_nb_bits_symb-1));
83  int re_idx = 0;
84  int im_idx = 0;
85  for (n=1; n<const_size; n++)
86  {
87  temp = bin_constellation.get_row(n);
88  buffer = constellation(n).real();
89  if (itpp::prod(re_part-buffer)>min_diff)
90  {
91  re_idx++;
92  re_part(re_idx) = buffer;
93  re_bin_part.set_row(re_idx, temp(select_half*half_nb_bits_symb,(1+select_half)*half_nb_bits_symb-1));
94  }
95  buffer = constellation(n).imag();
96  if (itpp::prod(im_part-buffer)>min_diff)
97  {
98  im_idx++;
99  im_part(im_idx) = buffer;
100  im_bin_part.set_row(im_idx, temp((1-select_half)*half_nb_bits_symb,(2-select_half)*half_nb_bits_symb-1));
101  }
102  }
103 }
104 
105 void SISO::EquivRecSig(itpp::vec &x_eq, const itpp::cmat &rec_sig)
106 //finds equivalent received signal with real coefficients
107 //the equivalent received signal follows the model of Hassibi's paper
108 //ouput:
109 //x_eq - equivalent received signal with real coefficients
110 //inputs:
111 //rec_sig - received signal
112 {
113  for (int k=0; k<nb_rec_ant; k++)
114  {
115  x_eq.set_subvector(k*2*block_duration, itpp::real(rec_sig.get_col(k)));
116  x_eq.set_subvector(k*2*block_duration+block_duration, itpp::imag(rec_sig.get_col(k)));
117  }
118 }
119 
120 void SISO::EquivCh(itpp::mat &H_eq, const itpp::cvec &H)
121 //finds equivalent channel with real coefficients following the model of Hassibi's paper
122 //output:
123 //H_eq - equivalent channel
124 //input:
125 //H - channel matrix
126 {
127  itpp::mat Aq(2*block_duration,2*nb_em_ant);
128  itpp::mat Bq(2*block_duration,2*nb_em_ant);
129  itpp::cmat temp(block_duration,nb_em_ant);
130  itpp::vec h(2*nb_em_ant);
131  itpp::mat AhBh(2*block_duration,2);
132  register int n,k;
133  for (k=0; k<symbols_block; k++)
134  {
135  temp = ST_gen1.get(k*block_duration,k*block_duration+block_duration-1,0,nb_em_ant-1);
136  Aq.set_submatrix(0, 0, itpp::real(temp));
137  Aq.set_submatrix(0, nb_em_ant, -itpp::imag(temp));
138  Aq.set_submatrix(block_duration, 0, itpp::imag(temp));
139  Aq.set_submatrix(block_duration, nb_em_ant, itpp::real(temp));
140  temp = ST_gen2.get(k*block_duration,k*block_duration+block_duration-1,0,nb_em_ant-1);
141  Bq.set_submatrix(0, 0, -itpp::imag(temp));
142  Bq.set_submatrix(0, nb_em_ant, -itpp::real(temp));
143  Bq.set_submatrix(block_duration, 0, itpp::real(temp));
144  Bq.set_submatrix(block_duration, nb_em_ant, -itpp::imag(temp));
145  for (n=0; n<nb_rec_ant; n++)
146  {
147  h.set_subvector(0, real(H.mid(n*nb_em_ant,nb_em_ant)));
148  h.set_subvector(nb_em_ant, imag(H.mid(n*nb_em_ant,nb_em_ant)));
149  AhBh.set_col(0, Aq*h);
150  AhBh.set_col(1, Bq*h);
151  H_eq.set_submatrix(2*block_duration*n, 2*k, AhBh);
152  }
153  }
154 }
155 
156 void SISO::Hassibi_maxlogMAP(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
157 //maxlogMAP algorithm for ST block codes using Hassibi's model
158 {
159  //general parameters
160  int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks of ST matrices/int period
161  double N0 = 2*sigma2;//noise DSP
162  int nb_all_symb = itpp::pow2i(nb_bits_symb*symbols_block);//nb. of all possible input symbols as a binary vector
163  double nom,denom;//nominator and denominator of extrinsic information
164  itpp::bvec bin_frame(nb_bits_symb*symbols_block);//binary frame at channel input
165  itpp::bmat mat_bin_frame(nb_bits_symb, symbols_block);
166  itpp::vec symb_frame_eq(2*symbols_block);//frame of symbols at equivalent channel input
167  double temp;
168  itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);//equivalent channel matrix
169  itpp::vec x_eq(2*block_duration*nb_rec_ant);//equivalent received signal
170  register int ns,q,nb,n,k;
171  int index;
172  extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
173  //main loop
174  for (ns=0; ns<nb_subblocks; ns++)//for each subblock
175  {
176  //find equivalent channel matrix
177  EquivCh(H_eq, c_impulse_response.get_col(ns));
178  //find equivalent received signal
179  EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
180  //compute the LLR of each bit in a frame of symbols_block symbols
181  for (q=0; q<symbols_block; q++)//for each symbol in a subblock
182  {
183  for (nb=0; nb<nb_bits_symb; nb++)//for a given bit try all possible sollutions for the input symbol vector
184  {
185  nom = -INFINITY;
186  denom = -INFINITY;
187  for (n=0; n<nb_all_symb; n++)//all possible symbols
188  {
189  bin_frame = itpp::dec2bin(nb_bits_symb*symbols_block, n);
190  mat_bin_frame = itpp::reshape(bin_frame, nb_bits_symb, symbols_block);
191  for (k=0; k<symbols_block; k++)
192  {
193  symb_frame_eq(2*k) = constellation(itpp::bin2dec(mat_bin_frame.get_col(k))).real();
194  symb_frame_eq(1+2*k) = constellation(itpp::bin2dec(mat_bin_frame.get_col(k))).imag();
195  }
196  temp = -itpp::sum_sqr(x_eq-H_eq*symb_frame_eq)/N0+\
197  itpp::to_vec(bin_frame)*apriori_data.mid(ns*nb_bits_symb*symbols_block,nb_bits_symb*symbols_block);
198  if (bin_frame(nb+q*nb_bits_symb))
199  nom = std::max(nom, temp);
200  else
201  denom = std::max(denom, temp);
202  }
203  index = nb+q*nb_bits_symb+ns*nb_bits_symb*symbols_block;
204  extrinsic_data(index) = (nom-denom)-apriori_data(index);
205  }//bits/symbol
206  }//symbols/subblock
207  }//subblocks
208 }
209 
210 void SISO::GA(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
211 // Gaussian Approximation algorithm for ST codes using Hassibi's model
212 {
213  //general parameters
214  int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks
215  int half_nb_bits_symb = nb_bits_symb/2;
216  int half_len = itpp::pow2i(half_nb_bits_symb);//number of values of real(imaginary) part
217 
218  //correspondence between real and imaginary part of symbols and their binary representations
219  int select_half;
220  itpp::vec re_part;
221  itpp::bmat re_bin_part;
222  itpp::vec im_part;
223  itpp::bmat im_bin_part;
224  find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
225 
226  //equivalent channel
227  itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
228  itpp::vec E_re_s(symbols_block);
229  itpp::vec E_im_s(symbols_block);
230  itpp::vec Var_re_s(symbols_block);
231  itpp::vec Var_im_s(symbols_block);
232  itpp::vec Ey(2*block_duration*nb_rec_ant);
233  itpp::mat Cy(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
234  itpp::mat Cy_inv(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
235  itpp::vec x_eq(2*block_duration*nb_rec_ant);
236  itpp::vec EZeta(2*block_duration*nb_rec_ant);
237  itpp::mat CZeta_inv(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
238  double nom,denom;
239  double temp;
240  register int ns,q,k,p,cs;
241  int index;
242  extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
243  for (ns=0; ns<nb_subblocks; ns++)//subblock by subblock
244  {
245  //mean and variance of real and imaginary parts of emitted symbols
246  E_re_s.zeros();
247  E_im_s.zeros();
248  Var_re_s.zeros();
249  Var_im_s.zeros();
250  for (q=0; q<symbols_block; q++)
251  {
252  index = q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
253  for (k=0; k<half_len; k++)
254  {
255  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)),\
256  apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))),\
257  1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))));
258  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)),\
259  apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))),\
260  1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))));
261  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)),\
262  apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))),\
263  1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))));
264  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)),\
265  apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))),\
266  1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))));
267  }
268  Var_re_s(q) -= itpp::sqr(E_re_s(q));
269  Var_im_s(q) -= itpp::sqr(E_im_s(q));
270  }
271 
272  //find equivalent channel
273  EquivCh(H_eq, c_impulse_response.get_col(ns));
274 
275  //compute E[y] and Cov[y]
276  Ey.zeros();
277  Cy = sigma2*itpp::eye(2*block_duration*nb_rec_ant);
278  for (q=0; q<symbols_block; q++)
279  {
280  //real & imaginary
281  Ey += (H_eq.get_col(2*q)*E_re_s(q)+H_eq.get_col(1+2*q)*E_im_s(q));
282  Cy += (itpp::outer_product(H_eq.get_col(2*q), H_eq.get_col(2*q)*Var_re_s(q))+\
283  itpp::outer_product(H_eq.get_col(1+2*q), H_eq.get_col(1+2*q)*Var_im_s(q)));
284  }
285 
286  //inverse of Cov[y]
287  Cy_inv = itpp::inv(Cy);
288 
289  //find equivalent received signal
290  EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
291 
292  //compute extrinsic information of coded bits
293  for (q=0; q<symbols_block; q++)
294  {
295  //real part
296  EZeta = Ey-H_eq.get_col(2*q)*E_re_s(q);
297  CZeta_inv = Cy_inv+itpp::outer_product(Cy_inv*\
298  ((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)))*\
299  H_eq.get_col(2*q)), Cy_inv.transpose()*H_eq.get_col(2*q));
300  index = select_half*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
301  for (p=0; p<half_nb_bits_symb; p++)
302  {
303  nom = -INFINITY;
304  denom = -INFINITY;
305  for (cs=0; cs<half_len; cs++)
306  {
307  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)+\
308  itpp::to_vec(re_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
309  if (re_bin_part(cs,p))
310  nom = std::max(nom, temp);
311  else
312  denom = std::max(denom, temp);
313  }
314  extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
315  }
316  //imaginary part
317  EZeta = Ey-H_eq.get_col(1+2*q)*E_im_s(q);
318  CZeta_inv = Cy_inv+itpp::outer_product(Cy_inv*\
319  ((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)))*\
320  H_eq.get_col(1+2*q)), Cy_inv.transpose()*H_eq.get_col(1+2*q));
321  index = (1-select_half)*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
322  for (p=0; p<half_nb_bits_symb; p++)
323  {
324  nom = -INFINITY;
325  denom = -INFINITY;
326  for (cs=0; cs<half_len; cs++)
327  {
328  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)+\
329  itpp::to_vec(im_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
330  if (im_bin_part(cs,p))
331  nom = std::max(nom, temp);
332  else
333  denom = std::max(denom, temp);
334  }
335  extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
336  }
337  }
338  }//subblock by subblock
339 }
340 
341 void SISO::sGA(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
342 //simplified Gaussian Approximation algorithm for ST codes using Hassibi's model
343 {
344  //general parameters
345  int nb_subblocks = (int)(rec_sig.rows()/block_duration);//number of subblocks
346  int half_nb_bits_symb = (int)(nb_bits_symb/2);
347  int half_len = itpp::pow2i(half_nb_bits_symb);//number of values of real(imaginary) part
348 
349  //correspondence between real and imaginary part of symbols and their binary representations
350  int select_half;
351  itpp::vec re_part;
352  itpp::bmat re_bin_part;
353  itpp::vec im_part;
354  itpp::bmat im_bin_part;
355  find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
356 
357  //equivalent channel
358  itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
359 
360  itpp::vec E_re_s(symbols_block);
361  itpp::vec E_im_s(symbols_block);
362  itpp::vec Var_re_s(symbols_block);
363  itpp::vec Var_im_s(symbols_block);
364  itpp::vec Ey(2*block_duration*nb_rec_ant);
365  itpp::mat Cy(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
366  itpp::vec x_eq(2*block_duration*nb_rec_ant);
367  itpp::vec EZeta(2*block_duration*nb_rec_ant);
368  itpp::vec CZeta(2*block_duration*nb_rec_ant);
369  double nom,denom;
370  double temp;
371  register int ns,q,k,p,cs;
372  int index;
373  extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
374  for (ns=0; ns<nb_subblocks; ns++)//subblock by subblock
375  {
376  //mean and variance of real and imaginary parts of emitted symbols
377  E_re_s.zeros();
378  E_im_s.zeros();
379  Var_re_s.zeros();
380  Var_im_s.zeros();
381  for (q=0; q<symbols_block; q++)
382  {
383  index = q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
384  for (k=0; k<half_len; k++)
385  {
386  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)),\
387  apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))),\
388  1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))));
389  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)),\
390  apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))),\
391  1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))));
392  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)),\
393  apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))),\
394  1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))));
395  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)),\
396  apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))),\
397  1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))));
398  }
399  Var_re_s(q) -= itpp::sqr(E_re_s(q));
400  Var_im_s(q) -= itpp::sqr(E_im_s(q));
401  }
402 
403  //find equivalent channel
404  EquivCh(H_eq, c_impulse_response.get_col(ns));
405 
406  //compute E[y] and Cov[y]
407  Ey.zeros();
408  Cy = sigma2*itpp::eye(2*block_duration*nb_rec_ant);
409  for (q=0; q<symbols_block; q++)
410  {
411  //real & imaginary
412  Ey += (H_eq.get_col(2*q)*E_re_s(q)+H_eq.get_col(1+2*q)*E_im_s(q));
413  Cy += (itpp::outer_product(H_eq.get_col(2*q), H_eq.get_col(2*q)*Var_re_s(q))+\
414  itpp::outer_product(H_eq.get_col(1+2*q), H_eq.get_col(1+2*q)*Var_im_s(q)));
415  }
416 
417  //find equivalent received signal
418  EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
419 
420  //compute extrinsic INFINITYormation of coded bits
421  for (q=0; q<symbols_block; q++)
422  {
423  //real part
424  EZeta = Ey-H_eq.get_col(2*q)*E_re_s(q);
425  CZeta = diag(Cy-itpp::outer_product(H_eq.get_col(2*q), H_eq.get_col(2*q)*Var_re_s(q)));
426  index = select_half*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
427  for (p=0; p<half_nb_bits_symb; p++)
428  {
429  nom = -INFINITY;
430  denom = -INFINITY;
431  for (cs=0; cs<half_len; cs++)
432  {
433  temp = -0.5*itpp::sum(itpp::elem_div(sqr(x_eq-H_eq.get_col(2*q)*re_part(cs)-EZeta), CZeta))+\
434  itpp::to_vec(re_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
435  if (re_bin_part(cs,p))
436  nom = std::max(nom, temp);
437  else
438  denom = std::max(denom, temp);
439  }
440  extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
441  }
442  //imaginary part
443  EZeta = Ey-H_eq.get_col(1+2*q)*E_im_s(q);
444  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)));
445  index = (1-select_half)*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
446  for (p=0; p<half_nb_bits_symb; p++)
447  {
448  nom = -INFINITY;
449  denom = -INFINITY;
450  for (cs=0; cs<half_len; cs++)
451  {
452  temp = -0.5*itpp::sum(itpp::elem_div(sqr(x_eq-H_eq.get_col(1+2*q)*im_part(cs)-EZeta), CZeta))+\
453  itpp::to_vec(im_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
454  if (im_bin_part(cs,p))
455  nom = std::max(nom, temp);
456  else
457  denom = std::max(denom, temp);
458  }
459  extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
460  }
461  }
462  }//subblock by subblock
463 }
464 
465 void SISO::mmsePIC(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
466 //MMSE Parallel Interference Canceller
467 {
468  //general parameters
469  int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks
470  int half_nb_bits_symb = nb_bits_symb/2;
471  int half_const_len = itpp::pow2i(half_nb_bits_symb);
472  int nb_bits_subblock = nb_bits_symb*symbols_block;//number of coded bits in an ST block
473  itpp::vec Es(2*symbols_block);
474  itpp::vec Vs(2*symbols_block);
475  itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
476  itpp::mat K(2*nb_rec_ant*block_duration,2*nb_rec_ant*block_duration);
477  itpp::mat K_inv(2*nb_rec_ant*block_duration,2*nb_rec_ant*block_duration);
478  itpp::vec x_eq(2*nb_rec_ant*block_duration);
479  itpp::vec interf(2*symbols_block);
480  itpp::vec temp(2*nb_rec_ant*block_duration);
481  itpp::vec w(2*nb_rec_ant*block_duration);//filter impulse response
482  double s_tilde;
483  double mu_res;
484  double sigma2_res;
485  double nom,denom;
486  double tmp;
487  register int ns,q,k,s;
488  int index;
489 
490  //correspondence between real and imaginary part of symbols and their binary representations
491  int select_half;
492  itpp::vec re_part;
493  itpp::bmat re_bin_part;
494  itpp::vec im_part;
495  itpp::bmat im_bin_part;
496  find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
497  double part_var = 1/(double)(2*nb_em_ant);//real and imaginary part variance
498 
499  extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
500  for (ns=0; ns<nb_subblocks; ns++)//compute block by block
501  {
502  //mean and variance of real and imaginary parts of emitted symbols
503  Es.zeros();
504  Vs.zeros();
505  for (q=0; q<symbols_block; q++)
506  {
507  index = q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
508  for (k=0; k<half_const_len; k++)
509  {
510  Es(2*q) += re_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)), \
511  apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb))), \
512  (1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb)))));
513  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)), \
514  apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb))), \
515  (1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb)))));
516  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)), \
517  apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb))), \
518  (1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb)))));
519  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)), \
520  apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb))), \
521  (1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb)))));
522  }
523  Vs(2*q) -= itpp::sqr(Es(2*q));
524  Vs(1+2*q) -= itpp::sqr(Es(1+2*q));
525  }
526 
527  //find equivalent channel matrix
528  EquivCh(H_eq, c_impulse_response.get_col(ns));
529  //compute invariant inverse
530  K = H_eq*diag(Vs)*H_eq.transpose()+sigma2*itpp::eye(2*block_duration*nb_rec_ant);
531  K_inv = itpp::inv(K);
532  //find equivalent received signal
533  EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
534  for (q=0; q<symbols_block; q++)//symbols/block
535  {
536  //compute the extrinsic information of coded bits
537  //real part
538  //IC + filtering (real and imaginary parts of one symbol)
539  interf = Es;
540  interf(2*q) = 0;//this is the symbol to recover
541  temp = H_eq.get_col(2*q);
542  w = (part_var*temp.transpose())*(K_inv-itpp::outer_product(K_inv*(((part_var-Vs(2*q))/ \
543  (1+((temp.transpose()*K_inv)*(temp*(part_var-Vs(2*q))))(0)))*temp), K_inv.transpose()*temp));
544  s_tilde = w*(x_eq-H_eq*interf);
545  mu_res = w*temp;//mean of the filtered signal
546  index = select_half*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
547  for (k=0; k<half_nb_bits_symb; k++)
548  {
549  nom = -INFINITY;
550  denom = -INFINITY;
551  for (s=0; s<half_const_len; s++)
552  {
553  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
554  tmp = -itpp::sqr(s_tilde-mu_res*re_part(s))/(2*sigma2_res)+ \
555  itpp::to_vec(re_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
556  if (re_bin_part(s,k))
557  nom = std::max(nom, tmp);
558  else
559  denom = std::max(denom, tmp);
560  }
561  extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
562  }
563  //end real part
564  //imaginary part
565  //IC + filtering (real and imaginary parts of one symbol)
566  interf = Es;
567  interf(2*q+1) = 0;//this is the symbol to recover
568  temp = H_eq.get_col(2*q+1);
569  w = (part_var*temp.transpose())*(K_inv-itpp::outer_product(K_inv*(((part_var-Vs(1+2*q))/ \
570  (1+((temp.transpose()*K_inv)*(temp*(part_var-Vs(1+2*q))))(0)))*temp), K_inv.transpose()*temp));
571  s_tilde = w*(x_eq-H_eq*interf);
572  mu_res = w*temp;//mean of the filtered signal
573  index = (1-select_half)*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
574  for (k=0; k<half_nb_bits_symb; k++)
575  {
576  nom = -INFINITY;
577  denom = -INFINITY;
578  for (s=0; s<half_const_len; s++)
579  {
580  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
581  tmp = -itpp::sqr(s_tilde-mu_res*im_part(s))/(2*sigma2_res)+ \
582  itpp::to_vec(im_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
583  if (im_bin_part(s,k))
584  nom = std::max(nom, tmp);
585  else
586  denom = std::max(denom, tmp);
587  }
588  extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
589  }
590  //end imaginary part
591  }//symbols/block
592  }//block by block
593 }
594 
595 void SISO::zfPIC(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
596 //ZF Parallel Interference Canceller
597 {
598  //general parameters
599  int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks
600  int half_nb_bits_symb = nb_bits_symb/2;
601  int half_const_len = itpp::pow2i(half_nb_bits_symb);
602  int nb_bits_subblock = nb_bits_symb*symbols_block;//number of coded bits in an ST block
603  itpp::vec Es(2*symbols_block);
604  itpp::vec Vs(2*symbols_block);
605  itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
606  itpp::mat K(2*nb_rec_ant*block_duration,2*nb_rec_ant*block_duration);
607  itpp::vec x_eq(2*nb_rec_ant*block_duration);
608  itpp::vec interf(2*symbols_block);
609  itpp::vec temp(2*nb_rec_ant*block_duration);
610  itpp::vec w(2*nb_rec_ant*block_duration);//filter impulse response
611  double s_tilde;
612  double mu_res;
613  double sigma2_res;
614  double nom,denom;
615  double tmp;
616  register int ns,q,k,s;
617  int index;
618 
619  //correspondence between real and imaginary part of symbols and their binary representations
620  int select_half;
621  itpp::vec re_part;
622  itpp::bmat re_bin_part;
623  itpp::vec im_part;
624  itpp::bmat im_bin_part;
625  find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
626 
627  extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
628  for (ns=0; ns<nb_subblocks; ns++)//compute block by block
629  {
630  //mean and variance of real and imaginary parts of emitted symbols
631  Es.zeros();
632  Vs.zeros();
633  for (q=0; q<symbols_block; q++)
634  {
635  index = q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
636  for (k=0; k<half_const_len; k++)
637  {
638  Es(2*q) += re_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)), \
639  apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb))), \
640  (1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb)))));
641  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)), \
642  apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb))), \
643  (1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb)))));
644  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)), \
645  apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb))), \
646  (1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index, half_nb_bits_symb)))));
647  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)), \
648  apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb))), \
649  (1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index, half_nb_bits_symb)))));
650  }
651  Vs(2*q) -= itpp::sqr(Es(2*q));
652  Vs(1+2*q) -= itpp::sqr(Es(1+2*q));
653  }
654 
655  //find equivalent channel matrix
656  EquivCh(H_eq, c_impulse_response.get_col(ns));
657  //compute invariant inverse
658  K = H_eq*itpp::diag(Vs)*H_eq.transpose()+sigma2*itpp::eye(2*block_duration*nb_rec_ant);
659  //find equivalent received signal
660  EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
661  for (q=0; q<symbols_block; q++)//symbols/block
662  {
663  //compute the extrinsic information of coded bits
664  //real part
665  //IC + filtering (real and imaginary parts of one symbol)
666  interf = Es;
667  interf(2*q) = 0;//this is the symbol to recover
668  temp = H_eq.get_col(2*q);
669  w = temp/(temp*temp);//filter impulse response
670  s_tilde = w*(x_eq-H_eq*interf);
671  mu_res = w*temp;//mean of the filtered signal
672  index = select_half*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
673  for (k=0; k<half_nb_bits_symb; k++)
674  {
675  nom = -INFINITY;
676  denom = -INFINITY;
677  for (s=0; s<half_const_len; s++)
678  {
679  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);
680  tmp = -itpp::sqr(s_tilde-mu_res*re_part(s))/(2*sigma2_res)+ \
681  itpp::to_vec(re_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
682  if (re_bin_part(s,k))
683  nom = std::max(nom, tmp);
684  else
685  denom = std::max(denom, tmp);
686  }
687  extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
688  }
689  //end real part
690  //imaginary part
691  //IC + filtering (real and imaginary parts of one symbol)
692  interf = Es;
693  interf(2*q+1) = 0;//this is the symbol to recover
694  temp = H_eq.get_col(2*q+1);
695  w = temp/(temp*temp);//filter impulse response
696  s_tilde = w*(x_eq-H_eq*interf);
697  mu_res = w*temp;//mean of the filtered signal
698  index = (1-select_half)*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
699  for (k=0; k<half_nb_bits_symb; k++)
700  {
701  nom = -INFINITY;
702  denom = -INFINITY;
703  for (s=0; s<half_const_len; s++)
704  {
705  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);
706  tmp = -itpp::sqr(s_tilde-mu_res*im_part(s))/(2*sigma2_res)+ \
707  itpp::to_vec(im_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
708  if (im_bin_part(s,k))
709  nom = std::max(nom, tmp);
710  else
711  denom = std::max(denom, tmp);
712  }
713  extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
714  }
715  //end imaginary part
716  }//symbols/block
717  }//block by block
718 }
719 
720 void SISO::Alamouti_maxlogMAP(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
721 //maxlogMAP algorithm for Alamouti ST code
722 {
723  //matched filter
724  int int_len = apriori_data.length();//interleaver length
725  int nb_symb = (int)(int_len/nb_bits_symb);//number of symbols/block
726  itpp::cvec comb_sig(nb_symb);
727  comb_sig.zeros();
728  itpp::cmat conj_H = itpp::conj(c_impulse_response);
729  itpp::cmat conj_X = itpp::conj(rec_sig);
730  register int nr,n,cs;
731  for (nr=0; nr<nb_rec_ant; nr++)
732  {
733  for (n=0; n<(nb_symb/2); n++)
734  {
735  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));
736  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));
737  }
738  }
739 
740  //extrinsic information of coded bits
741  int const_size = itpp::pow2i(nb_bits_symb);//constellation size
742  double buffer;
743  double nom,denom;
744  double temp;
745  int index;
746  extrinsic_data.set_size(nb_bits_symb*nb_symb);
747  for (n=0; n<nb_symb; n++)
748  {
749  buffer = itpp::sum_sqr(itpp::abs(c_impulse_response.get_col(n/2)));
750  for (nr=0; nr<nb_bits_symb; nr++)
751  {
752  nom = -INFINITY;
753  denom = -INFINITY;
754  for (cs=0; cs<const_size; cs++)
755  {
756  temp = -itpp::sqr(comb_sig(n)-buffer*constellation(cs))/(2*buffer*sigma2)+ \
757  itpp::to_vec(bin_constellation.get_row(cs))*apriori_data.mid(n*nb_bits_symb, nb_bits_symb);
758  if (bin_constellation(cs,nr))
759  nom = std::max(nom, temp);
760  else
761  denom = std::max(denom, temp);
762  }
763  index = n*nb_bits_symb+nr;
764  extrinsic_data(index) = (nom-denom)-apriori_data(index);//extrinsic information
765  }
766  }
767 }
768 
769 void SISO::demodulator_logMAP(itpp::vec &extrinsic_data, const itpp::cvec &rec_sig, const itpp::vec &apriori_data)
771 {
772  int nb_symb = rec_sig.length();
773  int const_size = itpp::pow2i(nb_bits_symb);
774  double nom,denom,temp;
775  register int k,i,cs;
776  int index;
777  extrinsic_data.set_size(nb_bits_symb*nb_symb);
778  for (k=0; k<nb_symb; k++)
779  {
780  for (i=0; i<nb_bits_symb; i++)
781  {
782  nom = 0;
783  denom = 0;
784  for (cs=0; cs<const_size; cs++)
785  {
786  temp = -itpp::sqr(rec_sig(k)-c_impulse_response(0,k)*constellation(cs))/(2*sigma2)+\
787  itpp::to_vec(bin_constellation.get_row(cs))*apriori_data.mid(k*nb_bits_symb, nb_bits_symb);
788  if (bin_constellation(cs,i))
789  nom += std::exp(temp);
790  else
791  denom += std::exp(temp);
792  }
793  index = k*nb_bits_symb+i;
794  extrinsic_data(index) = std::log(nom/denom)-apriori_data(index);//extrinsic information
795  }
796  }
797 }
798 
799 void SISO::demodulator_maxlogMAP(itpp::vec &extrinsic_data, const itpp::cvec &rec_sig, const itpp::vec &apriori_data)
801 {
802  int nb_symb = rec_sig.length();
803  int const_size = itpp::pow2i(nb_bits_symb);
804  double nom,denom,temp;
805  register int k,i,cs;
806  int index;
807  extrinsic_data.set_size(nb_bits_symb*nb_symb);
808  for (k=0; k<nb_symb; k++)
809  {
810  for (i=0; i<nb_bits_symb; i++)
811  {
812  nom = -INFINITY;
813  denom = -INFINITY;
814  for (cs=0; cs<const_size; cs++)
815  {
816  temp = -itpp::sqr(rec_sig(k)-c_impulse_response(0,k)*constellation(cs))/(2*sigma2)+\
817  itpp::to_vec(bin_constellation.get_row(cs))*apriori_data.mid(k*nb_bits_symb, nb_bits_symb);
818  if (bin_constellation(cs,i))
819  nom = std::max(nom, temp);
820  else
821  denom = std::max(denom, temp);
822  }
823  index = k*nb_bits_symb+i;
824  extrinsic_data(index) = (nom-denom)-apriori_data(index);//extrinsic information
825  }
826  }
827 }
828 }//end namespace tr
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:13 for IT++ by Doxygen 1.8.1.2