1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
/*!
* \file
* \brief EXIT class test program
* \author Bogdan Cristea
*
* -------------------------------------------------------------------------
*
* Copyright (C) 1995-2010 (see AUTHORS file for a list of contributors)
*
* This file is part of IT++ - a C++ library of mathematical, signal
* processing, speech processing, and communications classes and functions.
*
* IT++ is free software: you can redistribute it and/or modify it under the
* terms of the GNU General Public License as published by the Free Software
* Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* IT++ is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
* details.
*
* You should have received a copy of the GNU General Public License along
* with IT++. If not, see <http://www.gnu.org/licenses/>.
*
* -------------------------------------------------------------------------
*/
#include "itpp/itcomm.h"
#include <iomanip>
using namespace itpp;
using std::cout;
using std::endl;
using std::string;
using std::fixed;
int main(void)
{
//general parameters
vec sigmaA = "0.01 2 4 7";//standard deviation (sqrt(variance)) of the mutual a priori information
double threshold_value = 50;
string map_metric="maxlogMAP";
ivec gen = "07 05";//octal form
int constraint_length = 3;
int nb_blocks_lim = 10;
int perm_len = int(itpp::pow10(3.0));//total number of bits in a block (with tail)
double EbN0_dB = 0.8;
double R = 1.0/3.0;//coding rate of PCCC
double Ec = 1.0;//coded bit energy
//other parameters
vec sigma2A = sqr(sigmaA);
int sigma2A_len = sigma2A.length();
int nb_bits = perm_len-(constraint_length-1);//number of bits in a block (without tail)
double sigma2 = (0.5*Ec/R)*pow(inv_dB(EbN0_dB), -1.0);//N0/2
double Lc = -2/sigma2;//normalisation factor for intrinsic information (take into account the BPSK modulation)
bvec bits(nb_bits);
bvec tail;
bvec bits_tail(perm_len);
bmat parity_bits;
int coded_bits_len = 2*perm_len;
bvec coded_bits(coded_bits_len);
vec mod_bits(coded_bits_len);
vec rec_sig(coded_bits_len);
vec intrinsic_coded(coded_bits_len);
vec intrinsic_coded_p(2*nb_bits);
intrinsic_coded_p.zeros();
vec apriori_data(perm_len);
vec extrinsic_coded;
vec extrinsic_data;
vec apriori_mutual_info(sigma2A_len);
vec extrinsic_mutual_info(sigma2A_len);
vec extrinsic_mutual_info_p(sigma2A_len);
extrinsic_mutual_info.zeros();
extrinsic_mutual_info_p.zeros();
register int en,n,nb_blocks;
//Recursive Systematic Convolutional Code
Rec_Syst_Conv_Code rsc;
rsc.set_generator_polynomials(gen, constraint_length);//initial state should be the zero state
//BPSK modulator
BPSK bpsk;
//AWGN channel
AWGN_Channel channel;
channel.set_noise(sigma2);
//SISO module
SISO siso;
siso.set_generators(gen, constraint_length);
siso.set_map_metric(map_metric);
//EXIT chart
EXIT exit;
//Randomize generators
//RNG_randomize();
RNG_reset(12345);
cout << "=============================================" << endl;
cout << " Starting Simulation " << endl;
cout << " Extrinsic Mutual Information (IE) as a function of A priori Mutual Information (IA) " << endl;
cout << "=============================================" << endl;
cout << " Block length = " << perm_len << endl;
cout << " Generator polynomials = " << std::oct << gen << std::dec << endl;
cout << " Turbo encoder rate 1/3 (plus tail bits)" << endl;
cout << "=============================================" << endl;
//main loop
for (en=0; en<sigma2A_len; en++)
{
apriori_mutual_info(en) = exit.apriori_mutual_info(sigma2A(en));//a priori mutual info
for (nb_blocks=0; nb_blocks<nb_blocks_lim; nb_blocks++)
{
//bits generation
bits = randb(nb_bits);
//RSC code
rsc.encode_tail(bits, tail, parity_bits);//tail is added
//form coder output
bits_tail = concat(bits, tail);
for (n=0; n<perm_len; n++)
{
coded_bits(2*n) = bits_tail(n);//systematic output
coded_bits(2*n+1) = parity_bits(n,0);//parity output
}
//BPSK modulation (1->-1,0->+1)
mod_bits = bpsk.modulate_bits(coded_bits);
//AWGN channel
rec_sig = channel(mod_bits);
//first SISO RSC module (tail is added)
//intrinsic info. of coded bits
intrinsic_coded = Lc*rec_sig;
//a priori info. of data bits
apriori_data = exit.generate_apriori_info(bits_tail);
//SISO RSC module
siso.rsc(extrinsic_coded, extrinsic_data, intrinsic_coded, apriori_data, true);
//threshold
extrinsic_data = SISO::threshold(extrinsic_data, threshold_value);
//extrinsic mutual info
extrinsic_mutual_info(en) += exit.extrinsic_mutual_info(extrinsic_data.left(nb_bits), bits);
//second SISO RSC module (no tail added)
//intrinsic info. of coded bits
for (n=0; n<nb_bits; n++)
intrinsic_coded_p(2*n+1) = Lc*rec_sig(2*n+1);//parity bits only
//a priori info. of data bits
apriori_data = exit.generate_apriori_info(bits);
//SISO RSC module
siso.rsc(extrinsic_coded, extrinsic_data, intrinsic_coded_p, apriori_data, false);
//threshold
extrinsic_data = SISO::threshold(extrinsic_data, threshold_value);
//extrinsic mutual info
extrinsic_mutual_info_p(en) += exit.extrinsic_mutual_info(extrinsic_data, bits);
}//end blocks (while loop)
//mean extrinsic mutual info over all blocks
extrinsic_mutual_info(en) /= nb_blocks_lim;
extrinsic_mutual_info_p(en) /= nb_blocks_lim;
}
cout << "IA = " << fixed << round_to_infty(apriori_mutual_info) << endl;
cout << "IE = " << round_to_infty(extrinsic_mutual_info) << endl;
cout << "IE^p = " << round_to_infty(extrinsic_mutual_info_p) << endl;
return 0;
}
|
