| /* |
| * SnsQuantization.cpp |
| * |
| * Copyright 2021 HIMSA II K/S - www.himsa.com. Represented by EHIMA - |
| * www.ehima.com |
| * |
| * Licensed under the Apache License, Version 2.0 (the "License"); |
| * you may not use this file except in compliance with the License. |
| * You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| |
| #include "SnsQuantization.hpp" |
| |
| #include <cmath> |
| #include <limits> |
| |
| #include "SnsQuantizationTables.hpp" |
| |
| namespace Lc3Enc { |
| |
| SnsQuantization::SnsQuantization() |
| : ind_LF(0), |
| ind_HF(0), |
| shape_j(0), |
| Gind(0), |
| LS_indA(0), |
| LS_indB(0), |
| index_joint_j(0), |
| idxA(0), |
| idxB(0), |
| datapoints(nullptr) {} |
| |
| SnsQuantization::~SnsQuantization() {} |
| |
| void SnsQuantization::run(const double* const scf) { |
| // 3.3.7.3.2 Stage 1 (d09r02_F2F) |
| if (!std::numeric_limits<double>::has_infinity) { |
| // we never should reach this line, but really want to make sure to |
| // go further when the assumed infinity() feature is not supported |
| return; |
| } |
| double dMSE_LF_min = std::numeric_limits<double>::infinity(); |
| double dMSE_HF_min = std::numeric_limits<double>::infinity(); |
| for (uint8_t i = 0; i < 32; i++) { |
| double dMSE_LF_i = 0; |
| double dMSE_HF_i = 0; |
| for (uint8_t n = 0; n < 8; n++) { |
| dMSE_LF_i += (scf[n] - LFCB[i][n]) * (scf[n] - LFCB[i][n]); |
| dMSE_HF_i += (scf[n + 8] - HFCB[i][n]) * (scf[n + 8] - HFCB[i][n]); |
| } |
| if (dMSE_LF_i < dMSE_LF_min) { |
| ind_LF = i; |
| dMSE_LF_min = dMSE_LF_i; |
| } |
| if (dMSE_HF_i < dMSE_HF_min) { |
| ind_HF = i; |
| dMSE_HF_min = dMSE_HF_i; |
| } |
| } |
| |
| // The first stage vector is composed as: |
| //π π‘1(π) = πΏπΉπΆπ΅πππ_πΏπΉ (π), πππ π = [0 β¦ 7], # (33) |
| //π π‘1(π + 8) = π»πΉπΆπ΅πππ_π»πΉ(π), πππ π = [0 β¦ 7], # (34) |
| double st1[16]; |
| for (uint8_t n = 0; n < 8; n++) { |
| st1[n] = LFCB[ind_LF][n]; |
| st1[n + 8] = HFCB[ind_HF][n]; |
| } |
| double r1[16]; |
| for (uint8_t n = 0; n < 16; n++) { |
| r1[n] = scf[n] - st1[n]; |
| } |
| |
| // 3.3.7.3.3 Stage 2 (d09r02_F2F) |
| // 3.3.7.3.3.3 Stage 2 target preparation (d09r02_F2F) |
| for (uint8_t n = 0; n < 16; n++) { |
| t2rot[n] = 0; |
| for (uint8_t row = 0; row < 16; row++) { |
| t2rot[n] += r1[row] * D[row][n]; |
| } |
| } |
| |
| // 3.3.7.3.3.4 Shape candidates (d09r02_F2F) |
| // 3.3.7.3.3.5 Stage 2 shape search (d09r02_F2F) |
| |
| // (3) Shape name: βoutlier_farβ |
| // Scale factor set A: {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15} |
| // Scale factor set B: Empty set |
| // Pulse configuration, Set A, PVQ(NA, KA): PVQ(16, 6) |
| // Pulse configuration, Set B, PVQ(NB, KB): Empty |
| |
| // step 1 -> y3,start |
| // Project to or below pyramid N=16, K=6, |
| // (and update energy energyy and correlation corrxy terms to reflect |
| // the pulses present in y3, start ) |
| uint8_t K = 6; |
| uint8_t N = 16; |
| double proj_fac = 0; |
| double abs_x[N]; |
| for (uint8_t n = 0; n < N; n++) { |
| abs_x[n] = fabs(t2rot[n]); |
| proj_fac += abs_x[n]; |
| } |
| proj_fac = (K - 1) / proj_fac; |
| |
| uint8_t k = 0; |
| double corr_xy = 0; |
| double energy_y = 0; |
| for (uint8_t n = 0; n < N; n++) { |
| sns_Y3[n] = floor(abs_x[n] * proj_fac); |
| if (sns_Y3[n] != 0) { |
| k += sns_Y3[n]; |
| corr_xy += sns_Y3[n] * abs_x[n]; |
| energy_y += sns_Y3[n] * sns_Y3[n]; |
| } |
| } |
| |
| // step 2 -> y3 -> y2,start |
| // Add unit pulses until you reach K=6 over N=16 samples, save y3 |
| double corr_xy_last = corr_xy; |
| double energy_y_last = energy_y; |
| for (/*start with existing k intentionally*/; k < K; k++) { |
| uint8_t n_c = 0; |
| uint8_t n_best = 0; |
| corr_xy = corr_xy_last + 1 * abs_x[n_c]; |
| double bestCorrSq = corr_xy * corr_xy; |
| double bestEn = energy_y_last + 2 * 1 * sns_Y3[n_c] + 1 * 1; |
| for (n_c = 1; n_c < N; n_c++) { |
| corr_xy = corr_xy_last + 1 * abs_x[n_c]; |
| energy_y = energy_y_last + 2 * 1 * sns_Y3[n_c] + 1 * 1; |
| if (corr_xy * corr_xy * bestEn > bestCorrSq * energy_y) { |
| n_best = n_c; |
| bestCorrSq = corr_xy * corr_xy; |
| bestEn = energy_y; |
| } |
| } |
| corr_xy_last += 1 * abs_x[n_best]; |
| energy_y_last += 2 * 1 * sns_Y3[n_best] + 1 * 1; |
| sns_Y3[n_best]++; |
| } |
| |
| // (2) Shape name: βoutlier_nearβ |
| // Scale factor set A: {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15} |
| // Scale factor set B: Empty set |
| // Pulse configuration, Set A, PVQ(NA, KA): PVQ(16, 8) |
| // Pulse configuration, Set B, PVQ(NB, KB): Empty |
| |
| // step 3 -> y2 -> y1,pre-start |
| // Add unit pulses until you reach K=8 over N=16 samples, save y2 |
| N = 16; |
| K = 8; |
| for (uint8_t n = 0; n < N; n++) { |
| sns_Y2[n] = sns_Y3[n]; |
| } |
| for (/*start with existing k intentionally*/; k < K; k++) { |
| uint8_t n_c = 0; |
| uint8_t n_best = 0; |
| corr_xy = corr_xy_last + 1 * abs_x[n_c]; |
| double bestCorrSq = corr_xy * corr_xy; |
| double bestEn = energy_y_last + 2 * 1 * sns_Y2[n_c] + 1 * 1; |
| for (n_c = 1; n_c < N; n_c++) { |
| corr_xy = corr_xy_last + 1 * abs_x[n_c]; |
| energy_y = energy_y_last + 2 * 1 * sns_Y2[n_c] + 1 * 1; |
| if (corr_xy * corr_xy * bestEn > bestCorrSq * energy_y) { |
| n_best = n_c; |
| bestCorrSq = corr_xy * corr_xy; |
| bestEn = energy_y; |
| } |
| } |
| corr_xy_last += 1 * abs_x[n_best]; |
| energy_y_last += 2 * 1 * sns_Y2[n_best] + 1 * 1; |
| sns_Y2[n_best]++; |
| } |
| |
| // (1) Shape name: βregular_lfβ |
| // Scale factor set A: {0,1,2,3,4,5,6,7,8,9} |
| // Scale factor set B: {10,11,12,13,14,15} |
| // Pulse configuration, Set A, PVQ(NA, KA): PVQ(10, 10) |
| // Pulse configuration, Set B, PVQ(NB, KB): Zeroed |
| |
| // step 4 -> y1,start |
| // Remove any unit pulses in y1,pre-start that are not part of set A to |
| // yield y1, start |
| N = 10; |
| K = 10; |
| for (uint8_t n = 0; n < N; n++) { |
| sns_Y1[n] = sns_Y2[n]; |
| } |
| |
| // step 5 |
| // Update energy energyy and correlation corrxy terms to reflect the |
| // pulses present in y1, start |
| for (uint8_t n = N; n < 16; n++) { |
| if (sns_Y2[n] != 0) { |
| k -= sns_Y2[n]; |
| corr_xy_last -= sns_Y2[n] * abs_x[n]; |
| energy_y_last -= sns_Y2[n] * sns_Y2[n]; |
| } |
| } |
| |
| // step 6 -> y1 -> y0,start |
| // Add unit pulses until you reach K=10 over N=10 samples (in set A), |
| // save y1 |
| for (/*start with existing k intentionally*/; k < K; k++) { |
| uint8_t n_c = 0; |
| uint8_t n_best = 0; |
| corr_xy = corr_xy_last + 1 * abs_x[n_c]; |
| double bestCorrSq = corr_xy * corr_xy; |
| double bestEn = energy_y_last + 2 * 1 * sns_Y1[n_c] + 1 * 1; |
| // std::cout << static_cast<uint16_t>(k) << " n_c: " << |
| // static_cast<uint16_t>(n_c) << " -> " << (bestCorrSq/bestEn) << " x: " << |
| // abs_x[n_c] << std::endl; |
| for (n_c = 1; n_c < N; n_c++) { |
| corr_xy = corr_xy_last + 1 * abs_x[n_c]; |
| energy_y = energy_y_last + 2 * 1 * sns_Y1[n_c] + 1 * 1; |
| // std::cout << static_cast<uint16_t>(k) << " n_c: " << |
| // static_cast<uint16_t>(n_c) << " -> " << (corr_xy*corr_xy/energy_y) << |
| // " x: " << abs_x[n_c] << std::endl; |
| if (corr_xy * corr_xy * bestEn > bestCorrSq * energy_y) { |
| n_best = n_c; |
| bestCorrSq = corr_xy * corr_xy; |
| bestEn = energy_y; |
| } |
| } |
| corr_xy_last += 1 * abs_x[n_best]; |
| energy_y_last += 2 * 1 * sns_Y1[n_best] + 1 * 1; |
| sns_Y1[n_best]++; |
| // std::cout << static_cast<uint16_t>(k) << " n_best: " << |
| // static_cast<uint16_t>(n_best) << " -> " << |
| // (corr_xy_last*corr_xy_last/energy_y_last) << std::endl; |
| } |
| |
| // (0) Shape name: βregularβ |
| // Scale factor set A: {0,1,2,3,4,5,6,7,8,9} |
| // Scale factor set B: {10,11,12,13,14,15} |
| // Pulse configuration, Set A, PVQ(NA, KA): PVQ(10, 10) |
| // Pulse configuration, Set B, PVQ(NB, KB): PVQ(6, 1) |
| |
| // step 7 -> y0 |
| // Add unit pulses to y0,start until you reach K=1 over N=6 samples |
| // (in set B), save y0 |
| // A |
| N = 10; |
| K = 10; |
| for (uint8_t n = 0; n < N; n++) { |
| sns_Y0[n] = sns_Y1[n]; |
| } |
| // B |
| N = 6; |
| K = 1; |
| double max_abs_x = 0; |
| uint8_t n_best = 0; |
| for (uint8_t n_c = 10; n_c < N + 10; n_c++) { |
| sns_Y0[n_c] = 0; |
| if (abs_x[n_c] > max_abs_x) { |
| max_abs_x = abs_x[n_c]; |
| n_best = n_c; |
| } |
| } |
| sns_Y0[n_best] = 1; |
| |
| // step 8 -> y3, y2, y1, y0 |
| // Add signs to non-zero positions of each yj vector from the target |
| // vector x, save y3, y2, y1, y0 as shape vector candidates (and for |
| // subsequent indexing of one of them) |
| for (uint8_t n = 0; n < 10; n++) { |
| if (t2rot[n] < 0) { |
| sns_Y0[n] *= -1; |
| sns_Y1[n] *= -1; |
| sns_Y2[n] *= -1; |
| sns_Y3[n] *= -1; |
| } |
| } |
| for (uint8_t n = 10; n < 16; n++) { |
| if (t2rot[n] < 0) { |
| sns_Y0[n] *= -1; |
| sns_Y2[n] *= -1; |
| sns_Y3[n] *= -1; |
| } |
| } |
| |
| // step 9 |
| // Unit energy normalize each yj vector to candidate vector xq_j |
| normalizeCandidate(sns_Y0, sns_XQ0, 16); |
| normalizeCandidate(sns_Y1, sns_XQ1, 10); |
| normalizeCandidate(sns_Y2, sns_XQ2, 16); |
| normalizeCandidate(sns_Y3, sns_XQ3, 16); |
| |
| // 3.3.7.3.3.6 Adjustment gain candidates (d09r02_F2F) |
| // 3.3.7.3.3.7 Shape and gain combination determination (d09r02_F2F) |
| shape_j = 0; |
| Gind = 0; |
| double G_gain_i_shape_j = 0; |
| double* sns_XQ_shape_j = nullptr; |
| double dMSE_min = std::numeric_limits<double>::infinity(); |
| for (uint8_t j = 0; j < 4; j++) { |
| uint8_t Gmaxind_j; |
| double* sns_vq_gains; |
| double* sns_XQ; |
| switch (j) { |
| default:; // intentionally fall through to case 0; just to avoid compiler |
| // warning |
| case 0: { |
| Gmaxind_j = 1; |
| sns_vq_gains = sns_vq_reg_adj_gains; |
| sns_XQ = sns_XQ0; |
| break; |
| } |
| case 1: { |
| Gmaxind_j = 3; |
| sns_vq_gains = sns_vq_reg_lf_adj_gains; |
| sns_XQ = sns_XQ1; |
| break; |
| } |
| case 2: { |
| Gmaxind_j = 3; |
| sns_vq_gains = sns_vq_near_adj_gains; |
| sns_XQ = sns_XQ2; |
| break; |
| } |
| case 3: { |
| Gmaxind_j = 7; |
| sns_vq_gains = sns_vq_far_adj_gains; |
| sns_XQ = sns_XQ3; |
| break; |
| } |
| } |
| for (uint8_t i = 0; i <= Gmaxind_j; i++) { |
| double dMSE = 0; |
| for (uint8_t n = 0; n < Nscales; n++) { |
| double diff = t2rot[n] - sns_vq_gains[i] * sns_XQ[n]; |
| dMSE += diff * diff; |
| } |
| if (dMSE < dMSE_min) { |
| shape_j = j; |
| Gind = i; |
| dMSE_min = dMSE; |
| G_gain_i_shape_j = sns_vq_gains[Gind]; |
| sns_XQ_shape_j = sns_XQ; |
| } |
| } |
| } |
| uint8_t LSB_gain = Gind & 1; |
| |
| // 3.3.7.3.3.8 Enumeration of the selected PVQ pulse configurations |
| // (d09r02_F2F) |
| switch (shape_j) { |
| case 0: { |
| MPVQenum(idxA, LS_indA, 10, /* i : dimension of vec_in */ |
| sns_Y0 /* i : PVQ integer pulse train */ |
| ); |
| MPVQenum(idxB, LS_indB, 6, /* i : dimension of vec_in */ |
| &sns_Y0[10] /* i : PVQ integer pulse train */ |
| ); |
| const uint32_t SZ_shapeA_0 = 2390004; |
| index_joint_j = (2 * idxB + LS_indB + 2) * SZ_shapeA_0 + idxA; |
| break; |
| } |
| case 1: { |
| MPVQenum(idxA, LS_indA, 10, /* i : dimension of vec_in */ |
| sns_Y1 /* i : PVQ integer pulse train */ |
| ); |
| const uint32_t SZ_shapeA_1 = 2390004; |
| index_joint_j = LSB_gain * SZ_shapeA_1 + idxA; |
| break; |
| } |
| case 2: { |
| MPVQenum(idxA, LS_indA, 16, /* i : dimension of vec_in */ |
| sns_Y2 /* i : PVQ integer pulse train */ |
| ); |
| index_joint_j = idxA; |
| break; |
| } |
| case 3: { |
| MPVQenum(idxA, LS_indA, 16, /* i : dimension of vec_in */ |
| sns_Y3 /* i : PVQ integer pulse train */ |
| ); |
| const uint32_t SZ_shapeA_2 = 15158272; |
| index_joint_j = SZ_shapeA_2 + LSB_gain + (2 * idxA); |
| break; |
| } |
| } |
| |
| // 3.3.7.3.4 Multiplexing of SNS VQ codewords (d09r02_F2F) |
| // *** final multiplexing skipped for now **** |
| // computation of index_joint_j already done above |
| |
| // 5.3.7.3.5 Synthesis of the Quantized SNS scale factor vector (d09r01) |
| // [3.3.7.3.4 got lost] Synthesis of the Quantized SNS scale factor vector |
| // (d09r02_F2F) |
| for (uint8_t n = 0; n < Nscales; n++) { |
| double factor = 0; |
| for (uint8_t col = 0; col < 16; col++) { |
| factor += sns_XQ_shape_j[col] * D[n][col]; |
| } |
| scfQ[n] = st1[n] + G_gain_i_shape_j * factor; |
| } |
| } |
| |
| void SnsQuantization::normalizeCandidate(int8_t* y, double* XQ, uint8_t N) { |
| double normValue = 0; |
| for (uint8_t n = 0; n < N; n++) { |
| if (y[n] != 0) // try to save some computations |
| { |
| normValue += y[n] * y[n]; |
| } |
| } |
| normValue = sqrt(normValue); |
| for (uint8_t n = 0; n < N; n++) { |
| XQ[n] = y[n]; |
| if (y[n] != 0) // try to save some computations |
| { |
| XQ[n] /= normValue; |
| } |
| } |
| // ensure trailing zero values for N<Nscales (e.g. for XQ1) |
| for (uint8_t n = N; n < Nscales; n++) { |
| XQ[n] = 0; |
| } |
| } |
| |
| void SnsQuantization::MPVQenum(int32_t& index, int32_t& lead_sign_ind, |
| uint8_t dim_in, /* i : dimension of vec_in */ |
| int8_t vec_in[] /* i : PVQ integer pulse train */ |
| ) { |
| /* init */ |
| int32_t next_sign_ind = 0x80000000U; /* sentinel for first sign */ |
| int32_t k_val_acc = 0; |
| int8_t pos = dim_in; |
| index = 0; |
| uint8_t n = 0; |
| unsigned int* row_ptr = &(MPVQ_offsets[n][0]); |
| /* MPVQ-index composition loop */ |
| unsigned int tmp_h_row = row_ptr[0]; |
| for (pos--; pos >= 0; pos--) { |
| int8_t tmp_val = vec_in[pos]; |
| |
| //[index, next_sign_ind] = encPushSign(tmp_val, next_sign_ind, index); |
| index = encPushSign(tmp_val, next_sign_ind, index); |
| index += tmp_h_row; |
| k_val_acc += (tmp_val < 0) ? -tmp_val : tmp_val; |
| if (pos != 0) { |
| n += 1; /* switch row in offset table MPVQ_offsets(n, k) */ |
| } |
| row_ptr = &(MPVQ_offsets[n][0]); |
| tmp_h_row = row_ptr[k_val_acc]; |
| } |
| lead_sign_ind = next_sign_ind; |
| // return [ index, lead_sign_ind ] ; |
| } |
| |
| //[ index, next_sign_ind ] = |
| // encPushSign( val, next_sign_ind_in, index_in) |
| uint32_t SnsQuantization::encPushSign(int8_t val, int32_t& next_sign_ind_in, |
| int32_t index_in) { |
| int32_t index = index_in; |
| if (((next_sign_ind_in & 0x80000000U) == 0) && (val != 0)) { |
| index = 2 * index_in + next_sign_ind_in; |
| } |
| // next_sign_ind = next_sign_ind_in; |
| if (val < 0) { |
| next_sign_ind_in = 1; |
| } |
| if (val > 0) { |
| next_sign_ind_in = 0; |
| } |
| /* if val==0, there is no new sign information to βpushβ, |
| i.e. next_sign_ind is not changed */ |
| // return [ index, next_sign_ind ]; |
| return index; |
| } |
| |
| void SnsQuantization::registerDatapoints(DatapointContainer* datapoints_) { |
| datapoints = datapoints_; |
| if (nullptr != datapoints) { |
| datapoints->addDatapoint("ind_LF", &ind_LF, sizeof(ind_LF)); |
| datapoints->addDatapoint("ind_HF", &ind_HF, sizeof(ind_HF)); |
| datapoints->addDatapoint("shape_j", &shape_j, sizeof(shape_j)); |
| datapoints->addDatapoint("Gind", &Gind, sizeof(Gind)); |
| datapoints->addDatapoint("LS_indA", &LS_indA, sizeof(LS_indA)); |
| datapoints->addDatapoint("LS_indB", &LS_indB, sizeof(LS_indB)); |
| datapoints->addDatapoint("idxA", &idxA, sizeof(idxA)); |
| datapoints->addDatapoint("idxB", &idxB, sizeof(idxB)); |
| |
| datapoints->addDatapoint("t2rot", &t2rot[0], sizeof(double) * Nscales); |
| datapoints->addDatapoint("scfQ", &scfQ[0], sizeof(double) * Nscales); |
| |
| datapoints->addDatapoint("sns_Y0", &sns_Y0[0], sizeof(int8_t) * 16); |
| datapoints->addDatapoint("sns_Y1", &sns_Y1[0], sizeof(int8_t) * 10); |
| datapoints->addDatapoint("sns_Y2", &sns_Y2[0], sizeof(int8_t) * 16); |
| datapoints->addDatapoint("sns_Y3", &sns_Y3[0], sizeof(int8_t) * 16); |
| |
| datapoints->addDatapoint("sns_XQ0", &sns_XQ0[0], sizeof(double) * 16); |
| datapoints->addDatapoint("sns_XQ1", &sns_XQ1[0], sizeof(double) * 10); |
| datapoints->addDatapoint("sns_XQ2", &sns_XQ2[0], sizeof(double) * 16); |
| datapoints->addDatapoint("sns_XQ3", &sns_XQ3[0], sizeof(double) * 16); |
| } |
| } |
| |
| } // namespace Lc3Enc |