#include <vector>

namespace gr {

#define MAXLOG 1e7

std::vector<int> d_rowpolys;

std::vector<int> d_colpolys;

unsigned int d_krow;

unsigned int d_kcol;

int d_bval;

int d_qval;

int d_max_iter;

int d_decoder_type;

// store the state transitions & outputs

int rowNumStates;

std::vector< std::vector<int> > rowOutputs;

int colNumStates;

std::vector< std::vector<int> > colOutputs;

std::vector< std::vector<int> > colNextStates;

int rowEncoder_K;

int rowEncoder_n;

int rowEncoder_m;

int colEncoder_m;

int outputSize;

int inputSize;

uint32_t codeword_M;

uint32_t codeword_N;

int mInit, nInit;

// memory allocated for processing

int inputSizeWithPad;

std::vector< std::vector<float> > channel_llr;

std::vector< std::vector<float> > Z;

std::vector<float> extrinsic_info;

std::vector<float> output_u_cols;

std::vector<float> output_c_rows;

std::vector<float> output_c_cols;

uint32_t numInitLoadIter;

int numInitRemaining;

int output_c_col_idx;

bool earlyExit;

FILE *fp;

// soft input soft output decoding

int mm_row, max_states_row, num_symbols_row;

std::vector< std::vector<float> > beta_row;

std::vector<float> num_llr_c_row;

std::vector<float> den_llr_c_row;

void siso_decode_row();

int mm_col, max_states_col, num_symbols_col;

std::vector< std::vector<float> > beta_col;

std::vector<float> alpha_prime_col;

// Computes the branch metric used for decoding, returning a metric between the

// hypothetical symbol and received vector

float gamma(const std::vector<float> rec_array, const int symbol);

float (tpc_decoder::*max_star)(const float, const float);

float linear_log_map(const float delta1, const float delta2);

float max_log_map(const float delta1, const float delta2);

float constant_log_map(const float delta1, const float delta2);

float log_map_lut_correction(const float delta1, const float delta2);

float log_map_cfunction_correction(const float delta1, const float delta2);

template <typename T> static int sgn(T val);

public:

#include <gnuradio/fec/maxstar.h>

generic_decoder::sptr

tpc_decoder::make(std::vector<int> row_polys, std::vector<int> col_polys, int krow, int kcol, int bval, int qval, int max_iter, int decoder_type)

// first we operate on data chunks of get_input_size()

// TODO: should we verify this and throw an error if it doesn't match? YES

// hwo do we do that?

rowEncoder_K = ceil(log(d_rowpolys[0])/log(2)); // rowEncoder_K is the constraint length of the row encoder polynomial

rowEncoder_n = d_rowpolys.size();

rowEncoder_m = rowEncoder_K - 1;

colEncoder_K = ceil(log(d_colpolys[0])/log(2)); // colEncoder_K is the constraint length of the col encoder polynomial

colEncoder_n = d_colpolys.size();

colEncoder_m = colEncoder_K - 1;

// calculate the input and output sizes

inputSize = ((d_krow+rowEncoder_m)*rowEncoder_n)*((d_kcol+colEncoder_m)*colEncoder_n) - d_bval;

outputSize = (d_krow*d_kcol - (d_bval+d_qval));

//DEBUG_PRINT("inputSize=%d outputSize=%d\n", inputSize, outputSize);

//fp = fopen("c_decoder_output.txt", "w");

rowNumStates = 1 << (rowEncoder_m); // 2^(row_mm)

colNumStates = 1 << (colEncoder_m); // 2^(col_mm)

rowOutputs.resize(2, std::vector<int>(rowNumStates,0));

rowNextStates.resize(2, std::vector<int>(rowNumStates,0));

colOutputs.resize(2, std::vector<int>(colNumStates,0));

colNextStates.resize(2, std::vector<int>(colNumStates,0));;

// precalculate the state transition matrix for the row polynomial

rowOutputs, rowNextStates);

// precalculate the state transition matrix for the column polynomial

colOutputs, colNextStates);

codeword_M = d_kcol + colEncoder_K - 1;

codeword_N = d_krow + rowEncoder_K - 1;

// pre-allocate memory we use for encoding

inputSizeWithPad = inputSize + d_bval;

channel_llr.resize(codeword_M, std::vector<float>(codeword_N, 0));

Z.resize(codeword_M, std::vector<float>(codeword_N, 0));

extrinsic_info.resize(codeword_M*codeword_N, 0);

input_u_rows.resize(d_krow, 0);

input_u_cols.resize(d_kcol, 0);

input_c_rows.resize(codeword_N, 0);

output_u_rows.resize(d_krow, 0);

output_u_cols.resize(d_kcol, 0);

output_c_rows.resize(codeword_N, 0);

output_c_cols.resize(codeword_M*codeword_N, 0);

output_c_col_idx = 0;

// setup the max_star function based on decoder type

switch(d_decoder_type) {

case 0:

max_star = &tpc_decoder::linear_log_map;

break;

}

// declare the reverse sweep trellis

// the beta vector is logically layed out in memory as follows, assuming the

// following values (for educational purposes)

// as well as some comments being added to help understand

// exactly what is going on w/ the siso decoder

void tpc_decoder::siso_decode_row() {

int LL, state, k, ii, symbol, mask;

float app_in, delta1, delta2;

LL = input_u_rows.size(); // code length

float num_llr_u;

// log-likelihood ratio of the uncoded bit being a 0

float den_llr_u;

// initialize beta_row trellis

// this initialization is saying that the likelihood that the reverse sweep

// starts at state=0 is 100%, b/c state 1, 2, 3's likelihood's are basically -inf

// }

// filling w/ 0xCC yields a value close to -MAXLOG, and memset is faster than for loops

memset(&beta_row[LL+rowEncoder_K-1][1], 0xCC, sizeof(float)*(max_states_row-1));

// initialize alpha_prime_row (current time instant), alpha_prime_row then gets updated

// by shifting in alpha_row at the end of each time instant of processing

// alpha_row needs to get initialized at the beginning of each processing loop, so we

// alpha_prime_row[state] = -MAXLOG;

// }

memset(&alpha_prime_row[1], 0xCC, sizeof(float)*(max_states_row-1));

// compute the beta_row matrix first, which is the reverse sweep (hence we start at the last

// time instant-1 and work our way back to t=0). we start at last time instant-1 b/c

// we already filled in beta_row values for the last time instant, forcing the trellis to

else {

app_in = 0;

}

// get the input associated w/ this time instant

memcpy(&rec_array_row[0], &input_c_rows[rowEncoder_n*k], sizeof(float)*rowEncoder_n);

// DEBUG_PRINT("k=%d\n", k);

// DEBUG_PRINT_F(fp, "k=%d\n", k);

// DEBUG_PRINT("rec_array -->\n");

// DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT(&rec_array_row[0], rec_array_row.size());

// DEBUG_PRINT_F(fp, "rec_array -->\n");

// DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT_F(fp, &rec_array_row[0], rec_array_row.size());

// represents the likelihood that this input corresponds to

// each of the possible symbols

for(symbol=0; symbol<num_symbols_row; symbol++) {

metric_c_row[symbol] = gamma(rec_array_row, symbol);

}

// DEBUG_PRINT("metric_c -->\n");

// DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT(&metric_c_row[0], metric_c_row.size());

// DEBUG_PRINT_F(fp, "metric_c -->\n");

// DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT_F(fp, &metric_c_row[0], metric_c_row.size());

// step through all states -- populating the beta_row values for each node

// in the trellis diagram as shown above w/ the maximum-likelihood

// of this current node coming from the previous node

delta2 = beta_row[k+1][rowNextStates[1][state]] + metric_c_row[rowOutputs[1][state]] + app_in;

// update beta_row

beta_row[k][state] = (*this.*max_star)(delta1, delta2);

// DEBUG_PRINT("delta1=%f delta2=%f beta=%f\n", delta1, delta2, beta_row[k][state]);

// DEBUG_PRINT_F(fp, "delta1=%f delta2=%f beta=%f\n", delta1, delta2, beta_row[k][state]);

}

// normalize beta_row

for (state=1;state<max_states_row;state++) {

beta_row[k][state] = beta_row[k][state] - beta_row[k][0];

}

beta_row[k][0] = 0;

}

// compute the forward sweep (alpha_row values), and update the llr's

// notice that we start at time index=1, b/c time index=0 has already been

// initialized, and we are forcing the trellis to start from state=0

// initialize the llr's for the uncoded bits

num_llr_u = -MAXLOG;

den_llr_u = -MAXLOG;

// intialize alpha_row

// for (state=0;state<max_states_row;state++) {

// alpha_row[state] = -MAXLOG;

num_llr_c_row[ii] = -MAXLOG;

rec_array_row[ii] = input_c_rows[rowEncoder_n*(k-1)+ii];

}

// represents the likelihood that this input corresponds to

// each of the possible symbols

for(symbol=0; symbol<num_symbols_row; symbol++) {

metric_c_row[ii] = gamma(rec_array_row, symbol);

}

// compute the alpha_row vector

// to understand the loop below, we need to think about the forward trellis.

// we know that any node, at any time instant in the trellis diagram can have

// previous time instance can transition into this current node, based on the

// polynomial). SO, in the loop below, delta1 represents the transition from

// the previous time instant for input (either 0 or 1) and the current state

delta2 = alpha_row[rowNextStates[1][state]];

alpha_row[rowNextStates[1][state]] = (*this.*max_star)(delta1, delta2);

}

// compute the llr's

for (state=0;state<max_states_row;state++) {

// data 0 branch (departing)

}

mask = mask>>1;

}

// data 1 branch (departing)

delta1 = alpha_prime_row[state] + metric_c_row[rowOutputs[1][state]] + beta_row[k][rowNextStates[1][state]] + app_in;

// the information bit

// shift alpha_row back to alpha_prime_row

alpha_prime_row[state] = alpha_row[state] - alpha_row[0];

}

if (k-1<LL) {

output_u_rows[k-1] = num_llr_u - den_llr_u;

}

float tpc_decoder::linear_log_map(const float delta1, const float delta2) {

float diff;

diff = delta2 - delta1;

if ( diff > TJIAN )

float tpc_decoder::constant_log_map(const float delta1, const float delta2) {

// Return maximum of delta1 and delta2

// and in correction value if |delta1-delta2| < TVALUE

diff = delta2 - delta1;

if ( diff > TVALUE )

float tpc_decoder::log_map_lut_correction(const float delta1, const float delta2) {

float diff;

diff = (float) fabs( delta2 - delta1 );

if (delta1 > delta2) {

if (diff > BOUNDARY8 )

return( delta1 );

return( delta1 + VALUE5 + SLOPE5*(diff-BOUNDARY5) );

else

return( delta1 + VALUE4 + SLOPE4*(diff-BOUNDARY4) );

} else {

if (diff > BOUNDARY2 ) {

if ( diff > BOUNDARY3 )

return( delta2 + VALUE5 + SLOPE5*(diff-BOUNDARY5) );

else

return( delta2 + VALUE4 + SLOPE4*(diff-BOUNDARY4) );

} else {

if (diff > BOUNDARY2 ) {

if ( diff > BOUNDARY3 )

float tpc_decoder::log_map_cfunction_correction(const float delta1, const float delta2) {

// Use C-function calls to compute the correction function

if (delta1 > delta2) {

} else {

}

}

int ii;

int mask;

int nn = rx_array.size();

mask = 1;

for (ii=0;ii<nn;ii++) {

if (symbol&mask)

// DEBUG_PRINT("nn=%d symbol=%d rm = %f\n", nn, symbol, rm);

// DEBUG_PRINT_F(fp, "nn=%d symbol=%d rm = %f\n", nn, symbol, rm);

return(rm);

}

void tpc_decoder::generic_work(void *inBuffer, void *outBuffer) {

const float *inPtr = (const float *) inBuffer;

unsigned char *out = (unsigned char *) outBuffer;

unsigned int m, n, ii;

int iter;

for(ii=0; ii<numInitLoadIter; ii++) {

memset(&channel_llr[ii][0], 0, sizeof(float)*codeword_N);

memset(&Z[ii][0], 0, sizeof(float)*codeword_N);

memset(&extrinsic_info[0], 0, sizeof(float)*extrinsic_info.size());

//DEBUG_PRINT("Starting TURBO Decoding\n");

for(iter=0; iter<d_max_iter; iter++) {

//DEBUG_PRINT("Turbo Iter=%d\n", iter+1);

//DEBUG_PRINT_F(fp, "Turbo Iter=%d\n", iter+1);

//DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT(&input_c_rows[0], codeword_N);

//DEBUG_PRINT_F(fp, "input_c_rows -->\n");

//DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT_F(fp, &input_c_rows[0], codeword_N);

// call siso decode

siso_decode_row();

//DEBUG_PRINT("output_u_rows -->\n");

//DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT(&output_u_rows[0], output_u_rows.size());

//DEBUG_PRINT_F(fp, "output_u_rows -->\n");

//DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT(&output_c_rows[0], output_c_rows.size());

//DEBUG_PRINT_F(fp, "output_c_rows -->\n");

//DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT_F(fp, &output_c_rows[0], output_c_rows.size());

// copy the output coded back into Z, so we can feed it back through the decoder

// for more iterations

volk_32f_x2_subtract_32f(&Z[m][0], &output_c_rows[0], &extrinsic_info[m*codeword_N], codeword_N);

}

// decode each col

earlyExit = true;

output_c_col_idx = 0;

for(ii=0; ii<codeword_M; ii++) {

input_c_cols[ii] = Z[ii][n];

}

//DEBUG_PRINT("input_c_cols -->\n");

//DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT(&input_c_cols[0], codeword_M);

//DEBUG_PRINT_F(fp, "input_c_cols -->\n");

//DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT_F(fp, &input_c_cols[0], codeword_M);

// call siso decode

siso_decode_col();

//DEBUG_PRINT("output_u_cols -->\n");

//DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT(&output_u_cols[0], output_u_cols.size());

//DEBUG_PRINT_F(fp, "output_u_cols -->\n");

//DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT(&output_c_cols[output_c_col_idx], codeword_M);

//DEBUG_PRINT_F(fp, "output_c_cols -->\n");

//DEBUG_PRINT_FLOAT_ARRAY_AS_FLOAT_F(fp, &output_c_cols[output_c_col_idx], codeword_M);

// create the extrinsic_info vector, which subtracts from input_c_cols to prevent feedback

//DEBUG_PRINT("extrinsic_info -->\n");

//DEBUG_PRINT_F(fp, "extrinsic_info -->\n");

for(ii=0; ii<codeword_M; ii++) {

extrinsic_info[ii*codeword_N+n] = output_c_cols[output_c_col_idx+ii] - input_c_cols[ii];

if(earlyExit) {

if(sgn(output_c_cols[output_c_col_idx+ii])!=sgn(input_c_cols[ii])) {

earlyExit = false;

}

//DEBUG_PRINT("\n");

//DEBUG_PRINT_F(fp, "\n");

output_c_col_idx += codeword_M;

}

if(earlyExit) break;

}

}

}

//DEBUG_PRINT(("Output\n"));

//DEBUG_PRINT_UCHAR_ARRAY(out, outputSize);

//DEBUG_PRINT_F(fp, "Output\n");

return "none";

}

fclose(fp);

}