siso_rsc.cpp
Go to the documentation of this file.
#include <itpp/comm/siso.h>
#include <itpp/base/itcompat.h>
#include <limits>
#ifndef INFINITY
#define INFINITY std::numeric_limits<double>::infinity()
#endif
namespace itpp
{
void SISO::gen_rsctrellis(void)
//generates 1/2 RSC trellis structure for binary symbols
//the states are numbered from 0
{
int mem_len = gen.cols()-1;
register int n,k,j;
itpp::bin feedback,out;
int buffer;
rsctrellis.numStates = (1<<mem_len);
rsctrellis.prevStates = new int[2*rsctrellis.numStates];
rsctrellis.nextStates = new int[2*rsctrellis.numStates];
rsctrellis.PARout = new double[2*rsctrellis.numStates];
rsctrellis.fm = new itpp::bin[rsctrellis.numStates];
itpp::bvec cases(mem_len);
for (n=0; n<2; n++)
{
#pragma omp parallel for private(k, j, cases, feedback, out, buffer)
for (k=0; k<rsctrellis.numStates; k++)
{
cases = itpp::dec2bin(mem_len, k);
//feedback
feedback = (itpp::bin)n;
for (j=1; j<(mem_len+1); j++)
{
feedback ^= (gen(0,j)*cases[j-1]);
}
//out
out = feedback*gen(1,0);
for (j=1; j<(mem_len+1); j++)
{
out ^= (gen(1,j)*cases[j-1]);
}
rsctrellis.PARout[k+n*rsctrellis.numStates] = (out?1.0:0.0);//parity bit
rsctrellis.fm[k] = itpp::bin(n)^out;
//shift
for (j=mem_len-1; j>0; j--)
{
cases[j] = cases[j-1];
}
cases[0] = feedback;
//next and previous state
buffer = itpp::bin2dec(cases, true);
rsctrellis.nextStates[k+n*rsctrellis.numStates] = buffer;//next state
rsctrellis.prevStates[buffer+n*rsctrellis.numStates] = k;//previous state
}
}
}
void SISO::rsc_logMAP(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data,
const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data)
/*
logMAP (SISO) decoder for RSC of rate 1/2
extrinsic_coded - extrinsic information of coded bits
extrinsic_data - extrinsic information of data (informational) bits
intrinsic_coded - intrinsic information of coded (systematic and parity) bits
apriori_data - a priori information of data (informational) bits
Reference: Steven S. Pietrobon and Adrian S. Barbulescu, "A simplification of
the modified Bahl decoding algorithm for systematic convolutional codes", Proc. ISITA, 1994
*/
{
//get parameters
int N = apriori_data.length();
//other parameters
register int n,k;
double buffer;
int index;
double A_min, A_max;
double sum0, sum1;
//trellis generation
gen_rsctrellis();
//parameter initialization
double* Lc1I = new double[N];
double* Lc2I = new double[N];
#pragma omp parallel for private(n)
for (n=0; n<N; n++)
{
Lc1I[n] = intrinsic_coded[2*n];
Lc2I[n] = intrinsic_coded[2*n+1];
}
double* A0 = new double[rsctrellis.numStates*N];
double* A1 = new double[rsctrellis.numStates*N];
double* A_mid = new double[N];
double* B0 = new double[rsctrellis.numStates*N];
double* B1 = new double[rsctrellis.numStates*N];
buffer = (tail?-INFINITY:0);//log(buffer)
#pragma omp parallel for private(n,k)
for (n=0; n<N; n++)
{
for (k=0; k<rsctrellis.numStates; k++)
{
A0[k+n*rsctrellis.numStates] = -INFINITY;
A1[k+n*rsctrellis.numStates] = -INFINITY;
B0[k+n*rsctrellis.numStates] = buffer;
B1[k+n*rsctrellis.numStates] = buffer;
}
A_mid[n] = 0;
}
//A
A0[0] = Lc2I[0]*rsctrellis.PARout[0];//i=0
A1[0] = Lc1I[0]+apriori_data[0]+Lc2I[0]*rsctrellis.PARout[rsctrellis.numStates];//i=1
for (n=1; n<N; n++)
{
A_min = INFINITY;
A_max = 0;
for (k=0; k<rsctrellis.numStates; k++)
{
buffer = itpp::log_add(A0[rsctrellis.prevStates[k]+(n-1)*rsctrellis.numStates],
A1[rsctrellis.prevStates[k+rsctrellis.numStates]+(n-1)*rsctrellis.numStates]);//log(alpha0+alpha1)
A0[k+rsctrellis.numStates*n] = Lc2I[n]*rsctrellis.PARout[k]+buffer;
A1[k+rsctrellis.numStates*n] = Lc1I[n]+apriori_data[n]+
Lc2I[n]*rsctrellis.PARout[k+rsctrellis.numStates]+buffer;
//find min A(:,n)
A_min = std::min(A_min, A0[k+rsctrellis.numStates*n]);
//find max A(:,n)
A_max = std::max(A_max, A0[k+rsctrellis.numStates*n]);
}
//normalization
A_mid[n] = (A_min+A_max)/2;
if (std::isinf(A_mid[n]))
{
continue;
}
for (k=0; k<rsctrellis.numStates; k++)
{
A0[k+rsctrellis.numStates*n] -= A_mid[n];
A1[k+rsctrellis.numStates*n] -= A_mid[n];
}
}
//B
B0[rsctrellis.prevStates[0]+(N-1)*rsctrellis.numStates] = 0;
B1[rsctrellis.prevStates[rsctrellis.numStates]+(N-1)*rsctrellis.numStates] = 0;
for (n=N-2; n>=0; n--)
{
for (k=0; k<rsctrellis.numStates; k++)
{
index = rsctrellis.nextStates[k];
B0[k+rsctrellis.numStates*n] = itpp::log_add(B0[index+(n+1)*rsctrellis.numStates]+
Lc2I[n+1]*rsctrellis.PARout[index],
B1[index+(n+1)*rsctrellis.numStates]+
Lc1I[n+1]+apriori_data[n+1]+Lc2I[n+1]*rsctrellis.PARout[index+rsctrellis.numStates]);
index = rsctrellis.nextStates[k+rsctrellis.numStates];
B1[k+rsctrellis.numStates*n] = itpp::log_add(B0[index+(n+1)*rsctrellis.numStates]+
Lc2I[n+1]*rsctrellis.PARout[index],
B1[index+(n+1)*rsctrellis.numStates]+
Lc1I[n+1]+apriori_data[n+1]+Lc2I[n+1]*rsctrellis.PARout[index+rsctrellis.numStates]);
}
if (std::isinf(A_mid[n+1]))
{
continue;
}
for (k=0; k<rsctrellis.numStates; k++)
{
B0[k+rsctrellis.numStates*n] -= A_mid[n+1];
B1[k+rsctrellis.numStates*n] -= A_mid[n+1];
}
}
//updated LLR for information bits
extrinsic_data.set_size(N);
extrinsic_coded.set_size(2*N);
#pragma omp parallel for private(n, k, sum0, sum1)
for (n=0; n<N; n++)
{
sum0 = 0;
sum1 = 0;
for (k=0; k<rsctrellis.numStates; k++)
{
sum1 += std::exp(A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]);
sum0 += std::exp(A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]);
}
extrinsic_data[n] = std::log(sum1/sum0)-apriori_data[n];//updated information must be independent of input LLR
extrinsic_coded[2*n] = std::log(sum1/sum0)-Lc1I[n];//this information is used in serial concatenations
}
//updated LLR for coded (parity) bits
#pragma omp parallel for private(n, k, sum0, sum1)
for (n=0; n<N; n++)
{
sum0 = 0;
sum1 = 0;
for (k=0; k<rsctrellis.numStates; k++)
{
if (rsctrellis.fm[k])
{
sum1 += std::exp(A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]);
sum0 += std::exp(A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]);
}
else
{
sum0 += std::exp(A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]);
sum1 += std::exp(A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]);
}
}
extrinsic_coded[2*n+1] = std::log(sum0/sum1)-Lc2I[n];//updated information must be independent of input LLR
}
//destroy trellis
delete[] rsctrellis.prevStates;
delete[] rsctrellis.nextStates;
delete[] rsctrellis.PARout;
delete[] rsctrellis.fm;
//destroy MAP parameters
delete[] Lc1I;
delete[] Lc2I;
delete[] A0;
delete[] A1;
delete[] A_mid;
delete[] B0;
delete[] B1;
}
void SISO::rsc_maxlogMAP(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data,
const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data)
/*
maxlogMAP (SISO) decoder for RSC of rate 1/2
extrinsic_coded - extrinsic information of coded bits
extrinsic_data - extrinsic information of data (informational) bits
intrinsic_coded - intrinsic information of coded (systematic and parity) bits
apriori_data - a priori information of data (informational) bits
Reference: Steven S. Pietrobon and Adrian S. Barbulescu, "A simplification of
the modified Bahl decoding algorithm for systematic convolutional codes", Proc. ISITA, 1994
*/
{
//get parameters
int N = apriori_data.length();
//other parameters
register int n,k;
double buffer;
int index;
double A_min, A_max;
double sum0, sum1;
//trellis generation
gen_rsctrellis();
//parameter initialization
double* Lc1I = new double[N];
double* Lc2I = new double[N];
#pragma omp parallel for private(n)
for (n=0; n<N; n++)
{
Lc1I[n] = intrinsic_coded[2*n];
Lc2I[n] = intrinsic_coded[2*n+1];
}
double* A0 = new double[rsctrellis.numStates*N];
double* A1 = new double[rsctrellis.numStates*N];
double* A_mid = new double[N];
double* B0 = new double[rsctrellis.numStates*N];
double* B1 = new double[rsctrellis.numStates*N];
buffer = (tail?-INFINITY:0);//log(buffer)
#pragma omp parallel for private(n,k)
for (n=0; n<N; n++)
{
for (k=0; k<rsctrellis.numStates; k++)
{
A0[k+n*rsctrellis.numStates] = -INFINITY;
A1[k+n*rsctrellis.numStates] = -INFINITY;
B0[k+n*rsctrellis.numStates] = buffer;
B1[k+n*rsctrellis.numStates] = buffer;
}
A_mid[n] = 0;
}
//A
A0[0] = Lc2I[0]*rsctrellis.PARout[0];//i=0
A1[0] = Lc1I[0]+apriori_data[0]+Lc2I[0]*rsctrellis.PARout[rsctrellis.numStates];//i=1
for (n=1; n<N; n++)
{
A_min = INFINITY;
A_max = 0;
for (k=0; k<rsctrellis.numStates; k++)
{
buffer = std::max(A0[rsctrellis.prevStates[k]+(n-1)*rsctrellis.numStates],
A1[rsctrellis.prevStates[k+rsctrellis.numStates]+(n-1)*rsctrellis.numStates]);//log(alpha0+alpha1)
A0[k+rsctrellis.numStates*n] = Lc2I[n]*rsctrellis.PARout[k]+buffer;
A1[k+rsctrellis.numStates*n] = Lc1I[n]+apriori_data[n]+
Lc2I[n]*rsctrellis.PARout[k+rsctrellis.numStates]+buffer;
//find min A(:,n)
A_min = std::min(A_min, A0[k+rsctrellis.numStates*n]);
//find max A(:,n)
A_max = std::max(A_max, A0[k+rsctrellis.numStates*n]);
}
//normalization
A_mid[n] = (A_min+A_max)/2;
if (std::isinf(A_mid[n]))
continue;
for (k=0; k<rsctrellis.numStates; k++)
{
A0[k+rsctrellis.numStates*n] -= A_mid[n];
A1[k+rsctrellis.numStates*n] -= A_mid[n];
}
}
//B
B0[rsctrellis.prevStates[0]+(N-1)*rsctrellis.numStates] = 0;
B1[rsctrellis.prevStates[rsctrellis.numStates]+(N-1)*rsctrellis.numStates] = 0;
for (n=N-2; n>=0; n--)
{
for (k=0; k<rsctrellis.numStates; k++)
{
index = rsctrellis.nextStates[k];
B0[k+rsctrellis.numStates*n] = std::max(B0[index+(n+1)*rsctrellis.numStates]+
Lc2I[n+1]*rsctrellis.PARout[index],
B1[index+(n+1)*rsctrellis.numStates]+
Lc1I[n+1]+apriori_data[n+1]+Lc2I[n+1]*rsctrellis.PARout[index+rsctrellis.numStates]);
index = rsctrellis.nextStates[k+rsctrellis.numStates];
B1[k+rsctrellis.numStates*n] = std::max(B0[index+(n+1)*rsctrellis.numStates]+
Lc2I[n+1]*rsctrellis.PARout[index],
B1[index+(n+1)*rsctrellis.numStates]+
Lc1I[n+1]+apriori_data[n+1]+Lc2I[n+1]*rsctrellis.PARout[index+rsctrellis.numStates]);
}
if (std::isinf(A_mid[n+1]))
continue;
for (k=0; k<rsctrellis.numStates; k++)
{
B0[k+rsctrellis.numStates*n] -= A_mid[n+1];
B1[k+rsctrellis.numStates*n] -= A_mid[n+1];
}
}
//updated LLR for information bits
extrinsic_data.set_size(N);
extrinsic_coded.set_size(2*N);
#pragma omp parallel for private(n, k, sum0, sum1)
for (n=0; n<N; n++)
{
sum0 = -INFINITY;
sum1 = -INFINITY;
for (k=0; k<rsctrellis.numStates; k++)
{
sum1 = std::max(sum1, A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]);
sum0 = std::max(sum0, A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]);
}
extrinsic_data[n] = (sum1-sum0)-apriori_data[n];//updated information must be independent of input LLR
extrinsic_coded[2*n] = (sum1-sum0)-Lc1I[n];
}
//updated LLR for coded (parity) bits
#pragma omp parallel for private(n, k, sum0, sum1)
for (n=0; n<N; n++)
{
sum0 = -INFINITY;
sum1 = -INFINITY;
for (k=0; k<rsctrellis.numStates; k++)
{
if (rsctrellis.fm[k])
{
sum1 = std::max(sum1, A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]);
sum0 = std::max(sum0, A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]);
}
else
{
sum0 = std::max(sum0, A1[k+n*rsctrellis.numStates]+B1[k+n*rsctrellis.numStates]);
sum1 = std::max(sum1, A0[k+n*rsctrellis.numStates]+B0[k+n*rsctrellis.numStates]);
}
}
extrinsic_coded[2*n+1] = (sum0-sum1)-Lc2I[n];//updated information must be independent of input LLR
}
//destroy trellis
delete[] rsctrellis.prevStates;
delete[] rsctrellis.nextStates;
delete[] rsctrellis.PARout;
delete[] rsctrellis.fm;
//destroy MAP parameters
delete[] Lc1I;
delete[] Lc2I;
delete[] A0;
delete[] A1;
delete[] A_mid;
delete[] B0;
delete[] B1;
}
void SISO::rsc_sova(itpp::vec &extrinsic_data, const itpp::vec &intrinsic_coded,
const itpp::vec &apriori_data, const int &win_len)
/* Soft Output Viterbi Algorithm (SOVA) optimized for 1/N encoders
* Output: extrinsic_data - extrinsic information of data bits
* Inputs: intrinsic_coded - intrinsic information of data bits
* apriori_data - a priori information of data bits
* win_len - window length used to represent the code trellis
*
* The original version has been written by Adrian Bohdanowicz (2003).
* It is assumed that the BPSK mapping is: 0 -> +1, 1 -> -1.
* Changes have been made to adapt the code for RSC codes of rate 1/2
* and for soft input informations.
* Improved SOVA has been implemented using a scaling factor and threshold for
* the reliability information (see Wang [2003]). Even so, PCCC performance
* are close to the original SOVA.
*/
{
//number of code outputs
int nb_outputs = gen.rows();
//setup internal variables based on RSC trellis
register int i,j,s;
gen_rsctrellis();//trellis generation for 1/2 RSC codes
int nb_states = rsctrellis.numStates;
itpp::Array<itpp::mat> bin_out(2);//contains code output for each initial state and code input
itpp::imat next_state(nb_states,2);//next state in the trellis
for (i=0; i<2; i++)
{
bin_out(i).set_size(nb_states, nb_outputs);
for (j=0; j<nb_states; j++)
{
bin_out(i)(j,0) = double(i);//systematic bit
bin_out(i)(j,1) = rsctrellis.PARout[j+i*nb_states];//parity bit
next_state(j,i) = rsctrellis.nextStates[j+i*nb_states];
}
}
itpp::vec bin_inp("0 1");//binary code inputs
int len = apriori_data.length();//number of decoding steps (total)
//allocate memory for the trellis window
itpp::mat metr(nb_states,win_len+1);//path metric buffer
metr.zeros();
metr += -INFINITY;
metr(0,0) = 0;//initial state => (0,0)
itpp::mat surv(nb_states,win_len+1);//survivor state buffer
surv.zeros();
itpp::mat inpt(nb_states,win_len+1);//survivor input buffer (dec. output)
inpt.zeros();
itpp::mat diff(nb_states,win_len+1);//path metric difference
diff.zeros();
itpp::mat comp(nb_states,win_len+1);//competitor state buffer
comp.zeros();
itpp::mat inpc(nb_states,win_len+1);//competitor input buffer
inpc.zeros();
//soft output (sign with reliability)
itpp::vec sft(len);
sft.zeros();
sft += INFINITY;
//decode all the bits
int Cur,Nxt,nxt,sur,b,tmp,idx;
itpp::vec buf(nb_outputs);
double llb,mtr,dif,cmp,inc,srv,inp;
itpp::vec bin(nb_outputs);
itpp::ivec cyclic_buffer(win_len);
extrinsic_data.set_size(len);
for (i = 0; i < len; i++)
{
//indices + precalculations
Cur = i%(win_len+1);//curr trellis (cycl. buf) position
Nxt = (i+1)%(win_len+1);//next trellis (cycl. buf) position
buf = intrinsic_coded(i*nb_outputs,(i+1)*nb_outputs-1);//intrinsic_info portion to be processed
llb = apriori_data(i);//SOVA: apriori_info portion to be processed
metr.set_col(Nxt, -INFINITY*itpp::ones(nb_states));
//forward recursion
for (s = 0; s<nb_states; s++)
{
for (j = 0; j<2; j++)
{
nxt = next_state(s,j);//state after transition
bin = bin_out(j).get_row(s);//transition output (encoder)
mtr = bin*buf+metr(s,Cur);//transition metric
mtr += bin_inp(j)*llb;//SOVA
if (metr(nxt,Nxt) < mtr)
{
diff(nxt,Nxt) = mtr-metr(nxt,Nxt);//SOVA
comp(nxt,Nxt) = surv(nxt,Nxt);//SOVA
inpc(nxt,Nxt) = inpt(nxt,Nxt);//SOVA
metr(nxt,Nxt) = mtr;//store the metric
surv(nxt,Nxt) = s;//store the survival state
inpt(nxt,Nxt) = j;//store the survival input
}
else
{
dif = metr(nxt,Nxt)-mtr;
if (dif <= diff(nxt,Nxt))
{
diff(nxt,Nxt) = dif;//SOVA
comp(nxt,Nxt) = s;//SOVA
inpc(nxt,Nxt) = j;//SOVA
}
}
}
}
//trace backwards
if (i < (win_len-1))
{
continue;
}//proceed if the buffer has been filled
mtr = itpp::max(metr.get_col(Nxt), sur);//find the initial state (max metric)
b = i;//temporary bit index
for (j=0; j<win_len; j++)//indices in a 'cyclic buffer' operation
{
cyclic_buffer(j) = (Nxt-j)%(win_len+1);
cyclic_buffer(j) = (cyclic_buffer(j)<0)?(cyclic_buffer(j)+win_len+1):cyclic_buffer(j);
}
for (j=0; j<win_len; j++)//for all the bits in the buffer
{
inp = inpt(sur,cyclic_buffer(j));//current bit-decoder output (encoder input)
extrinsic_data(b) = inp;//store the hard output
tmp = cyclic_buffer(j);
cmp = comp(sur, tmp);//SOVA: competitor state (previous)
inc = inpc(sur, tmp);//SOVA: competitor bit output
dif = diff(sur, tmp);//SOVA: corresp. path metric difference
srv = surv(sur, tmp);//SOVA: temporary survivor path state
for (s=j+1; s<=win_len; s++)//check all buffer bits srv and cmp paths
{
if (inp != inc)
{
idx = b-((s-1)-j);//calculate index: [enc.k*(b-(k-1)+j-1)+1:enc.k*(b-(k-1)+j)]
sft(idx) = std::min(sft(idx), dif);//update LLRs for bits that are different
}
if (srv == cmp)
{
break;
}//stop if surv and comp merge (no need to continue)
if (s == win_len)
{
break;
}//stop if the end (otherwise: error)
tmp = cyclic_buffer(s);
inp = inpt(int(srv), tmp);//previous surv bit
inc = inpt(int(cmp), tmp);//previous comp bit
srv = surv(int(srv), tmp);//previous surv state
cmp = surv(int(cmp), tmp);//previous comp state
}
sur = int(surv(sur, cyclic_buffer(j)));//state for the previous surv bit
b--;//update bit index
}
}
// provide soft output with +/- sign:
sft = threshold(sft, SOVA_threshold);
extrinsic_data =
itpp::elem_mult((2.0*extrinsic_data-1.0), SOVA_scaling_factor*sft)-apriori_data;
//free trellis memory
delete[] rsctrellis.prevStates;
delete[] rsctrellis.nextStates;
delete[] rsctrellis.PARout;
delete[] rsctrellis.fm;
}
void SISO::rsc_viterbi(itpp::vec &extrinsic_coded, itpp::vec &extrinsic_data,
const itpp::vec &intrinsic_coded, const itpp::vec &apriori_data, const int &win_len)
/* Soft Input Soft Output module based on Viterbi algorithm
* Output: extrinsic_data - extrinsic information of data bits
* Inputs: intrinsic_coded - intrinsic information of data bits
* apriori_data - a priori information of data bits
* win_len - window length used to represent the code trellis
*
* The implemented algorithm follows M. Kerner and O. Amrani, ''Iterative Decoding
* Using Optimum Soft Input - Hard Output Module`` (2009), in:
* IEEE Transactions on Communications, 57:7(1881-1885)
*/
{
//number of code outputs
int nb_outputs = gen.rows();
//setup internal variables based on RSC trellis
register int i,j,s;
gen_rsctrellis();//trellis generation for 1/2 RSC codes
int nb_states = rsctrellis.numStates;
itpp::Array<itpp::mat> bin_out(2);//contains code output for each initial state and code input
itpp::imat next_state(nb_states,2);//next state in the trellis
for (i=0; i<2; i++)
{
bin_out(i).set_size(nb_states, nb_outputs);
for (j=0; j<nb_states; j++)
{
bin_out(i)(j,0) = double(i);//systematic bit
bin_out(i)(j,1) = rsctrellis.PARout[j+i*nb_states];//parity bit
next_state(j,i) = rsctrellis.nextStates[j+i*nb_states];
}
}
itpp::vec bin_inp("0 1");//binary code inputs
int len = apriori_data.length();//number of decoding steps (total)
//allocate memory for the trellis window
itpp::mat metr(nb_states,win_len+1);//path metric buffer
metr.zeros();
metr += -INFINITY;
metr(0,0) = 0;//initial state => (0,0)
itpp::mat surv(nb_states,win_len+1);//survivor state buffer
surv.zeros();
itpp::mat inpt(nb_states,win_len+1);//survivor input bits buffer (dec. output)
inpt.zeros();
itpp::mat parity(nb_states,win_len+1);//survivor parity bits buffer (dec. output)
parity.zeros();
//decode all the bits
int Cur,Nxt,nxt,sur,b;
itpp::vec buf(nb_outputs);
double llb,mtr;
itpp::vec bin(nb_outputs);
itpp::ivec cyclic_buffer(win_len);
extrinsic_coded.set_size(2*len);//initialize memory for output
extrinsic_data.set_size(len);
for (i = 0; i < len; i++)
{
//indices + precalculations
Cur = i%(win_len+1);//curr trellis (cycl. buf) position
Nxt = (i+1)%(win_len+1);//next trellis (cycl. buf) position
buf = intrinsic_coded(i*nb_outputs,(i+1)*nb_outputs-1);//intrinsic_info portion to be processed
llb = apriori_data(i);//SOVA: apriori_info portion to be processed
metr.set_col(Nxt, -INFINITY*itpp::ones(nb_states));
//forward recursion
for (s = 0; s<nb_states; s++)
{
for (j = 0; j<2; j++)
{
nxt = next_state(s,j);//state after transition
bin = bin_out(j).get_row(s);//transition output (encoder)
mtr = bin*buf+metr(s,Cur);//transition metric
mtr += bin_inp(j)*llb;//add a priori info contribution
if (metr(nxt,Nxt) < mtr)
{
metr(nxt,Nxt) = mtr;//store the metric
surv(nxt,Nxt) = s;//store the survival state
inpt(nxt,Nxt) = bin(0);//j;//store the survival input
parity(nxt,Nxt) = bin(1);//store survival parity bit
}
}
}
//trace backwards
if (i < (win_len-1))
{
continue;
}//proceed if the buffer has been filled
mtr = itpp::max(metr.get_col(Nxt), sur);//find the initial state (max metric)
b = i;//temporary bit index
for (j=0; j<win_len; j++)//indices in a 'cyclic buffer' operation
{
cyclic_buffer(j) = (Nxt-j)%(win_len+1);
cyclic_buffer(j) = (cyclic_buffer(j)<0)?(cyclic_buffer(j)+win_len+1):cyclic_buffer(j);
}
for (j=0; j<win_len; j++)//for all the bits in the buffer
{
extrinsic_data(b) = inpt(sur,cyclic_buffer(j));//current bit-decoder output (encoder input)
extrinsic_coded(2*b) = extrinsic_data(b);
extrinsic_coded(2*b+1) = parity(sur,cyclic_buffer(j));
sur = int(surv(sur, cyclic_buffer(j)));//state for the previous surv bit
b--;//update bit index
}
}
//free trellis memory
delete[] rsctrellis.prevStates;
delete[] rsctrellis.nextStates;
delete[] rsctrellis.PARout;
delete[] rsctrellis.fm;
//output only the hard output if the flag is true
if (Viterbi_hard_output_flag==true)
{
return;
}
// provide soft output (extrinsic information) for data bits
//compute flipped and non-flipped bits positions
itpp::bvec tmp(len);
for (i=0; i<len; i++)
{
tmp(i) = ((apriori_data(i)+intrinsic_coded(2*i))>=0)^(extrinsic_data(i)==1.0);
}
itpp::ivec *ptr_idx_matching = new itpp::ivec;
itpp::ivec *ptr_idx_nonmatching = new itpp::ivec;
*ptr_idx_matching = itpp::find(tmp==itpp::bin(0));
*ptr_idx_nonmatching = itpp::find(tmp==itpp::bin(1));
//Estimated Bit Error Rate
int idx_nonmatching_len = ptr_idx_nonmatching->length();
double error_rate = double(idx_nonmatching_len)/double(len);
//Logarithm of Likelihood Ratio
double LLR;
if (error_rate==0.0)
{
LLR = std::log(double(len));
} else if (error_rate==1.0)
{
LLR = -std::log(double(len));
} else
{
LLR = std::log((1.0-error_rate)/error_rate);
}
for (i=0; i<ptr_idx_matching->length(); i++)
{
extrinsic_data(ptr_idx_matching->get(i)) = Viterbi_scaling_factor[0]*
(2.0*extrinsic_data(ptr_idx_matching->get(i))-1.0)*LLR;
}
for (i=0; i<idx_nonmatching_len; i++)
{
extrinsic_data(ptr_idx_nonmatching->get(i)) = Viterbi_scaling_factor[1]*
(2.0*extrinsic_data(ptr_idx_nonmatching->get(i))-1.0)*LLR;
}
//extrinsic information
extrinsic_data -= apriori_data;
//provide soft output for coded bits
tmp.set_length(2*len);
for (i=0; i<(2*len); i++)
{
tmp(i) = (intrinsic_coded(i)>=0)^(extrinsic_coded(i)==1.0);
}
delete ptr_idx_matching;//free old memory
delete ptr_idx_nonmatching;
ptr_idx_matching = new itpp::ivec;//allocate memory for new vectors
ptr_idx_nonmatching = new itpp::ivec;
*ptr_idx_matching = itpp::find(tmp==itpp::bin(0));
*ptr_idx_nonmatching = itpp::find(tmp==itpp::bin(1));
//Estimated Bit Error Rate
idx_nonmatching_len = ptr_idx_nonmatching->length();
error_rate = double(idx_nonmatching_len)/double(2*len);
//Logarithm of Likelihood Ratio
if (error_rate==0.0)
{
LLR = std::log(double(2*len));
} else if (error_rate==1.0)
{
LLR = -std::log(double(2*len));
} else
{
LLR = std::log((1.0-error_rate)/error_rate);
}
for (i=0; i<ptr_idx_matching->length(); i++)
{
extrinsic_coded(ptr_idx_matching->get(i)) = Viterbi_scaling_factor[0]*
(2.0*extrinsic_coded(ptr_idx_matching->get(i))-1.0)*LLR;
}
for (i=0; i<idx_nonmatching_len; i++)
{
extrinsic_coded(ptr_idx_nonmatching->get(i)) = Viterbi_scaling_factor[1]*
(2.0*extrinsic_coded(ptr_idx_nonmatching->get(i))-1.0)*LLR;
}
delete ptr_idx_matching;
delete ptr_idx_nonmatching;
//extrinsic information
extrinsic_coded -= intrinsic_coded;
}
}
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
itcompat.h
IT++ compatibility types and functions.
itpp::find
ivec find(const bvec &invector)
Return a integer vector with indicies where bvec == 1.
Definition: specmat.cpp:40
itpp
itpp namespace
Definition: itmex.h:36
itpp::bin
Binary arithmetic (boolean) class.
Definition: binary.h:56
itpp::elem_mult
Mat< Num_T > elem_mult(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise multiplication of two matrices.
Definition: mat.h:1582
itpp::bin2dec
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::Array
General array class.
Definition: array.h:42
itpp::log_add
double log_add(double log_a, double log_b)
Safe substitute for log(exp(log_a) + exp(log_b))
Definition: log_exp.cpp:40
itpp::min
T min(const Vec< T > &in)
Minimum value of vector.
Definition: min_max.h:125
itpp::exp
vec exp(const vec &x)
Exp of the elements of a vector x.
Definition: log_exp.h:155
itpp::max
T max(const Vec< T > &v)
Maximum value of vector.
Definition: min_max.h:45
siso.h
Definitions for Soft Input Soft Output (SISO) modules class.
itpp::dec2bin
ITPP_EXPORT bvec dec2bin(int length, int index)
Convert a decimal int index to bvec using length bits in the representation.
itpp::ones
ITPP_EXPORT vec ones(int size)
A float vector of ones.
itpp::log
vec log(const vec &x)
The natural logarithm of the elements.
Definition: log_exp.h:241

Generated on Mon Mar 23 2020 06:18:03 for IT++ by Doxygen 1.8.17