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::gen_chtrellis(void) 00038 // generate trellis for precoded FIR channels with real coefficients 00039 { 00040 //get parameters 00041 int mem_len = impulse_response.cols()-1;//memory length of the channel 00042 int p_order = prec_gen.length()-1;//order of the precoder polynomial 00043 00044 //other variables 00045 register int n,k,j; 00046 double inputs[] = {1.0,-1.0};//1->-1, 0->+1 00047 int index; 00048 double feedback[2]; 00049 00050 //create channel trellis 00051 int equiv_ch_mem_len = std::max(mem_len, p_order); 00052 chtrellis.stateNb = (1<<equiv_ch_mem_len); 00053 chtrellis.output = new double[2*chtrellis.stateNb]; 00054 chtrellis.nextState = new int[2*chtrellis.stateNb]; 00055 chtrellis.prevState = new int[2*chtrellis.stateNb]; 00056 chtrellis.input = new int[2*chtrellis.stateNb]; 00057 00058 //initialize trellis 00059 itpp::ivec enc_mem(equiv_ch_mem_len); 00060 #pragma omp parallel for private(n,enc_mem,k,feedback,j) 00061 for (n=0; n<chtrellis.stateNb; n++) //initial state 00062 { 00063 enc_mem = 1-2*itpp::to_ivec(itpp::dec2bin(equiv_ch_mem_len, n));//1->-1, 0->+1 00064 //output 00065 for (k=0; k<2; k++) 00066 { 00067 feedback[k] = inputs[k]; 00068 for (j=1; j<=p_order; j++) 00069 if (prec_gen[j]) 00070 feedback[k] *= enc_mem[j-1];//xor truth table must remain the same 00071 chtrellis.output[n+k*chtrellis.stateNb] = feedback[k]*impulse_response(0,0); 00072 for (j=1; j<=mem_len; j++) 00073 chtrellis.output[n+k*chtrellis.stateNb] += (enc_mem[j-1]*impulse_response(0,j)); 00074 } 00075 //next state 00076 for (j=(equiv_ch_mem_len-1); j>0; j--) 00077 enc_mem[j] = enc_mem[j-1]; 00078 for (k=0; k<2; k++) 00079 { 00080 enc_mem[0] = int(feedback[k]); 00081 chtrellis.nextState[n+k*chtrellis.stateNb] = itpp::bin2dec(itpp::to_bvec((1-enc_mem)/2), true);//-1->1, +1->0 00082 } 00083 } 00084 00085 #pragma omp parallel for private(j,index,n,k) 00086 for (j=0; j<chtrellis.stateNb; j++) 00087 { 00088 index = 0; 00089 for (n=0; n<chtrellis.stateNb; n++) 00090 { 00091 for (k=0; k<2; k++) 00092 { 00093 if (chtrellis.nextState[n+k*chtrellis.stateNb]==j) 00094 { 00095 chtrellis.prevState[j+index*chtrellis.stateNb] = n; 00096 chtrellis.input[j+index*chtrellis.stateNb] = k;//this is an index to the channel input 00097 index++; 00098 } 00099 } 00100 } 00101 } 00102 } 00103 00104 void SISO::equalizer_logMAP(itpp::vec &extrinsic_data, const itpp::vec &rec_sig, const itpp::vec &apriori_data) 00105 /* 00106 extrinsic_data - extrinsic information of channel input symbols 00107 rec_sig - received symbols 00108 apriori_data - a priori information of channel input symbols 00109 */ 00110 { 00111 //get parameters 00112 int N = rec_sig.length();//length of the received frame 00113 //other parameters 00114 register int n,k,m; 00115 int pstates[2]; 00116 int nstates[2]; 00117 int inputs[2]; 00118 double C[2];//log(gamma) 00119 double sum; 00120 double sumbis; 00121 double buffer; 00122 00123 //initialize trellis 00124 gen_chtrellis(); 00125 //initialize log(alpha) and log(beta) 00126 double* A = new double[chtrellis.stateNb*(N+1)]; 00127 double* B = new double[chtrellis.stateNb*(N+1)]; 00128 A[0] = 0; 00129 B[N*chtrellis.stateNb] = 0; 00130 sum = (tail?-INFINITY:0); 00131 #pragma omp parallel for private(n) 00132 for (n=1; n<chtrellis.stateNb; n++) 00133 { 00134 A[n] = -INFINITY; 00135 B[n+N*chtrellis.stateNb] = sum;//if tail==false the final state is not known 00136 } 00137 00138 #pragma omp parallel sections private(n,sum,m,k,C) 00139 { 00140 //forward recursion 00141 for (n=1; n<=N; n++) 00142 { 00143 sum = 0;//normalisation factor 00144 for (m=0; m<chtrellis.stateNb; m++) //final state 00145 { 00146 for (k=0; k<2; k++) 00147 { 00148 pstates[k] = chtrellis.prevState[m+chtrellis.stateNb*k];//determine previous states 00149 inputs[k] = chtrellis.input[m+chtrellis.stateNb*k];//determine input 00150 C[k] = (inputs[k])*apriori_data[n-1]-itpp::sqr(rec_sig[n-1]-chtrellis.output[pstates[k]+chtrellis.stateNb*inputs[k]])/(2*sigma2);//compute log of gamma 00151 } 00152 A[m+n*chtrellis.stateNb] = itpp::log_add(A[pstates[0]+(n-1)*chtrellis.stateNb]+C[0], A[pstates[1]+(n-1)*chtrellis.stateNb]+C[1]); 00153 sum += std::exp(A[m+n*chtrellis.stateNb]); 00154 } 00155 //normalisation 00156 sum = std::log(sum); 00157 for (m=0; m<chtrellis.stateNb; m++) 00158 { 00159 A[m+n*chtrellis.stateNb] -= sum; 00160 } 00161 } 00162 00163 //backward recursion 00164 #pragma omp section 00165 for (n=N-1; n>=0; n--) 00166 { 00167 sum = 0;//normalisation factor 00168 for (m=0; m<chtrellis.stateNb; m++) //initial state 00169 { 00170 for (k=0; k<2; k++) 00171 { 00172 nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states 00173 C[k] = (k)*apriori_data[n]-itpp::sqr(rec_sig[n]-chtrellis.output[m+k*chtrellis.stateNb])/(2*sigma2);//compute log of gamma 00174 } 00175 B[m+n*chtrellis.stateNb] = itpp::log_add(B[nstates[0]+(n+1)*chtrellis.stateNb]+C[0], B[nstates[1]+(n+1)*chtrellis.stateNb]+C[1]); 00176 sum += std::exp(B[m+n*chtrellis.stateNb]); 00177 } 00178 //normalisation 00179 sum = std::log(sum); 00180 for (m=0; m<chtrellis.stateNb; m++) 00181 { 00182 B[m+n*chtrellis.stateNb] -= sum; 00183 } 00184 } 00185 } 00186 00187 //compute extrinsic_data 00188 extrinsic_data.set_size(N); 00189 #pragma omp parallel for private(n,sum,sumbis,m,k,nstates,C,buffer) 00190 for (n=1; n<=N; n++) 00191 { 00192 sum = 0;//could be replaced by a vector 00193 sumbis = 0; 00194 for (m=0; m<chtrellis.stateNb; m++) //initial state 00195 { 00196 for (k=0; k<2; k++) //input index 00197 { 00198 nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states 00199 C[k] = (k)*apriori_data[n-1]-itpp::sqr(rec_sig[n-1]-chtrellis.output[m+k*chtrellis.stateNb])/(2*sigma2);//compute log of gamma 00200 buffer = std::exp(A[m+(n-1)*chtrellis.stateNb]+C[k]+B[nstates[k]+n*chtrellis.stateNb]); 00201 if (k) 00202 sum += buffer;//1 00203 else 00204 sumbis += buffer;//0 00205 } 00206 } 00207 extrinsic_data[n-1] = std::log(sum/sumbis)-apriori_data[n-1]; 00208 } 00209 00210 //free memory 00211 delete[] chtrellis.output; 00212 delete[] chtrellis.nextState; 00213 delete[] chtrellis.prevState; 00214 delete[] chtrellis.input; 00215 delete[] A; 00216 delete[] B; 00217 } 00218 00219 void SISO::equalizer_maxlogMAP(itpp::vec &extrinsic_data, const itpp::vec &rec_sig, const itpp::vec &apriori_data) 00220 /* 00221 extrinsic_data - extrinsic information for channel input symbols 00222 rec_sig - received symbols 00223 apriori_data - a priori information for channel input symbols 00224 */ 00225 { 00226 //get parameters 00227 int N = rec_sig.length();//length of the received frame 00228 //other parameters 00229 register int n,k,m; 00230 int pstates[2]; 00231 int nstates[2]; 00232 int inputs[2]; 00233 double C[2];//log(gamma) 00234 double sum; 00235 double sumbis; 00236 double buffer; 00237 00238 //initialize trellis 00239 gen_chtrellis(); 00240 //initialize log(alpha) and log(beta) 00241 double* A = new double[chtrellis.stateNb*(N+1)]; 00242 double* B = new double[chtrellis.stateNb*(N+1)]; 00243 A[0] = 0; 00244 for (n=1; n<chtrellis.stateNb; n++) 00245 A[n] = -INFINITY; 00246 B[N*chtrellis.stateNb] = 0; 00247 sum = (tail?-INFINITY:0); 00248 for (n=1; n<chtrellis.stateNb; n++) 00249 B[n+N*chtrellis.stateNb] = sum;//if tail==false the final state is not known 00250 00251 #pragma omp parallel sections private(n,sum,m,k,C) 00252 { 00253 //forward recursion 00254 for (n=1; n<=N; n++) 00255 { 00256 sum = -INFINITY;//normalization factor 00257 for (m=0; m<chtrellis.stateNb; m++) //final state 00258 { 00259 for (k=0; k<2; k++) 00260 { 00261 pstates[k] = chtrellis.prevState[m+chtrellis.stateNb*k];//determine previous states 00262 inputs[k] = chtrellis.input[m+chtrellis.stateNb*k];//determine input 00263 C[k] = (inputs[k])*apriori_data[n-1]-itpp::sqr(rec_sig[n-1]-chtrellis.output[pstates[k]+chtrellis.stateNb*inputs[k]])/(2*sigma2);//compute log of gamma 00264 } 00265 A[m+n*chtrellis.stateNb] = std::max(A[pstates[0]+(n-1)*chtrellis.stateNb]+C[0], A[pstates[1]+(n-1)*chtrellis.stateNb]+C[1]); 00266 sum = std::max(sum, A[m+n*chtrellis.stateNb]); 00267 } 00268 for (m=0; m<chtrellis.stateNb; m++) 00269 A[m+n*chtrellis.stateNb] -= sum; 00270 } 00271 00272 //backward recursion 00273 #pragma omp section 00274 for (n=N-1; n>=0; n--) 00275 { 00276 sum = -INFINITY;//normalisation factor 00277 for (m=0; m<chtrellis.stateNb; m++) //initial state 00278 { 00279 for (k=0; k<2; k++) 00280 { 00281 nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states 00282 C[k] = (k)*apriori_data[n]-itpp::sqr(rec_sig[n]-chtrellis.output[m+k*chtrellis.stateNb])/(2*sigma2);//compute log of gamma 00283 } 00284 B[m+n*chtrellis.stateNb] = std::max(B[nstates[0]+(n+1)*chtrellis.stateNb]+C[0], B[nstates[1]+(n+1)*chtrellis.stateNb]+C[1]); 00285 sum = std::max(sum, B[m+n*chtrellis.stateNb]); 00286 } 00287 for (m=0; m<chtrellis.stateNb; m++) 00288 B[m+n*chtrellis.stateNb] -= sum; 00289 } 00290 } 00291 00292 //compute extrinsic_data 00293 extrinsic_data.set_size(N); 00294 for (n=1; n<=N; n++) 00295 { 00296 sum = -INFINITY; 00297 sumbis = -INFINITY; 00298 for (m=0; m<chtrellis.stateNb; m++) //initial state 00299 { 00300 for (k=0; k<2; k++) //input index 00301 { 00302 nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states 00303 C[k] = (k)*apriori_data[n-1]-itpp::sqr(rec_sig[n-1]-chtrellis.output[m+k*chtrellis.stateNb])/(2*sigma2);//compute log of gamma 00304 buffer = A[m+(n-1)*chtrellis.stateNb]+C[k]+B[nstates[k]+n*chtrellis.stateNb]; 00305 if (k) 00306 sum = std::max(sum, buffer);//1 00307 else 00308 sumbis = std::max(sumbis, buffer);//0 00309 } 00310 } 00311 extrinsic_data[n-1] = (sum-sumbis)-apriori_data[n-1]; 00312 } 00313 00314 //free memory 00315 delete[] chtrellis.output; 00316 delete[] chtrellis.nextState; 00317 delete[] chtrellis.prevState; 00318 delete[] chtrellis.input; 00319 delete[] A; 00320 delete[] B; 00321 } 00322 }//end namespace tr
Generated on Sat Jul 9 2011 15:21:32 for IT++ by Doxygen 1.7.4