IT++ Logo
siso_eq.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::gen_chtrellis(void)
38 // generate trellis for precoded FIR channels with real coefficients
39 {
40  //get parameters
41  int mem_len = impulse_response.cols()-1;//memory length of the channel
42  int p_order = prec_gen.length()-1;//order of the precoder polynomial
43 
44  //other variables
45  register int n,k,j;
46  double inputs[] = {1.0,-1.0};//1->-1, 0->+1
47  int index;
48  double feedback[2];
49 
50  //create channel trellis
51  int equiv_ch_mem_len = std::max(mem_len, p_order);
52  chtrellis.stateNb = (1<<equiv_ch_mem_len);
53  chtrellis.output = new double[2*chtrellis.stateNb];
54  chtrellis.nextState = new int[2*chtrellis.stateNb];
55  chtrellis.prevState = new int[2*chtrellis.stateNb];
56  chtrellis.input = new int[2*chtrellis.stateNb];
57 
58  //initialize trellis
59  itpp::ivec enc_mem(equiv_ch_mem_len);
60 #pragma omp parallel for private(n,enc_mem,k,feedback,j)
61  for (n=0; n<chtrellis.stateNb; n++) //initial state
62  {
63  enc_mem = 1-2*itpp::to_ivec(itpp::dec2bin(equiv_ch_mem_len, n));//1->-1, 0->+1
64  //output
65  for (k=0; k<2; k++)
66  {
67  feedback[k] = inputs[k];
68  for (j=1; j<=p_order; j++)
69  if (prec_gen[j])
70  feedback[k] *= enc_mem[j-1];//xor truth table must remain the same
71  chtrellis.output[n+k*chtrellis.stateNb] = feedback[k]*impulse_response(0,0);
72  for (j=1; j<=mem_len; j++)
73  chtrellis.output[n+k*chtrellis.stateNb] += (enc_mem[j-1]*impulse_response(0,j));
74  }
75  //next state
76  for (j=(equiv_ch_mem_len-1); j>0; j--)
77  enc_mem[j] = enc_mem[j-1];
78  for (k=0; k<2; k++)
79  {
80  enc_mem[0] = int(feedback[k]);
81  chtrellis.nextState[n+k*chtrellis.stateNb] = itpp::bin2dec(itpp::to_bvec((1-enc_mem)/2), true);//-1->1, +1->0
82  }
83  }
84 
85 #pragma omp parallel for private(j,index,n,k)
86  for (j=0; j<chtrellis.stateNb; j++)
87  {
88  index = 0;
89  for (n=0; n<chtrellis.stateNb; n++)
90  {
91  for (k=0; k<2; k++)
92  {
93  if (chtrellis.nextState[n+k*chtrellis.stateNb]==j)
94  {
95  chtrellis.prevState[j+index*chtrellis.stateNb] = n;
96  chtrellis.input[j+index*chtrellis.stateNb] = k;//this is an index to the channel input
97  index++;
98  }
99  }
100  }
101  }
102 }
103 
104 void SISO::equalizer_logMAP(itpp::vec &extrinsic_data, const itpp::vec &rec_sig, const itpp::vec &apriori_data)
105 /*
106  extrinsic_data - extrinsic information of channel input symbols
107  rec_sig - received symbols
108  apriori_data - a priori information of channel input symbols
109  */
110 {
111  //get parameters
112  int N = rec_sig.length();//length of the received frame
113  //other parameters
114  register int n,k,m;
115  int pstates[2];
116  int nstates[2];
117  int inputs[2];
118  double C[2];//log(gamma)
119  double sum;
120  double sumbis;
121  double buffer;
122 
123  //initialize trellis
124  gen_chtrellis();
125  //initialize log(alpha) and log(beta)
126  double* A = new double[chtrellis.stateNb*(N+1)];
127  double* B = new double[chtrellis.stateNb*(N+1)];
128  A[0] = 0;
129  B[N*chtrellis.stateNb] = 0;
130  sum = (tail?-INFINITY:0);
131 #pragma omp parallel for private(n)
132  for (n=1; n<chtrellis.stateNb; n++)
133  {
134  A[n] = -INFINITY;
135  B[n+N*chtrellis.stateNb] = sum;//if tail==false the final state is not known
136  }
137 
138 #pragma omp parallel sections private(n,sum,m,k,C)
139  {
140  //forward recursion
141  for (n=1; n<=N; n++)
142  {
143  sum = 0;//normalisation factor
144  for (m=0; m<chtrellis.stateNb; m++) //final state
145  {
146  for (k=0; k<2; k++)
147  {
148  pstates[k] = chtrellis.prevState[m+chtrellis.stateNb*k];//determine previous states
149  inputs[k] = chtrellis.input[m+chtrellis.stateNb*k];//determine input
150  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
151  }
152  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]);
153  sum += std::exp(A[m+n*chtrellis.stateNb]);
154  }
155  //normalisation
156  sum = std::log(sum);
157  for (m=0; m<chtrellis.stateNb; m++)
158  {
159  A[m+n*chtrellis.stateNb] -= sum;
160  }
161  }
162 
163  //backward recursion
164 #pragma omp section
165  for (n=N-1; n>=0; n--)
166  {
167  sum = 0;//normalisation factor
168  for (m=0; m<chtrellis.stateNb; m++) //initial state
169  {
170  for (k=0; k<2; k++)
171  {
172  nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states
173  C[k] = (k)*apriori_data[n]-itpp::sqr(rec_sig[n]-chtrellis.output[m+k*chtrellis.stateNb])/(2*sigma2);//compute log of gamma
174  }
175  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]);
176  sum += std::exp(B[m+n*chtrellis.stateNb]);
177  }
178  //normalisation
179  sum = std::log(sum);
180  for (m=0; m<chtrellis.stateNb; m++)
181  {
182  B[m+n*chtrellis.stateNb] -= sum;
183  }
184  }
185  }
186 
187  //compute extrinsic_data
188  extrinsic_data.set_size(N);
189 #pragma omp parallel for private(n,sum,sumbis,m,k,nstates,C,buffer)
190  for (n=1; n<=N; n++)
191  {
192  sum = 0;//could be replaced by a vector
193  sumbis = 0;
194  for (m=0; m<chtrellis.stateNb; m++) //initial state
195  {
196  for (k=0; k<2; k++) //input index
197  {
198  nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states
199  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
200  buffer = std::exp(A[m+(n-1)*chtrellis.stateNb]+C[k]+B[nstates[k]+n*chtrellis.stateNb]);
201  if (k)
202  sum += buffer;//1
203  else
204  sumbis += buffer;//0
205  }
206  }
207  extrinsic_data[n-1] = std::log(sum/sumbis)-apriori_data[n-1];
208  }
209 
210  //free memory
211  delete[] chtrellis.output;
212  delete[] chtrellis.nextState;
213  delete[] chtrellis.prevState;
214  delete[] chtrellis.input;
215  delete[] A;
216  delete[] B;
217 }
218 
219 void SISO::equalizer_maxlogMAP(itpp::vec &extrinsic_data, const itpp::vec &rec_sig, const itpp::vec &apriori_data)
220 /*
221  extrinsic_data - extrinsic information for channel input symbols
222  rec_sig - received symbols
223  apriori_data - a priori information for channel input symbols
224  */
225 {
226  //get parameters
227  int N = rec_sig.length();//length of the received frame
228  //other parameters
229  register int n,k,m;
230  int pstates[2];
231  int nstates[2];
232  int inputs[2];
233  double C[2];//log(gamma)
234  double sum;
235  double sumbis;
236  double buffer;
237 
238  //initialize trellis
239  gen_chtrellis();
240  //initialize log(alpha) and log(beta)
241  double* A = new double[chtrellis.stateNb*(N+1)];
242  double* B = new double[chtrellis.stateNb*(N+1)];
243  A[0] = 0;
244  for (n=1; n<chtrellis.stateNb; n++)
245  A[n] = -INFINITY;
246  B[N*chtrellis.stateNb] = 0;
247  sum = (tail?-INFINITY:0);
248  for (n=1; n<chtrellis.stateNb; n++)
249  B[n+N*chtrellis.stateNb] = sum;//if tail==false the final state is not known
250 
251 #pragma omp parallel sections private(n,sum,m,k,C)
252  {
253  //forward recursion
254  for (n=1; n<=N; n++)
255  {
256  sum = -INFINITY;//normalization factor
257  for (m=0; m<chtrellis.stateNb; m++) //final state
258  {
259  for (k=0; k<2; k++)
260  {
261  pstates[k] = chtrellis.prevState[m+chtrellis.stateNb*k];//determine previous states
262  inputs[k] = chtrellis.input[m+chtrellis.stateNb*k];//determine input
263  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
264  }
265  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]);
266  sum = std::max(sum, A[m+n*chtrellis.stateNb]);
267  }
268  for (m=0; m<chtrellis.stateNb; m++)
269  A[m+n*chtrellis.stateNb] -= sum;
270  }
271 
272  //backward recursion
273 #pragma omp section
274  for (n=N-1; n>=0; n--)
275  {
276  sum = -INFINITY;//normalisation factor
277  for (m=0; m<chtrellis.stateNb; m++) //initial state
278  {
279  for (k=0; k<2; k++)
280  {
281  nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states
282  C[k] = (k)*apriori_data[n]-itpp::sqr(rec_sig[n]-chtrellis.output[m+k*chtrellis.stateNb])/(2*sigma2);//compute log of gamma
283  }
284  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]);
285  sum = std::max(sum, B[m+n*chtrellis.stateNb]);
286  }
287  for (m=0; m<chtrellis.stateNb; m++)
288  B[m+n*chtrellis.stateNb] -= sum;
289  }
290  }
291 
292  //compute extrinsic_data
293  extrinsic_data.set_size(N);
294  for (n=1; n<=N; n++)
295  {
296  sum = -INFINITY;
297  sumbis = -INFINITY;
298  for (m=0; m<chtrellis.stateNb; m++) //initial state
299  {
300  for (k=0; k<2; k++) //input index
301  {
302  nstates[k] = chtrellis.nextState[m+k*chtrellis.stateNb];//determine next states
303  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
304  buffer = A[m+(n-1)*chtrellis.stateNb]+C[k]+B[nstates[k]+n*chtrellis.stateNb];
305  if (k)
306  sum = std::max(sum, buffer);//1
307  else
308  sumbis = std::max(sumbis, buffer);//0
309  }
310  }
311  extrinsic_data[n-1] = (sum-sumbis)-apriori_data[n-1];
312  }
313 
314  //free memory
315  delete[] chtrellis.output;
316  delete[] chtrellis.nextState;
317  delete[] chtrellis.prevState;
318  delete[] chtrellis.input;
319  delete[] A;
320  delete[] B;
321 }
322 }//end namespace tr
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:23 for IT++ by Doxygen 1.8.2