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