Simulation of PCCCs in an AWGN channel

This program simulates Parallel Concatenated Convolutional Codes (PCCCs) of coding rate 1/3 using a turbo decoder with two SISO RSC modules.

Reference: S. Benedetto, D. Divsalar, G. Motorsi and F. Pollara, "A Soft-Input Soft-Output Maximum A posteriori (MAP) Module to Decode Parallel and Serial Concatenated Codes", TDA Progress Report, nov. 1996

#include "itpp/itcomm.h"
using namespace itpp;
using std::cout;
using std::endl;
using std::string;
int main(void)
{
//general parameters
double threshold_value = 10;
string map_metric="maxlogMAP";
ivec gen = "07 05";//octal form, feedback first
int constraint_length = 3;
int nb_errors_lim = 3000;
int nb_bits_lim = int(1e6);
int perm_len = (1<<14);//total number of bits in a block (with tail)
int nb_iter = 10;//number of iterations in the turbo decoder
vec EbN0_dB = "0:0.1:5";
double R = 1.0/3.0;//coding rate (non punctured PCCC)
double Ec = 1.0;//coded bit energy
//other parameters
int nb_bits = perm_len-(constraint_length-1);//number of bits in a block (without tail)
vec sigma2 = (0.5*Ec/R)*pow(inv_dB(EbN0_dB), -1.0);//N0/2
double Lc;//scaling factor
int nb_blocks;//number of blocks
int nb_errors;
ivec perm(perm_len);
ivec inv_perm(perm_len);
bvec bits(nb_bits);
int cod_bits_len = perm_len*gen.length();
bmat cod1_bits;//tail is added
bvec tail;
bvec cod2_input;
bmat cod2_bits;
int rec_len = int(1.0/R)*perm_len;
bvec coded_bits(rec_len);
vec rec(rec_len);
vec dec1_intrinsic_coded(cod_bits_len);
vec dec2_intrinsic_coded(cod_bits_len);
vec apriori_data(perm_len);//a priori LLR for information bits
vec extrinsic_coded(perm_len);
vec extrinsic_data(perm_len);
bvec rec_bits(perm_len);
int snr_len = EbN0_dB.length();
mat ber(nb_iter,snr_len);
ber.zeros();
register int en,n;
//Recursive Systematic Convolutional Code
Rec_Syst_Conv_Code cc;
cc.set_generator_polynomials(gen, constraint_length);//initial state should be the zero state
//BPSK modulator
BPSK bpsk;
//AWGN channel
AWGN_Channel channel;
//SISO modules
SISO siso;
siso.set_generators(gen, constraint_length);
siso.set_map_metric(map_metric);
//BER
BERC berc;
//Randomize generators
RNG_randomize();
//main loop
for (en=0;en<snr_len;en++)
{
cout << "EbN0_dB = " << EbN0_dB(en) << endl;
channel.set_noise(sigma2(en));
Lc = -2/sigma2(en);//normalisation factor for intrinsic information (take into account the BPSK mapping)
nb_errors = 0;
nb_blocks = 0;
while ((nb_errors<nb_errors_lim) && (nb_blocks*nb_bits<nb_bits_lim))
{
//permutation
perm = sort_index(randu(perm_len));
//inverse permutation
inv_perm = sort_index(perm);
//bits generation
bits = randb(nb_bits);
//parallel concatenated convolutional code
cc.encode_tail(bits, tail, cod1_bits);//tail is added here to information bits to close the trellis
cod2_input = concat(bits, tail);
cc.encode(cod2_input(perm), cod2_bits);
for (n=0;n<perm_len;n++)//output with no puncturing
{
coded_bits(3*n) = cod2_input(n);//systematic output
coded_bits(3*n+1) = cod1_bits(n,0);//first parity output
coded_bits(3*n+2) = cod2_bits(n,0);//second parity output
}
//BPSK modulation (1->-1,0->+1) + AWGN channel
rec = channel(bpsk.modulate_bits(coded_bits));
//form input for SISO blocks
for (n=0;n<perm_len;n++)
{
dec1_intrinsic_coded(2*n) = Lc*rec(3*n);
dec1_intrinsic_coded(2*n+1) = Lc*rec(3*n+1);
dec2_intrinsic_coded(2*n) = 0.0;//systematic output of the CC is already used in decoder1
dec2_intrinsic_coded(2*n+1) = Lc*rec(3*n+2);
}
//turbo decoder
apriori_data.zeros();//a priori LLR for information bits
for (n=0;n<nb_iter;n++)
{
//first decoder (terminated trellis)
siso.rsc(extrinsic_coded, extrinsic_data, dec1_intrinsic_coded, apriori_data, true);
//interleave
apriori_data = extrinsic_data(perm);
//threshold
apriori_data = SISO::threshold(apriori_data, threshold_value);
//second decoder (unterminated trellis)
siso.rsc(extrinsic_coded, extrinsic_data, dec2_intrinsic_coded, apriori_data);
//decision
apriori_data += extrinsic_data;//a posteriori information
rec_bits = bpsk.demodulate_bits(-apriori_data(inv_perm));//take into account the BPSK mapping
//count errors
berc.clear();
berc.count(bits, rec_bits.left(nb_bits));
ber(n,en) += berc.get_errorrate();
//deinterleave for the next iteration
apriori_data = extrinsic_data(inv_perm);
}//end iterations
nb_errors += int(berc.get_errors());//get number of errors at the last iteration
nb_blocks++;
}//end blocks (while loop)
//compute BER over all tx blocks
ber.set_col(en, ber.get_col(en)/nb_blocks);
}
it_file ff("pccc_bersim_awgn.it");
ff << Name("EbN0_dB") << EbN0_dB;
ff << Name("BER") << ber;
ff.close();
return 0;
}
itpp::SISO::threshold
static double threshold(const double &x, const double &value)
Functions used to limit values at a given +- threshold.
Definition: siso.h:901
itpp::pow
vec pow(const double x, const vec &y)
Calculates x to the power of y (x^y)
Definition: log_exp.h:176
itpp::inv_dB
double inv_dB(double x)
Inverse of decibel of x.
Definition: log_exp.h:73
itpp::RNG_randomize
void RNG_randomize()
Set a random seed for all Random Number Generators in the current thread.
Definition: random.cpp:254
itpp::randu
double randu(void)
Generates a random uniform (0,1) number.
Definition: random.h:804
itpp::randb
bin randb(void)
Generates a random bit (equally likely 0s and 1s)
Definition: random.h:793
itcomm.h
Include file for the IT++ communications module.
bmat
Mat< bin > bmat
bin matrix
Definition: mat.h:508
itpp
itpp namespace
Definition: itmex.h:37
itpp::concat
const Array< T > concat(const Array< T > &a, const T &e)
Append element e to the end of the Array a.
Definition: array.h:486

When you run this program, the results (BER and EbN0_dB) are saved into pccc_bersim_awgn.it file. Using the following MATLAB script:

clear all
itload('pccc_bersim_awgn.it');
figure
semilogy(EbN0_dB, BER, 'o-')
grid on
xlabel('E_b/N_0 [dB]')
ylabel('BER')
itpp::dB
double dB(double x)
Decibel of x (10*log10(x))
Definition: log_exp.h:71
itpp::all
ITPP_EXPORT bool all(const bvec &testvec)
Returns true if all elements are ones and false otherwise.

the results can be displayed.

Similarly, the results can be displayed using the following Python script (pyitpp, numpy and matplotlib modules are required):

#!/usr/bin/env python
from pyitpp import itload
from matplotlib.pyplot import *
out = itload('pccc_bersim_awgn.it')
semilogy(out['EbN0_dB'], out['BER'].T, 'o-')
grid()
xlabel('$E_b/N_0 [dB]$')
ylabel('$BER$')
show()

Generated on Tue Aug 17 2021 10:59:15 for IT++ by Doxygen 1.9.1