IT++ Logo
siso_mud.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/comm/siso.h>
00030 #include <limits>
00031 #ifndef INFINITY
00032 #define INFINITY std::numeric_limits<double>::infinity()
00033 #endif
00034 
00035 namespace itpp
00036 {
00037 void SISO::descrambler(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data, const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data)
00038 /*
00039   inputs:
00040   intrinsic_coded - intrinsic information of coded bits (repetition code output)
00041   apriori_data - a priori information of informational bits (repetition code input)
00042   outputs:
00043   extrinsic_coded - extrinsic information of coded bits
00044   extrinsic_data - Logarithm of Likelihood Ratio of informational bits
00045 */
00046 {
00047     //get parameters
00048     int nb_bits = apriori_data.length();
00049     int Nc = scrambler_pattern.length();
00050     //implementation
00051     extrinsic_data.set_size(nb_bits);
00052     extrinsic_coded.set_size(nb_bits*Nc);
00053     register int n,k;
00054 #pragma omp parallel for private(n,k)
00055     for (k=0; k<nb_bits; k++)
00056     {
00057         extrinsic_data[k] = 0;//apriori_data[k];//add a priori info
00058         for (n=0; n<Nc; n++)
00059         {
00060             extrinsic_data[k] += (1-2*double(scrambler_pattern[n]))*intrinsic_coded[n+k*Nc];
00061         }
00062         for (n=0; n<Nc; n++)
00063         {
00064             extrinsic_coded[n+k*Nc] = (1-2*double(scrambler_pattern[n]))*extrinsic_data[k]-intrinsic_coded[n+k*Nc];
00065         }
00066     }
00067 }
00068 
00069 void SISO::zpFIRfilter(itpp::vec &filt, const itpp::vec &h, const itpp::vec &sig)
00070 //FIR filter for a zero padded signal (L zeros are added at the end of the signal)
00071 {
00072     //get parameters
00073     int L = h.length()-1;
00074     int N = sig.length();
00075     //implementation
00076     register int n,l;
00077 #pragma omp parallel for private(n,l)
00078     for (n=0; n<(N+L); n++)
00079     {
00080         filt[n] = 0;
00081         for (l=0; l<=L; l++)
00082         {
00083             if ((n-l)<0)
00084             {
00085                 break;//channel has state 0 at the beginning
00086             }
00087             if ((n-l)>=N)
00088             {
00089                 continue;//channel has state 0 in the end
00090             }
00091             filt[n] += (h[l]*sig[n-l]);
00092         }
00093     }
00094 }
00095 
00096 void SISO::gen_hyperTrellis(void)
00097 /* generates channel hyper trellis for binary symbols
00098  * the channel is a MISO system
00099  * BPSK mapping: 0->+1, 1->-1
00100  */
00101 {
00102     //get parameters
00103     int nb_usr = impulse_response.rows();
00104     int ch_order = impulse_response.cols()-1;
00105     int p_order = prec_gen.length()-1;
00106     int max_order = std::max(ch_order, p_order);
00107 
00108     //initialize hypertrellis
00109     chtrellis.numInputSymbols = itpp::pow2i(nb_usr);
00110     int mem_len = nb_usr*max_order;
00111     if (mem_len>=(int)(8*sizeof(int)))
00112     {
00113         std::string msg = "SISO::gen_hyperTrellis: memory length of the hyperchannel should be at most ";
00114         msg += itpp::to_str(8*sizeof(int)-1);
00115         msg += ". Try to lower the number of users, channel order or the order of the precoding polynomial (if any)";
00116         print_err_msg(msg);
00117         return;
00118     }
00119     chtrellis.stateNb = itpp::pow2i(mem_len);
00120     try
00121     {
00122         unsigned int len =  static_cast<unsigned int>(chtrellis.stateNb)*static_cast<unsigned int>(chtrellis.numInputSymbols);
00123         chtrellis.nextState = new int[len];
00124         chtrellis.prevState = new int[len];
00125         chtrellis.output = new double[len];
00126         chtrellis.input = new int[len];
00127     } catch (std::bad_alloc)
00128     {
00129         std::string msg = "SISO::gen_hyperTrellis: not enough memory for the channel trellis variables. The number of states is ";
00130         msg += itpp::to_str(chtrellis.stateNb);
00131         msg += " and the number of input symbols ";
00132         msg += itpp::to_str(chtrellis.numInputSymbols);
00133         print_err_msg(msg);
00134         return;
00135     }
00136     itpp::ivec index(chtrellis.stateNb);
00137     index.zeros();
00138     itpp::bvec hyper_ch_mem(mem_len);
00139     itpp::bvec hyper_ch_in(nb_usr);
00140     itpp::bvec hyper_states(mem_len);
00141     itpp::bin feedback;
00142 
00143     //create hypertrellis
00144     register int n,k,p,r;
00145     int buffer;
00146     double hyper_ch_out;
00147     for (k=0; k<chtrellis.stateNb; k++)
00148     {
00149         hyper_ch_mem = itpp::dec2bin(mem_len, k);//initial state
00150         for (n=0; n<chtrellis.numInputSymbols; n++)
00151         {
00152             hyper_ch_in = itpp::dec2bin(nb_usr, n);//MISO channel input
00153             hyper_ch_out = 0;
00154             for (r=0; r<nb_usr; r++)
00155             {
00156                 buffer = r*max_order;
00157                 //precoder
00158                 feedback = hyper_ch_in[r];
00159                 for (p=1; p<=p_order; p++)
00160                 {
00161                     if (prec_gen(p))
00162                     {
00163                         feedback ^= hyper_ch_mem[p-1+buffer];//xor
00164                     }
00165                 }
00166                 //FIR channel output
00167                 hyper_ch_out += (1-2*double(feedback))*impulse_response(r,0);
00168                 for (p=0; p<ch_order; p++)
00169                 {
00170                     hyper_ch_out += (1-2*double(hyper_ch_mem[p+buffer]))*impulse_response(r,p+1);//hyper channel output for user r
00171                 }
00172                 //(equivalent) channel next state
00173                 hyper_states[buffer] = feedback;
00174                 for (p=0; p<(max_order-1); p++)
00175                 {
00176                     hyper_states[p+buffer+1] = hyper_ch_mem[p+buffer];//next hyper state for user r
00177                 }
00178             }
00179             chtrellis.output[k+n*chtrellis.stateNb] = hyper_ch_out;
00180             buffer = itpp::bin2dec(hyper_states);//next state from an initial state and a given input
00181             chtrellis.nextState[k+n*chtrellis.stateNb] = buffer;
00182             chtrellis.prevState[buffer+index[buffer]*chtrellis.stateNb] = k;
00183             chtrellis.input[buffer+index[buffer]*chtrellis.stateNb] = n;
00184             index[buffer]++;
00185         }
00186     }
00187 }
00188 
00190 
00193 void SISO::mud_maxlogMAP(itpp::mat &extrinsic_data, const itpp::vec &rec_sig, const itpp::mat &apriori_data)
00194 /* output:
00195  * extrinsic_data - extrinsic information for the chips (usr_nb x block_len)
00196  * inputs:
00197  * rec_sig - received signal (1 x block_len)
00198  * apriori_data - a priori information for the chips (usr_nb x block_len)
00199  */
00200 {
00201     //get parameters
00202     int nb_usr = apriori_data.rows();
00203     int block_len = apriori_data.cols();
00204 
00205     //init trellis
00206     gen_hyperTrellis();
00207 
00208     //initial conditions for A = log(alpha) and B = log(beta)
00209     double *A = NULL,*B = NULL;
00210     try
00211     {
00212         A = new double[chtrellis.stateNb*(block_len+1)];
00213         B = new double[chtrellis.stateNb*(block_len+1)];
00214     } catch (std::bad_alloc)
00215     {
00216         std::string msg = "SISO::mud_maxlogMAP: Not enough memory for alphas and betas. The number of states is ";
00217         msg += itpp::to_str(chtrellis.stateNb);
00218         msg += " and the block length ";
00219         msg += itpp::to_str(block_len);
00220         print_err_msg(msg);
00221     }
00222     register int n;
00223     A[0] = 0;
00224     B[block_len*chtrellis.stateNb] = 0;
00225     double buffer = (tail?-INFINITY:0);
00226 #pragma omp parallel for private(n)
00227     for (n=1; n<chtrellis.stateNb; n++)
00228     {
00229         A[n] = -INFINITY;
00230         B[n+block_len*chtrellis.stateNb] = buffer;//if tail==false the final state is not known
00231     }
00232 
00233     //compute log(alpha) (forward recursion)
00234     register int s,k;
00235     int sp,i;
00236     itpp::bvec in_chips(nb_usr);
00237 #pragma omp parallel sections private(n,buffer,s,k,sp,in_chips)
00238     {
00239         for (n=1; n<=block_len; n++)
00240         {
00241             buffer = -INFINITY;//normalization factor
00242             for (s=0; s<chtrellis.stateNb; s++)
00243             {
00244                 A[s+n*chtrellis.stateNb] = -INFINITY;
00245                 for (k=0; k<chtrellis.numInputSymbols; k++)
00246                 {
00247                     sp = chtrellis.prevState[s+k*chtrellis.stateNb];
00248                     i = chtrellis.input[s+k*chtrellis.stateNb];
00249                     in_chips = itpp::dec2bin(nb_usr, i);
00250                     A[s+n*chtrellis.stateNb] = std::max(A[s+n*chtrellis.stateNb], \
00251                                                         A[sp+(n-1)*chtrellis.stateNb]-itpp::sqr(rec_sig[n-1]-chtrellis.output[sp+i*chtrellis.stateNb])/(2*sigma2)+\
00252                                                         itpp::to_vec(in_chips)*apriori_data.get_col(n-1));
00253                 }
00254                 buffer = std::max(buffer, A[s+n*chtrellis.stateNb]);
00255             }
00256             //normalization
00257             for (s=0; s<chtrellis.stateNb; s++)
00258             {
00259                 A[s+n*chtrellis.stateNb] -= buffer;
00260             }
00261         }
00262 
00263         //compute log(beta) (backward recursion)
00264 #pragma omp section
00265         for (n=block_len-1; n>=0; n--)
00266         {
00267             buffer = -INFINITY;//normalization factor
00268             for (s=0; s<chtrellis.stateNb; s++)
00269             {
00270                 B[s+n*chtrellis.stateNb] = -INFINITY;
00271                 for (k=0; k<chtrellis.numInputSymbols; k++)
00272                 {
00273                     sp = chtrellis.nextState[s+k*chtrellis.stateNb];
00274                     in_chips = itpp::dec2bin(nb_usr, k);
00275                     B[s+n*chtrellis.stateNb] = std::max(B[s+n*chtrellis.stateNb], \
00276                                                         B[sp+(n+1)*chtrellis.stateNb]-itpp::sqr(rec_sig[n]-chtrellis.output[s+k*chtrellis.stateNb])/(2*sigma2)+\
00277                                                         itpp::to_vec(in_chips)*apriori_data.get_col(n));
00278                 }
00279                 buffer = std::max(buffer, B[s+n*chtrellis.stateNb]);
00280             }
00281             //normalization
00282             for (s=0; s<chtrellis.stateNb; s++)
00283             {
00284                 B[s+n*chtrellis.stateNb] -= buffer;
00285             }
00286         }
00287     }
00288 
00289     //compute extrinsic information
00290     double nom, denom;
00291     extrinsic_data.set_size(nb_usr,block_len);
00292     register int u;
00293 #pragma omp parallel for private(u,n,s,k,nom,denom,in_chips,buffer)
00294     for (u=0; u<nb_usr; u++)
00295     {
00296         for (n=1; n<=block_len; n++)
00297         {
00298             nom = -INFINITY;
00299             denom = -INFINITY;
00300             for (s=0; s<chtrellis.stateNb; s++)
00301             {
00302                 for (k=0; k<chtrellis.numInputSymbols; k++)
00303                 {
00304                     in_chips = itpp::dec2bin(nb_usr, k);
00305                     buffer = A[s+(n-1)*chtrellis.stateNb]+B[chtrellis.nextState[s+k*chtrellis.stateNb]+n*chtrellis.stateNb]-\
00306                              itpp::sqr(rec_sig[n-1]-chtrellis.output[s+k*chtrellis.stateNb])/(2*sigma2)+\
00307                              itpp::to_vec(in_chips)*apriori_data.get_col(n-1);
00308                     if (in_chips[u])
00309                     {
00310                         nom = std::max(nom, buffer);
00311                     }
00312                     else
00313                     {
00314                         denom = std::max(denom, buffer);
00315                     }
00316                 }
00317             }
00318             extrinsic_data(u,n-1) = (nom-denom)-apriori_data(u,n-1);
00319         }
00320     }
00321     //free memory
00322     delete[] chtrellis.nextState;
00323     delete[] chtrellis.prevState;
00324     delete[] chtrellis.output;
00325     delete[] chtrellis.input;
00326     delete[] A;
00327     delete[] B;
00328 }
00329 
00331 
00333 void SISO::GCD(itpp::mat &extrinsic_data, const itpp::vec &rec_sig, const itpp::mat &apriori_data)
00334 /* Gaussian Chip Detector
00335  * output:
00336  * extrinsic_data - extrinsic information of emitted chips
00337  * inputs:
00338  * rec_sig - received signal
00339  * apriori_data - a priori information of emitted chips
00340  */
00341 {
00342     //get parameters
00343     int N = apriori_data.cols();//emitted frames of non-zero samples
00344     int K = apriori_data.rows();//number of users
00345     int L = impulse_response.cols()-1;//channel order
00346     //other parameters
00347     register int n,k;
00348 
00349     //mean and variance of each chip (NxK)
00350     itpp::mat Ex = -itpp::tanh(apriori_data/2.0);//take into account BPSK mapping
00351     itpp::mat Vx = 1.0-sqr(Ex);
00352 
00353     //expectation and variance of the received signal
00354     itpp::vec Er(N+L);
00355     Er.zeros();
00356     itpp::mat Covr;
00357     try
00358     {
00359         Covr.set_size(N+L,N+L);
00360     } catch (std::bad_alloc)
00361     {
00362         std::string msg = "SISO::GCD: not enough memory when allocating space for the covariance matrix. The interleaver length is ";
00363         msg += itpp::to_str(N);
00364         msg += ". Use sGCD instead or try to reduce the interleaver length";
00365         print_err_msg(msg);
00366         return;
00367     }
00368     //create eye function
00369     Covr.zeros();
00370     for (n=0; n<(N+L); n++)
00371         Covr(n,n) = 1;
00372     itpp::vec col(N+L);
00373     col.zeros();
00374     itpp::vec row(N);
00375     row.zeros();
00376     itpp::mat h_eq(N+L,N);
00377     for (k=0; k<K; k++)
00378     {
00379         col.replace_mid(0, impulse_response.get_row(k));
00380         row(0) = impulse_response(k,0);
00381         h_eq = itpp::toeplitz(col, row);
00382         Er += h_eq*Ex.get_row(k);
00383         Covr += (h_eq*itpp::diag(Vx.get_row(k)))*h_eq.T();
00384     }
00385 
00386     //extrinsic information
00387     double nom,denom;
00388     Er = rec_sig-Er;
00389     itpp::mat inv_Covr(N+L,N+L);
00390     inv_Covr = itpp::inv(Covr);
00391     itpp::mat inv_cov_zeta(N+L,N+L);
00392     itpp::vec rec_chip(N+L);
00393     extrinsic_data.set_size(K,N);
00394     for (k=0; k<K; k++)
00395     {
00396 #pragma omp parallel for private(n,inv_cov_zeta,rec_chip,nom,denom) firstprivate(col)
00397         for (n=0; n<N; n++)
00398         {
00399             col.replace_mid(n, impulse_response.get_row(k));
00400             inv_cov_zeta = inv_Covr+itpp::outer_product(inv_Covr*col, inv_Covr.T()*(col*Vx(k,0)))/(1-(col*Vx(k,0))*(inv_Covr*col));
00401             rec_chip = Er-col*(+1-Ex(k,n));
00402             nom = -0.5*rec_chip*(inv_cov_zeta*rec_chip);
00403             rec_chip = Er-col*(-1-Ex(k,n));
00404             denom = -0.5*rec_chip*(inv_cov_zeta*rec_chip);
00405             extrinsic_data(k,n) = denom-nom;//take into account BPSK mapping
00406             col(n) = 0;
00407         }
00408     }
00409 }
00410 
00412 
00415 void SISO::sGCD(itpp::mat &extrinsic_data, const itpp::vec &rec_sig, const itpp::mat &apriori_data)
00416 /* simplified GCD
00417  * output:
00418  * extrinsic_data - extrinsic information of emitted chips
00419  * inputs:
00420  * rec_sig - received signal
00421  * apriori_data - a priori information of emitted chips
00422  */
00423 {
00424     //get parameters
00425     int N = apriori_data.cols();//emitted frames of non-zero samples
00426     int K = apriori_data.rows();//number of users
00427     int L = impulse_response.cols()-1;//channel order
00428     //other parameters
00429     register int n,k;
00430 
00431     //mean and variance of each chip (NxK)
00432     itpp::mat Ex = -itpp::tanh(apriori_data/2.0);//take into account BPSK mapping
00433     itpp::mat Vx = 1.0-itpp::sqr(Ex);
00434 
00435     //mean and variance for the samples of the received signal
00436     itpp::vec Er(N+L);
00437     Er.zeros();
00438     itpp::vec Vr = sigma2*itpp::ones(N+L);
00439     itpp::vec buffer(N+L);
00440     for (k=0; k<K; k++)
00441     {
00442         zpFIRfilter(buffer, impulse_response.get_row(k), Ex.get_row(k));
00443         Er += buffer;
00444         zpFIRfilter(buffer, itpp::sqr(impulse_response.get_row(k)), Vx.get_row(k));
00445         Vr += buffer;
00446     }
00447 
00448     //extrinsic information for the samples of the received signal
00449     Er = rec_sig-Er;
00450     itpp::vec ch(L+1);
00451     extrinsic_data.set_size(K,N);
00452     for (k=0; k<K; k++)
00453     {
00454         ch = impulse_response.get_row(k);
00455 #pragma omp parallel for private(n)
00456         for (n=0; n<N; n++)
00457         {
00458             extrinsic_data(k,n) = -2*itpp::elem_div(ch, Vr.mid(n,L+1)-sqr(ch)*Vx(k,n))*(Er.mid(n,L+1)+ch*Ex(k,n));//take into account BPSK mapping
00459         }
00460     }
00461 }
00462 }//end namespace tr
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:32 for IT++ by Doxygen 1.7.4