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
35namespace itpp
36{
37void 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
104void 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
219void 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
double log_add(double log_a, double log_b)
Safe substitute for log(exp(log_a) + exp(log_b))
Definition: log_exp.cpp:40
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition: matfunc.h:59
vec sqr(const cvec &data)
Absolute square of elements.
Definition: elem_math.cpp:36
itpp namespace
Definition: itmex.h:37
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_EXPORT bvec dec2bin(int length, int index)
Convert a decimal int index to bvec using length bits in the representation.
ivec to_ivec(const Vec< T > &v)
Converts a Vec<T> to ivec.
Definition: converters.h:79
bvec to_bvec(const Vec< T > &v)
Converts a Vec<T> to bvec.
Definition: converters.h:51
Definitions for Soft Input Soft Output (SISO) modules class.
SourceForge Logo

Generated on Tue Jan 24 2023 00:00:00 for IT++ by Doxygen 1.9.6