siso_dem.cpp
Go to the documentation of this file.
#include <itpp/comm/siso.h>
#include <limits>
#ifndef INFINITY
#define INFINITY std::numeric_limits<double>::infinity()
#endif
namespace itpp
{
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)
/* finds real in imaginary parts of the constellation and its corresponding bits
* this approach is used for equivalent channel according to Hassibi's model
* the constellation must be quadratic and the number of bits per symbol must be a multiple of two
*/
{
//values needed for initializations
int const_size = itpp::pow2i(nb_bits_symb);//constellation size
int half_nb_bits_symb = nb_bits_symb/2;
int half_len = itpp::pow2i(half_nb_bits_symb);//number of values of real(imaginary) part
//initialize output variables
select_half = 0;
re_part.set_size(half_len);
re_bin_part.set_size(half_len, half_nb_bits_symb);
re_part.zeros();
re_part(0) = constellation(0).real();
im_part.set_size(half_len);
im_bin_part.set_size(half_len, half_nb_bits_symb);
im_part.zeros();
im_part(0) = constellation(0).imag();
//select half for real (imaginary) to binary correspondence
if (nb_bits_symb%2)
{
print_err_msg("SISO::find_half_const: number of bits per symbol must be a multiple of two");
return;
}
const double min_diff = 1e-3;
itpp::ivec idx = itpp::find(itpp::abs(itpp::real(constellation)-re_part(0))<min_diff);
if (idx.length()!=half_len)
{
print_err_msg("SISO::find_half_const: the constellation must be quadratic");
return;
}
itpp::bvec temp(nb_bits_symb);
register int n;
for (n=0; n<2; n++)
{
temp = bin_constellation.get_row(idx(n));
re_bin_part.set_row(n,temp(0,half_nb_bits_symb-1));
}
select_half = (re_bin_part.get_row(0)==re_bin_part.get_row(1))?0:1;
//algorithm
double buffer;
temp = bin_constellation.get_row(0);
re_bin_part.set_row(0,temp(select_half*half_nb_bits_symb,(1+select_half)*half_nb_bits_symb-1));
im_bin_part.set_row(0,temp((1-select_half)*half_nb_bits_symb,(2-select_half)*half_nb_bits_symb-1));
int re_idx = 0;
int im_idx = 0;
for (n=1; n<const_size; n++)
{
temp = bin_constellation.get_row(n);
buffer = constellation(n).real();
if (fabs(itpp::prod(re_part-buffer))>min_diff)
{
re_idx++;
re_part(re_idx) = buffer;
re_bin_part.set_row(re_idx, temp(select_half*half_nb_bits_symb,(1+select_half)*half_nb_bits_symb-1));
}
buffer = constellation(n).imag();
if (fabs(itpp::prod(im_part-buffer))>min_diff)
{
im_idx++;
im_part(im_idx) = buffer;
im_bin_part.set_row(im_idx, temp((1-select_half)*half_nb_bits_symb,(2-select_half)*half_nb_bits_symb-1));
}
}
}
void SISO::EquivRecSig(itpp::vec &x_eq, const itpp::cmat &rec_sig)
//finds equivalent received signal with real coefficients
//the equivalent received signal follows the model of Hassibi's paper
//ouput:
//x_eq - equivalent received signal with real coefficients
//inputs:
//rec_sig - received signal
{
for (int k=0; k<nb_rec_ant; k++)
{
x_eq.set_subvector(k*2*block_duration, itpp::real(rec_sig.get_col(k)));
x_eq.set_subvector(k*2*block_duration+block_duration, itpp::imag(rec_sig.get_col(k)));
}
}
void SISO::EquivCh(itpp::mat &H_eq, const itpp::cvec &H)
//finds equivalent channel with real coefficients following the model of Hassibi's paper
//output:
//H_eq - equivalent channel
//input:
//H - channel matrix
{
itpp::mat Aq(2*block_duration,2*nb_em_ant);
itpp::mat Bq(2*block_duration,2*nb_em_ant);
itpp::cmat temp(block_duration,nb_em_ant);
itpp::vec h(2*nb_em_ant);
itpp::mat AhBh(2*block_duration,2);
register int n,k;
for (k=0; k<symbols_block; k++)
{
temp = ST_gen1.get(k*block_duration,k*block_duration+block_duration-1,0,nb_em_ant-1);
Aq.set_submatrix(0, 0, itpp::real(temp));
Aq.set_submatrix(0, nb_em_ant, -itpp::imag(temp));
Aq.set_submatrix(block_duration, 0, itpp::imag(temp));
Aq.set_submatrix(block_duration, nb_em_ant, itpp::real(temp));
temp = ST_gen2.get(k*block_duration,k*block_duration+block_duration-1,0,nb_em_ant-1);
Bq.set_submatrix(0, 0, -itpp::imag(temp));
Bq.set_submatrix(0, nb_em_ant, -itpp::real(temp));
Bq.set_submatrix(block_duration, 0, itpp::real(temp));
Bq.set_submatrix(block_duration, nb_em_ant, -itpp::imag(temp));
for (n=0; n<nb_rec_ant; n++)
{
h.set_subvector(0, real(H.mid(n*nb_em_ant,nb_em_ant)));
h.set_subvector(nb_em_ant, imag(H.mid(n*nb_em_ant,nb_em_ant)));
AhBh.set_col(0, Aq*h);
AhBh.set_col(1, Bq*h);
H_eq.set_submatrix(2*block_duration*n, 2*k, AhBh);
}
}
}
void SISO::compute_symb_stats(itpp::vec &Es, itpp::vec &Vs,
int ns, int select_half, const itpp::vec &apriori_data,
const itpp::vec &re_part, const itpp::vec &im_part,
const itpp::bmat &re_bin_part, const itpp::bmat &im_bin_part)
{
int half_nb_bits_symb = nb_bits_symb/2;
int half_len = itpp::pow2i(half_nb_bits_symb);//number of values of real(imaginary) part
Es.zeros();
Vs.zeros();
for (int q = 0; q < symbols_block; q++)
{
int index = q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
for (int k = 0; k < half_len; k++)
{
Es(2*q) += re_part(k)*itpp::prod(itpp::elem_div(exp(itpp::elem_mult(itpp::to_vec(re_bin_part.get_row(k)),\
apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))),\
1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))));
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)),\
apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))),\
1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))));
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)),\
apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))),\
1+exp(apriori_data.mid(select_half*half_nb_bits_symb+index,half_nb_bits_symb))));
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)),\
apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))),\
1+exp(apriori_data.mid((1-select_half)*half_nb_bits_symb+index,half_nb_bits_symb))));
}
Vs(2*q) -= itpp::sqr(Es(2*q));
Vs(1+2*q) -= itpp::sqr(Es(1+2*q));
}
}
void SISO::Hassibi_maxlogMAP(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
//maxlogMAP algorithm for ST block codes using Hassibi's model
{
//general parameters
int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks of ST matrices/int period
double N0 = 2*sigma2;//noise DSP
int nb_all_symb = itpp::pow2i(nb_bits_symb*symbols_block);//nb. of all possible input symbols as a binary vector
double nom,denom;//nominator and denominator of extrinsic information
itpp::bvec bin_frame(nb_bits_symb*symbols_block);//binary frame at channel input
itpp::bmat mat_bin_frame(nb_bits_symb, symbols_block);
itpp::vec symb_frame_eq(2*symbols_block);//frame of symbols at equivalent channel input
double temp;
itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);//equivalent channel matrix
itpp::vec x_eq(2*block_duration*nb_rec_ant);//equivalent received signal
register int ns,q,nb,n,k;
int index;
extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
//main loop
for (ns=0; ns<nb_subblocks; ns++)//for each subblock
{
//find equivalent channel matrix
EquivCh(H_eq, c_impulse_response.get_col(ns));
//find equivalent received signal
EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
//compute the LLR of each bit in a frame of symbols_block symbols
for (q=0; q<symbols_block; q++)//for each symbol in a subblock
{
for (nb=0; nb<nb_bits_symb; nb++)//for a given bit try all possible sollutions for the input symbol vector
{
nom = -INFINITY;
denom = -INFINITY;
for (n=0; n<nb_all_symb; n++)//all possible symbols
{
bin_frame = itpp::dec2bin(nb_bits_symb*symbols_block, n);
mat_bin_frame = itpp::reshape(bin_frame, nb_bits_symb, symbols_block);
for (k=0; k<symbols_block; k++)
{
symb_frame_eq(2*k) = constellation(itpp::bin2dec(mat_bin_frame.get_col(k))).real();
symb_frame_eq(1+2*k) = constellation(itpp::bin2dec(mat_bin_frame.get_col(k))).imag();
}
temp = -itpp::sum_sqr(x_eq-H_eq*symb_frame_eq)/N0+\
itpp::to_vec(bin_frame)*apriori_data.mid(ns*nb_bits_symb*symbols_block,nb_bits_symb*symbols_block);
if (bin_frame(nb+q*nb_bits_symb))
nom = std::max(nom, temp);
else
denom = std::max(denom, temp);
}
index = nb+q*nb_bits_symb+ns*nb_bits_symb*symbols_block;
extrinsic_data(index) = (nom-denom)-apriori_data(index);
}//bits/symbol
}//symbols/subblock
}//subblocks
}
void SISO::GA(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
// Gaussian Approximation algorithm for ST codes using Hassibi's model
{
//general parameters
int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks
int half_nb_bits_symb = nb_bits_symb/2;
int half_len = itpp::pow2i(half_nb_bits_symb);//number of values of real(imaginary) part
//correspondence between real and imaginary part of symbols and their binary representations
int select_half;
itpp::vec re_part;
itpp::bmat re_bin_part;
itpp::vec im_part;
itpp::bmat im_bin_part;
find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
//equivalent channel
itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
itpp::vec Es(2*symbols_block);
itpp::vec Vs(2*symbols_block);
itpp::vec Ey(2*block_duration*nb_rec_ant);
itpp::mat Cy(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
itpp::mat Cy_inv(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
itpp::vec x_eq(2*block_duration*nb_rec_ant);
itpp::vec EZeta(2*block_duration*nb_rec_ant);
itpp::mat CZeta_inv(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
double nom,denom;
double temp;
register int ns,q,p,cs;
int index;
extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
for (ns=0; ns<nb_subblocks; ns++)//subblock by subblock
{
//mean and variance of real and imaginary parts of emitted symbols
compute_symb_stats(Es, Vs, ns, select_half, apriori_data,
re_part, im_part, re_bin_part, im_bin_part);
//find equivalent channel
EquivCh(H_eq, c_impulse_response.get_col(ns));
//compute E[y] and Cov[y]
Ey.zeros();
Cy = sigma2*itpp::eye(2*block_duration*nb_rec_ant);
for (q=0; q<symbols_block; q++)
{
//real & imaginary
Ey += (H_eq.get_col(2*q)*Es(2*q)+H_eq.get_col(1+2*q)*Es(1+2*q));
Cy += (itpp::outer_product(H_eq.get_col(2*q), H_eq.get_col(2*q)*Vs(2*q))+\
itpp::outer_product(H_eq.get_col(1+2*q), H_eq.get_col(1+2*q)*Vs(1+2*q)));
}
//inverse of Cov[y]
Cy_inv = itpp::inv(Cy);
//find equivalent received signal
EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
//compute extrinsic information of coded bits
for (q=0; q<symbols_block; q++)
{
//real part
EZeta = Ey-H_eq.get_col(2*q)*Es(2*q);
CZeta_inv = Cy_inv+itpp::outer_product(Cy_inv*\
((Vs(2*q)/(1-(((H_eq.get_col(2*q)).transpose()*Cy_inv)*(H_eq.get_col(2*q)*Vs(2*q)))(0)))*\
H_eq.get_col(2*q)), Cy_inv.transpose()*H_eq.get_col(2*q));
index = select_half*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
for (p=0; p<half_nb_bits_symb; p++)
{
nom = -INFINITY;
denom = -INFINITY;
for (cs=0; cs<half_len; cs++)
{
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)+\
itpp::to_vec(re_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
if (re_bin_part(cs,p))
nom = std::max(nom, temp);
else
denom = std::max(denom, temp);
}
extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
}
//imaginary part
EZeta = Ey-H_eq.get_col(1+2*q)*Es(1+2*q);
CZeta_inv = Cy_inv+itpp::outer_product(Cy_inv*\
((Vs(1+2*q)/(1-(((H_eq.get_col(1+2*q)).transpose()*Cy_inv)*(H_eq.get_col(1+2*q)*Vs(1+2*q)))(0)))*\
H_eq.get_col(1+2*q)), Cy_inv.transpose()*H_eq.get_col(1+2*q));
index = (1-select_half)*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
for (p=0; p<half_nb_bits_symb; p++)
{
nom = -INFINITY;
denom = -INFINITY;
for (cs=0; cs<half_len; cs++)
{
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)+\
itpp::to_vec(im_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
if (im_bin_part(cs,p))
nom = std::max(nom, temp);
else
denom = std::max(denom, temp);
}
extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
}
}
}//subblock by subblock
}
void SISO::sGA(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
//simplified Gaussian Approximation algorithm for ST codes using Hassibi's model
{
//general parameters
int nb_subblocks = (int)(rec_sig.rows()/block_duration);//number of subblocks
int half_nb_bits_symb = (int)(nb_bits_symb/2);
int half_len = itpp::pow2i(half_nb_bits_symb);//number of values of real(imaginary) part
//correspondence between real and imaginary part of symbols and their binary representations
int select_half;
itpp::vec re_part;
itpp::bmat re_bin_part;
itpp::vec im_part;
itpp::bmat im_bin_part;
find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
//equivalent channel
itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
itpp::vec Es(2*symbols_block);
itpp::vec Vs(2*symbols_block);
itpp::vec Ey(2*block_duration*nb_rec_ant);
itpp::mat Cy(2*block_duration*nb_rec_ant,2*block_duration*nb_rec_ant);
itpp::vec x_eq(2*block_duration*nb_rec_ant);
itpp::vec EZeta(2*block_duration*nb_rec_ant);
itpp::vec CZeta(2*block_duration*nb_rec_ant);
double nom,denom;
double temp;
register int ns,q,p,cs;
int index;
extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
for (ns=0; ns<nb_subblocks; ns++)//subblock by subblock
{
//mean and variance of real and imaginary parts of emitted symbols
compute_symb_stats(Es, Vs, ns, select_half, apriori_data,
re_part, im_part, re_bin_part, im_bin_part);
//find equivalent channel
EquivCh(H_eq, c_impulse_response.get_col(ns));
//compute E[y] and Cov[y]
Ey.zeros();
Cy = sigma2*itpp::eye(2*block_duration*nb_rec_ant);
for (q=0; q<symbols_block; q++)
{
//real & imaginary
Ey += (H_eq.get_col(2*q)*Es(2*q)+H_eq.get_col(1+2*q)*Es(1+2*q));
Cy += (itpp::outer_product(H_eq.get_col(2*q), H_eq.get_col(2*q)*Vs(2*q))+\
itpp::outer_product(H_eq.get_col(1+2*q), H_eq.get_col(1+2*q)*Vs(1+2*q)));
}
//find equivalent received signal
EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
//compute extrinsic INFINITYormation of coded bits
for (q=0; q<symbols_block; q++)
{
//real part
EZeta = Ey-H_eq.get_col(2*q)*Es(2*q);
CZeta = diag(Cy-itpp::outer_product(H_eq.get_col(2*q), H_eq.get_col(2*q)*Vs(2*q)));
index = select_half*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
for (p=0; p<half_nb_bits_symb; p++)
{
nom = -INFINITY;
denom = -INFINITY;
for (cs=0; cs<half_len; cs++)
{
temp = -0.5*itpp::sum(itpp::elem_div(sqr(x_eq-H_eq.get_col(2*q)*re_part(cs)-EZeta), CZeta))+\
itpp::to_vec(re_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
if (re_bin_part(cs,p))
nom = std::max(nom, temp);
else
denom = std::max(denom, temp);
}
extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
}
//imaginary part
EZeta = Ey-H_eq.get_col(1+2*q)*Es(1+2*q);
CZeta = itpp::diag(Cy-itpp::outer_product(H_eq.get_col(1+2*q), H_eq.get_col(1+2*q)*Vs(1+2*q)));
index = (1-select_half)*half_nb_bits_symb+q*nb_bits_symb+ns*symbols_block*nb_bits_symb;
for (p=0; p<half_nb_bits_symb; p++)
{
nom = -INFINITY;
denom = -INFINITY;
for (cs=0; cs<half_len; cs++)
{
temp = -0.5*itpp::sum(itpp::elem_div(sqr(x_eq-H_eq.get_col(1+2*q)*im_part(cs)-EZeta), CZeta))+\
itpp::to_vec(im_bin_part.get_row(cs))*apriori_data.mid(index,half_nb_bits_symb);
if (im_bin_part(cs,p))
nom = std::max(nom, temp);
else
denom = std::max(denom, temp);
}
extrinsic_data(index+p) = (nom-denom)-apriori_data(index+p);
}
}
}//subblock by subblock
}
void SISO::mmsePIC(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
//MMSE Parallel Interference Canceller
{
//general parameters
int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks
int half_nb_bits_symb = nb_bits_symb/2;
int half_const_len = itpp::pow2i(half_nb_bits_symb);
int nb_bits_subblock = nb_bits_symb*symbols_block;//number of coded bits in an ST block
itpp::vec Es(2*symbols_block);
itpp::vec Vs(2*symbols_block);
itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
itpp::mat K(2*nb_rec_ant*block_duration,2*nb_rec_ant*block_duration);
itpp::mat K_inv(2*nb_rec_ant*block_duration,2*nb_rec_ant*block_duration);
itpp::vec x_eq(2*nb_rec_ant*block_duration);
itpp::vec interf(2*symbols_block);
itpp::vec temp(2*nb_rec_ant*block_duration);
itpp::vec w(2*nb_rec_ant*block_duration);//filter impulse response
double s_tilde;
double mu_res;
double sigma2_res;
double nom,denom;
double tmp;
register int ns,q,k,s;
int index;
//correspondence between real and imaginary part of symbols and their binary representations
int select_half;
itpp::vec re_part;
itpp::bmat re_bin_part;
itpp::vec im_part;
itpp::bmat im_bin_part;
find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
double part_var = 1/(double)(2*nb_em_ant);//real and imaginary part variance
extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
for (ns=0; ns<nb_subblocks; ns++)//compute block by block
{
//mean and variance of real and imaginary parts of emitted symbols
compute_symb_stats(Es, Vs, ns, select_half, apriori_data,
re_part, im_part, re_bin_part, im_bin_part);
//find equivalent channel matrix
EquivCh(H_eq, c_impulse_response.get_col(ns));
//compute invariant inverse
K = H_eq*diag(Vs)*H_eq.transpose()+sigma2*itpp::eye(2*block_duration*nb_rec_ant);
K_inv = itpp::inv(K);
//find equivalent received signal
EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
for (q=0; q<symbols_block; q++)//symbols/block
{
//compute the extrinsic information of coded bits
//real part
//IC + filtering (real and imaginary parts of one symbol)
interf = Es;
interf(2*q) = 0;//this is the symbol to recover
temp = H_eq.get_col(2*q);
w = (part_var*temp.transpose())*(K_inv-itpp::outer_product(K_inv*(((part_var-Vs(2*q))/ \
(1+((temp.transpose()*K_inv)*(temp*(part_var-Vs(2*q))))(0)))*temp), K_inv.transpose()*temp));
s_tilde = w*(x_eq-H_eq*interf);
mu_res = w*temp;//mean of the filtered signal
index = select_half*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
for (k=0; k<half_nb_bits_symb; k++)
{
nom = -INFINITY;
denom = -INFINITY;
for (s=0; s<half_const_len; s++)
{
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
tmp = -itpp::sqr(s_tilde-mu_res*re_part(s))/(2*sigma2_res)+ \
itpp::to_vec(re_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
if (re_bin_part(s,k))
nom = std::max(nom, tmp);
else
denom = std::max(denom, tmp);
}
extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
}
//end real part
//imaginary part
//IC + filtering (real and imaginary parts of one symbol)
interf = Es;
interf(2*q+1) = 0;//this is the symbol to recover
temp = H_eq.get_col(2*q+1);
w = (part_var*temp.transpose())*(K_inv-itpp::outer_product(K_inv*(((part_var-Vs(1+2*q))/ \
(1+((temp.transpose()*K_inv)*(temp*(part_var-Vs(1+2*q))))(0)))*temp), K_inv.transpose()*temp));
s_tilde = w*(x_eq-H_eq*interf);
mu_res = w*temp;//mean of the filtered signal
index = (1-select_half)*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
for (k=0; k<half_nb_bits_symb; k++)
{
nom = -INFINITY;
denom = -INFINITY;
for (s=0; s<half_const_len; s++)
{
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
tmp = -itpp::sqr(s_tilde-mu_res*im_part(s))/(2*sigma2_res)+ \
itpp::to_vec(im_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
if (im_bin_part(s,k))
nom = std::max(nom, tmp);
else
denom = std::max(denom, tmp);
}
extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
}
//end imaginary part
}//symbols/block
}//block by block
}
void SISO::zfPIC(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
//ZF Parallel Interference Canceller
{
//general parameters
int nb_subblocks = rec_sig.rows()/block_duration;//number of subblocks
int half_nb_bits_symb = nb_bits_symb/2;
int half_const_len = itpp::pow2i(half_nb_bits_symb);
int nb_bits_subblock = nb_bits_symb*symbols_block;//number of coded bits in an ST block
itpp::vec Es(2*symbols_block);
itpp::vec Vs(2*symbols_block);
itpp::mat H_eq(2*nb_rec_ant*block_duration,2*symbols_block);
itpp::mat K(2*nb_rec_ant*block_duration,2*nb_rec_ant*block_duration);
itpp::vec x_eq(2*nb_rec_ant*block_duration);
itpp::vec interf(2*symbols_block);
itpp::vec temp(2*nb_rec_ant*block_duration);
itpp::vec w(2*nb_rec_ant*block_duration);//filter impulse response
double s_tilde;
double mu_res;
double sigma2_res;
double nom,denom;
double tmp;
register int ns,q,k,s;
int index;
//correspondence between real and imaginary part of symbols and their binary representations
int select_half;
itpp::vec re_part;
itpp::bmat re_bin_part;
itpp::vec im_part;
itpp::bmat im_bin_part;
find_half_const(select_half, re_part, re_bin_part, im_part, im_bin_part);
extrinsic_data.set_size(nb_bits_symb*nb_subblocks*symbols_block);
for (ns=0; ns<nb_subblocks; ns++)//compute block by block
{
//mean and variance of real and imaginary parts of emitted symbols
compute_symb_stats(Es, Vs, ns, select_half, apriori_data,
re_part, im_part, re_bin_part, im_bin_part);
//find equivalent channel matrix
EquivCh(H_eq, c_impulse_response.get_col(ns));
//compute invariant inverse
K = H_eq*itpp::diag(Vs)*H_eq.transpose()+sigma2*itpp::eye(2*block_duration*nb_rec_ant);
//find equivalent received signal
EquivRecSig(x_eq, rec_sig(ns*block_duration,(ns+1)*block_duration-1,0,nb_rec_ant-1));
for (q=0; q<symbols_block; q++)//symbols/block
{
//compute the extrinsic information of coded bits
//real part
//IC + filtering (real and imaginary parts of one symbol)
interf = Es;
interf(2*q) = 0;//this is the symbol to recover
temp = H_eq.get_col(2*q);
w = temp/(temp*temp);//filter impulse response
s_tilde = w*(x_eq-H_eq*interf);
mu_res = w*temp;//mean of the filtered signal
index = select_half*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
for (k=0; k<half_nb_bits_symb; k++)
{
nom = -INFINITY;
denom = -INFINITY;
for (s=0; s<half_const_len; s++)
{
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);
tmp = -itpp::sqr(s_tilde-mu_res*re_part(s))/(2*sigma2_res)+ \
itpp::to_vec(re_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
if (re_bin_part(s,k))
nom = std::max(nom, tmp);
else
denom = std::max(denom, tmp);
}
extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
}
//end real part
//imaginary part
//IC + filtering (real and imaginary parts of one symbol)
interf = Es;
interf(2*q+1) = 0;//this is the symbol to recover
temp = H_eq.get_col(2*q+1);
w = temp/(temp*temp);//filter impulse response
s_tilde = w*(x_eq-H_eq*interf);
mu_res = w*temp;//mean of the filtered signal
index = (1-select_half)*half_nb_bits_symb+nb_bits_symb*q+ns*nb_bits_subblock;
for (k=0; k<half_nb_bits_symb; k++)
{
nom = -INFINITY;
denom = -INFINITY;
for (s=0; s<half_const_len; s++)
{
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);
tmp = -itpp::sqr(s_tilde-mu_res*im_part(s))/(2*sigma2_res)+ \
itpp::to_vec(im_bin_part.get_row(s))*apriori_data.mid(index, half_nb_bits_symb);
if (im_bin_part(s,k))
nom = std::max(nom, tmp);
else
denom = std::max(denom, tmp);
}
extrinsic_data(index+k) = (nom-denom)-apriori_data(index+k);
}
//end imaginary part
}//symbols/block
}//block by block
}
void SISO::Alamouti_maxlogMAP(itpp::vec &extrinsic_data, const itpp::cmat &rec_sig, const itpp::vec &apriori_data)
//maxlogMAP algorithm for Alamouti ST code
{
//matched filter
int int_len = apriori_data.length();//interleaver length
int nb_symb = (int)(int_len/nb_bits_symb);//number of symbols/block
itpp::cvec comb_sig(nb_symb);
comb_sig.zeros();
itpp::cmat conj_H = itpp::conj(c_impulse_response);
itpp::cmat conj_X = itpp::conj(rec_sig);
register int nr,n,cs;
for (nr=0; nr<nb_rec_ant; nr++)
{
for (n=0; n<(nb_symb/2); n++)
{
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));
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));
}
}
//extrinsic information of coded bits
int const_size = itpp::pow2i(nb_bits_symb);//constellation size
double buffer;
double nom,denom;
double temp;
int index;
extrinsic_data.set_size(nb_bits_symb*nb_symb);
for (n=0; n<nb_symb; n++)
{
buffer = itpp::sum_sqr(itpp::abs(c_impulse_response.get_col(n/2)));
for (nr=0; nr<nb_bits_symb; nr++)
{
nom = -INFINITY;
denom = -INFINITY;
for (cs=0; cs<const_size; cs++)
{
temp = -itpp::sqr(comb_sig(n)-buffer*constellation(cs))/(2*buffer*sigma2)+ \
itpp::to_vec(bin_constellation.get_row(cs))*apriori_data.mid(n*nb_bits_symb, nb_bits_symb);
if (bin_constellation(cs,nr))
nom = std::max(nom, temp);
else
denom = std::max(denom, temp);
}
index = n*nb_bits_symb+nr;
extrinsic_data(index) = (nom-denom)-apriori_data(index);//extrinsic information
}
}
}
void SISO::demodulator_logMAP(itpp::vec &extrinsic_data, const itpp::cvec &rec_sig, const itpp::vec &apriori_data)
{
int nb_symb = rec_sig.length();
int const_size = itpp::pow2i(nb_bits_symb);
double nom,denom,temp;
register int k,i,cs;
int index;
extrinsic_data.set_size(nb_bits_symb*nb_symb);
for (k=0; k<nb_symb; k++)
{
for (i=0; i<nb_bits_symb; i++)
{
nom = 0;
denom = 0;
for (cs=0; cs<const_size; cs++)
{
temp = -itpp::sqr(rec_sig(k)-c_impulse_response(0,k)*constellation(cs))/(2*sigma2)+\
itpp::to_vec(bin_constellation.get_row(cs))*apriori_data.mid(k*nb_bits_symb, nb_bits_symb);
if (bin_constellation(cs,i))
nom += std::exp(temp);
else
denom += std::exp(temp);
}
index = k*nb_bits_symb+i;
extrinsic_data(index) = std::log(nom/denom)-apriori_data(index);//extrinsic information
}
}
}
void SISO::demodulator_maxlogMAP(itpp::vec &extrinsic_data, const itpp::cvec &rec_sig, const itpp::vec &apriori_data)
{
int nb_symb = rec_sig.length();
int const_size = itpp::pow2i(nb_bits_symb);
double nom,denom,temp;
register int k,i,cs;
int index;
extrinsic_data.set_size(nb_bits_symb*nb_symb);
for (k=0; k<nb_symb; k++)
{
for (i=0; i<nb_bits_symb; i++)
{
nom = -INFINITY;
denom = -INFINITY;
for (cs=0; cs<const_size; cs++)
{
temp = -itpp::sqr(rec_sig(k)-c_impulse_response(0,k)*constellation(cs))/(2*sigma2)+\
itpp::to_vec(bin_constellation.get_row(cs))*apriori_data.mid(k*nb_bits_symb, nb_bits_symb);
if (bin_constellation(cs,i))
nom = std::max(nom, temp);
else
denom = std::max(denom, temp);
}
index = k*nb_bits_symb+i;
extrinsic_data(index) = (nom-denom)-apriori_data(index);//extrinsic information
}
}
}
}//end namespace tr
itpp::Mat::transpose
Mat< Num_T > transpose() const
Matrix transpose.
Definition mat.h:1169
itpp::diag
Mat< T > diag(const Vec< T > &v, const int K=0)
Create a diagonal matrix using vector v as its diagonal.
Definition matfunc.h:557
itpp::inv
bool inv(const mat &X, mat &Y)
Inverse of real square matrix.
Definition inv.cpp:87
itpp::pow2i
int pow2i(int x)
Calculate two to the power of x (2^x); x is integer.
Definition log_exp.h:53
itpp::exp
vec exp(const vec &x)
Exp of the elements of a vector x.
Definition log_exp.h:155
itpp::sum
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition matfunc.h:59
itpp::transpose
void transpose(const Mat< T > &m, Mat< T > &out)
Transposition of the matrix m returning the transposed matrix in out.
Definition matfunc.h:308
itpp::sum_sqr
T sum_sqr(const Vec< T > &v)
Sum of square of the elements in a vector.
Definition matfunc.h:116
itpp::prod
T prod(const Vec< T > &v)
The product of all elements in the vector.
Definition matfunc.h:195
itpp::imag
vec imag(const cvec &data)
Imaginary part of complex values.
Definition elem_math.cpp:180
itpp::find
ivec find(const bvec &invector)
Return a integer vector with indicies where bvec == 1.
Definition specmat.cpp:40
itpp::real
vec real(const cvec &data)
Real part of complex values.
Definition elem_math.cpp:157
itpp::sqr
vec sqr(const cvec &data)
Absolute square of elements.
Definition elem_math.cpp:36
itpp::conj
cvec conj(const cvec &x)
Conjugate of complex value.
Definition elem_math.h:226
itpp::reshape
Mat< T > reshape(const Mat< T > &m, int rows, int cols)
Reshape the matrix into an rows*cols matrix.
Definition matfunc.h:822
itpp
itpp namespace
Definition itmex.h:37
itpp::eye
template void eye(int, mat &)
Template instantiation of eye.
itpp::bin2dec
ITPP_EXPORT int bin2dec(const bvec &inbvec, bool msb_first=true)
Convert a bvec to decimal int with the first bit as MSB if msb_first == true.
itpp::to_vec
vec to_vec(const Vec< T > &v)
Converts a Vec<T> to vec.
Definition converters.h:93
itpp::dec2bin
ITPP_EXPORT bvec dec2bin(int length, int index)
Convert a decimal int index to bvec using length bits in the representation.
itpp::elem_div
Mat< Num_T > elem_div(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise division of two matrices.
Definition mat.h:1688
itpp::elem_mult
Mat< Num_T > elem_mult(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise multiplication of two matrices.
Definition mat.h:1582
itpp::abs
bin abs(const bin &inbin)
absolute value of bin
Definition binary.h:174
itpp::outer_product
Mat< Num_T > outer_product(const Vec< Num_T > &v1, const Vec< Num_T > &v2, bool hermitian=false)
Outer product of two vectors v1 and v2.
Definition vec.h:1021
siso.h
Definitions for Soft Input Soft Output (SISO) modules class.

Generated on Mon Apr 22 2024 09:45:27 for IT++ by Doxygen 1.9.8