+
+/* Returns a random number drawn from a normal distribution
+ * with mean of 0 and variance of 1
+ * Marsaglia algorithm
+ */
+static double
+randn(int n)
+{
+ double S, Z, U1, U2, u, v, fac;
+
+ do {
+ U1 = (double)rand() / RAND_MAX;
+ U2 = (double)rand() / RAND_MAX;
+ u = 2. * U1 - 1.;
+ v = 2. * U2 - 1.;
+ S = u * u + v * v;
+ } while (S >= 1 || S == 0);
+ fac = sqrt(-2. * log(S) / S);
+ Z = (n % 2) ? u * fac : v * fac;
+ return Z;
+}
+
+static inline double
+maxstar(double A, double B)
+{
+ if (fabs(A - B) > 5)
+ return RTE_MAX(A, B);
+ else
+ return RTE_MAX(A, B) + log1p(exp(-fabs(A - B)));
+}
+
+/*
+ * Generate Qm LLRS for Qm==8
+ * Modulation, AWGN and LLR estimation from max log development
+ */
+static void
+gen_qm8_llr(int8_t *llrs, uint32_t i, double N0, double llr_max)
+{
+ int qm = 8;
+ int qam = 256;
+ int m, k;
+ double I, Q, p0, p1, llr_, b[qm], log_syml_prob[qam];
+ /* 5.1.4 of TS38.211 */
+ const double symbols_I[256] = {
+ 5, 5, 7, 7, 5, 5, 7, 7, 3, 3, 1, 1, 3, 3, 1, 1, 5,
+ 5, 7, 7, 5, 5, 7, 7, 3, 3, 1, 1, 3, 3, 1, 1, 11,
+ 11, 9, 9, 11, 11, 9, 9, 13, 13, 15, 15, 13, 13,
+ 15, 15, 11, 11, 9, 9, 11, 11, 9, 9, 13, 13, 15,
+ 15, 13, 13, 15, 15, 5, 5, 7, 7, 5, 5, 7, 7, 3, 3,
+ 1, 1, 3, 3, 1, 1, 5, 5, 7, 7, 5, 5, 7, 7, 3, 3, 1,
+ 1, 3, 3, 1, 1, 11, 11, 9, 9, 11, 11, 9, 9, 13, 13,
+ 15, 15, 13, 13, 15, 15, 11, 11, 9, 9, 11, 11, 9, 9,
+ 13, 13, 15, 15, 13, 13, 15, 15, -5, -5, -7, -7, -5,
+ -5, -7, -7, -3, -3, -1, -1, -3, -3, -1, -1, -5, -5,
+ -7, -7, -5, -5, -7, -7, -3, -3, -1, -1, -3, -3,
+ -1, -1, -11, -11, -9, -9, -11, -11, -9, -9, -13,
+ -13, -15, -15, -13, -13, -15, -15, -11, -11, -9,
+ -9, -11, -11, -9, -9, -13, -13, -15, -15, -13,
+ -13, -15, -15, -5, -5, -7, -7, -5, -5, -7, -7, -3,
+ -3, -1, -1, -3, -3, -1, -1, -5, -5, -7, -7, -5, -5,
+ -7, -7, -3, -3, -1, -1, -3, -3, -1, -1, -11, -11,
+ -9, -9, -11, -11, -9, -9, -13, -13, -15, -15, -13,
+ -13, -15, -15, -11, -11, -9, -9, -11, -11, -9, -9,
+ -13, -13, -15, -15, -13, -13, -15, -15};
+ const double symbols_Q[256] = {
+ 5, 7, 5, 7, 3, 1, 3, 1, 5, 7, 5, 7, 3, 1, 3, 1, 11,
+ 9, 11, 9, 13, 15, 13, 15, 11, 9, 11, 9, 13, 15, 13,
+ 15, 5, 7, 5, 7, 3, 1, 3, 1, 5, 7, 5, 7, 3, 1, 3, 1,
+ 11, 9, 11, 9, 13, 15, 13, 15, 11, 9, 11, 9, 13,
+ 15, 13, 15, -5, -7, -5, -7, -3, -1, -3, -1, -5,
+ -7, -5, -7, -3, -1, -3, -1, -11, -9, -11, -9, -13,
+ -15, -13, -15, -11, -9, -11, -9, -13, -15, -13,
+ -15, -5, -7, -5, -7, -3, -1, -3, -1, -5, -7, -5,
+ -7, -3, -1, -3, -1, -11, -9, -11, -9, -13, -15,
+ -13, -15, -11, -9, -11, -9, -13, -15, -13, -15, 5,
+ 7, 5, 7, 3, 1, 3, 1, 5, 7, 5, 7, 3, 1, 3, 1, 11,
+ 9, 11, 9, 13, 15, 13, 15, 11, 9, 11, 9, 13, 15,
+ 13, 15, 5, 7, 5, 7, 3, 1, 3, 1, 5, 7, 5, 7, 3, 1,
+ 3, 1, 11, 9, 11, 9, 13, 15, 13, 15, 11, 9, 11, 9,
+ 13, 15, 13, 15, -5, -7, -5, -7, -3, -1, -3, -1,
+ -5, -7, -5, -7, -3, -1, -3, -1, -11, -9, -11, -9,
+ -13, -15, -13, -15, -11, -9, -11, -9, -13, -15,
+ -13, -15, -5, -7, -5, -7, -3, -1, -3, -1, -5, -7,
+ -5, -7, -3, -1, -3, -1, -11, -9, -11, -9, -13, -15,
+ -13, -15, -11, -9, -11, -9, -13, -15, -13, -15};
+ /* Average constellation point energy */
+ N0 *= 170.0;
+ for (k = 0; k < qm; k++)
+ b[k] = llrs[qm * i + k] < 0 ? 1.0 : 0.0;
+ /* 5.1.4 of TS38.211 */
+ I = (1 - 2 * b[0]) * (8 - (1 - 2 * b[2]) *
+ (4 - (1 - 2 * b[4]) * (2 - (1 - 2 * b[6]))));
+ Q = (1 - 2 * b[1]) * (8 - (1 - 2 * b[3]) *
+ (4 - (1 - 2 * b[5]) * (2 - (1 - 2 * b[7]))));
+ /* AWGN channel */
+ I += sqrt(N0 / 2) * randn(0);
+ Q += sqrt(N0 / 2) * randn(1);
+ /*
+ * Calculate the log of the probability that each of
+ * the constellation points was transmitted
+ */
+ for (m = 0; m < qam; m++)
+ log_syml_prob[m] = -(pow(I - symbols_I[m], 2.0)
+ + pow(Q - symbols_Q[m], 2.0)) / N0;
+ /* Calculate an LLR for each of the k_64QAM bits in the set */
+ for (k = 0; k < qm; k++) {
+ p0 = -999999;
+ p1 = -999999;
+ /* For each constellation point */
+ for (m = 0; m < qam; m++) {
+ if ((m >> (qm - k - 1)) & 1)
+ p1 = maxstar(p1, log_syml_prob[m]);
+ else
+ p0 = maxstar(p0, log_syml_prob[m]);
+ }
+ /* Calculate the LLR */
+ llr_ = p0 - p1;
+ llr_ *= (1 << ldpc_llr_decimals);
+ llr_ = round(llr_);
+ if (llr_ > llr_max)
+ llr_ = llr_max;
+ if (llr_ < -llr_max)
+ llr_ = -llr_max;
+ llrs[qm * i + k] = (int8_t) llr_;
+ }
+}
+
+
+/*
+ * Generate Qm LLRS for Qm==6
+ * Modulation, AWGN and LLR estimation from max log development
+ */
+static void
+gen_qm6_llr(int8_t *llrs, uint32_t i, double N0, double llr_max)
+{
+ int qm = 6;
+ int qam = 64;
+ int m, k;
+ double I, Q, p0, p1, llr_, b[qm], log_syml_prob[qam];
+ /* 5.1.4 of TS38.211 */
+ const double symbols_I[64] = {
+ 3, 3, 1, 1, 3, 3, 1, 1, 5, 5, 7, 7, 5, 5, 7, 7,
+ 3, 3, 1, 1, 3, 3, 1, 1, 5, 5, 7, 7, 5, 5, 7, 7,
+ -3, -3, -1, -1, -3, -3, -1, -1, -5, -5, -7, -7,
+ -5, -5, -7, -7, -3, -3, -1, -1, -3, -3, -1, -1,
+ -5, -5, -7, -7, -5, -5, -7, -7};
+ const double symbols_Q[64] = {
+ 3, 1, 3, 1, 5, 7, 5, 7, 3, 1, 3, 1, 5, 7, 5, 7,
+ -3, -1, -3, -1, -5, -7, -5, -7, -3, -1, -3, -1,
+ -5, -7, -5, -7, 3, 1, 3, 1, 5, 7, 5, 7, 3, 1, 3, 1,
+ 5, 7, 5, 7, -3, -1, -3, -1, -5, -7, -5, -7,
+ -3, -1, -3, -1, -5, -7, -5, -7};
+ /* Average constellation point energy */
+ N0 *= 42.0;
+ for (k = 0; k < qm; k++)
+ b[k] = llrs[qm * i + k] < 0 ? 1.0 : 0.0;
+ /* 5.1.4 of TS38.211 */
+ I = (1 - 2 * b[0])*(4 - (1 - 2 * b[2]) * (2 - (1 - 2 * b[4])));
+ Q = (1 - 2 * b[1])*(4 - (1 - 2 * b[3]) * (2 - (1 - 2 * b[5])));
+ /* AWGN channel */
+ I += sqrt(N0 / 2) * randn(0);
+ Q += sqrt(N0 / 2) * randn(1);
+ /*
+ * Calculate the log of the probability that each of
+ * the constellation points was transmitted
+ */
+ for (m = 0; m < qam; m++)
+ log_syml_prob[m] = -(pow(I - symbols_I[m], 2.0)
+ + pow(Q - symbols_Q[m], 2.0)) / N0;
+ /* Calculate an LLR for each of the k_64QAM bits in the set */
+ for (k = 0; k < qm; k++) {
+ p0 = -999999;
+ p1 = -999999;
+ /* For each constellation point */
+ for (m = 0; m < qam; m++) {
+ if ((m >> (qm - k - 1)) & 1)
+ p1 = maxstar(p1, log_syml_prob[m]);
+ else
+ p0 = maxstar(p0, log_syml_prob[m]);
+ }
+ /* Calculate the LLR */
+ llr_ = p0 - p1;
+ llr_ *= (1 << ldpc_llr_decimals);
+ llr_ = round(llr_);
+ if (llr_ > llr_max)
+ llr_ = llr_max;
+ if (llr_ < -llr_max)
+ llr_ = -llr_max;
+ llrs[qm * i + k] = (int8_t) llr_;
+ }
+}
+
+/*
+ * Generate Qm LLRS for Qm==4
+ * Modulation, AWGN and LLR estimation from max log development
+ */
+static void
+gen_qm4_llr(int8_t *llrs, uint32_t i, double N0, double llr_max)
+{
+ int qm = 4;
+ int qam = 16;
+ int m, k;
+ double I, Q, p0, p1, llr_, b[qm], log_syml_prob[qam];
+ /* 5.1.4 of TS38.211 */
+ const double symbols_I[16] = {1, 1, 3, 3, 1, 1, 3, 3,
+ -1, -1, -3, -3, -1, -1, -3, -3};
+ const double symbols_Q[16] = {1, 3, 1, 3, -1, -3, -1, -3,
+ 1, 3, 1, 3, -1, -3, -1, -3};
+ /* Average constellation point energy */
+ N0 *= 10.0;
+ for (k = 0; k < qm; k++)
+ b[k] = llrs[qm * i + k] < 0 ? 1.0 : 0.0;
+ /* 5.1.4 of TS38.211 */
+ I = (1 - 2 * b[0]) * (2 - (1 - 2 * b[2]));
+ Q = (1 - 2 * b[1]) * (2 - (1 - 2 * b[3]));
+ /* AWGN channel */
+ I += sqrt(N0 / 2) * randn(0);
+ Q += sqrt(N0 / 2) * randn(1);
+ /*
+ * Calculate the log of the probability that each of
+ * the constellation points was transmitted
+ */
+ for (m = 0; m < qam; m++)
+ log_syml_prob[m] = -(pow(I - symbols_I[m], 2.0)
+ + pow(Q - symbols_Q[m], 2.0)) / N0;
+ /* Calculate an LLR for each of the k_64QAM bits in the set */
+ for (k = 0; k < qm; k++) {
+ p0 = -999999;
+ p1 = -999999;
+ /* For each constellation point */
+ for (m = 0; m < qam; m++) {
+ if ((m >> (qm - k - 1)) & 1)
+ p1 = maxstar(p1, log_syml_prob[m]);
+ else
+ p0 = maxstar(p0, log_syml_prob[m]);
+ }
+ /* Calculate the LLR */
+ llr_ = p0 - p1;
+ llr_ *= (1 << ldpc_llr_decimals);
+ llr_ = round(llr_);
+ if (llr_ > llr_max)
+ llr_ = llr_max;
+ if (llr_ < -llr_max)
+ llr_ = -llr_max;
+ llrs[qm * i + k] = (int8_t) llr_;
+ }
+}
+
+static void
+gen_qm2_llr(int8_t *llrs, uint32_t j, double N0, double llr_max)
+{
+ double b, b1, n;
+ double coeff = 2.0 * sqrt(N0);
+
+ /* Ignore in vectors rare quasi null LLRs not to be saturated */
+ if (llrs[j] < 8 && llrs[j] > -8)
+ return;
+
+ /* Note don't change sign here */
+ n = randn(j % 2);
+ b1 = ((llrs[j] > 0 ? 2.0 : -2.0)
+ + coeff * n) / N0;
+ b = b1 * (1 << ldpc_llr_decimals);
+ b = round(b);
+ if (b > llr_max)
+ b = llr_max;
+ if (b < -llr_max)
+ b = -llr_max;
+ llrs[j] = (int8_t) b;
+}
+
+/* Generate LLR for a given SNR */
+static void
+generate_llr_input(uint16_t n, struct rte_bbdev_op_data *inputs,
+ struct rte_bbdev_dec_op *ref_op)
+{
+ struct rte_mbuf *m;
+ uint16_t qm;
+ uint32_t i, j, e, range;
+ double N0, llr_max;
+
+ e = ref_op->ldpc_dec.cb_params.e;
+ qm = ref_op->ldpc_dec.q_m;
+ llr_max = (1 << (ldpc_llr_size - 1)) - 1;
+ range = e / qm;
+ N0 = 1.0 / pow(10.0, get_snr() / 10.0);
+
+ for (i = 0; i < n; ++i) {
+ m = inputs[i].data;
+ int8_t *llrs = rte_pktmbuf_mtod_offset(m, int8_t *, 0);
+ if (qm == 8) {
+ for (j = 0; j < range; ++j)
+ gen_qm8_llr(llrs, j, N0, llr_max);
+ } else if (qm == 6) {
+ for (j = 0; j < range; ++j)
+ gen_qm6_llr(llrs, j, N0, llr_max);
+ } else if (qm == 4) {
+ for (j = 0; j < range; ++j)
+ gen_qm4_llr(llrs, j, N0, llr_max);
+ } else {
+ for (j = 0; j < e; ++j)
+ gen_qm2_llr(llrs, j, N0, llr_max);
+ }
+ }
+}
+