1 /******************************************************************************
2 * *
3 * Copyright (C) 2023 The Android Open Source Project
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at:
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *****************************************************************************
18 * Originally developed and contributed by Ittiam Systems Pvt. Ltd, Bangalore
19 */
20
21 #include <math.h>
22 #include "ixheaac_type_def.h"
23 #include "iusace_bitbuffer.h"
24 #include "iusace_cnst.h"
25 #include "iusace_tns_usac.h"
26 #include "iusace_psy_mod.h"
27 #include "iusace_block_switch_const.h"
28 #include "iusace_rom.h"
29
iusace_lpc_eval_chebyshev_polyn(FLOAT32 x,FLOAT32 * coefs,WORD32 order)30 static FLOAT32 iusace_lpc_eval_chebyshev_polyn(FLOAT32 x, FLOAT32 *coefs, WORD32 order) {
31 WORD32 i;
32 FLOAT32 b0, b1, b2, x2;
33 x2 = 2.0f * x;
34 b2 = 1.0f;
35 b1 = x2 + coefs[1];
36 for (i = 2; i < order; i++) {
37 b0 = x2 * b1 - b2 + coefs[i];
38 b2 = b1;
39 b1 = b0;
40 }
41 return (x * b1 - b2 + 0.5f * coefs[order]);
42 }
43
iusace_lpc_2_lsp_conversion(FLOAT32 * lpc,FLOAT32 * lsp,FLOAT32 * prev_lsp)44 VOID iusace_lpc_2_lsp_conversion(FLOAT32 *lpc, FLOAT32 *lsp, FLOAT32 *prev_lsp) {
45 FLOAT32 sum_polyn[(ORDER_BY_2) + 1], diff_polyn[(ORDER_BY_2) + 1];
46 FLOAT32 *p1_lpc, *p2_lpc, *p_sum_polyn, *p_diff_polyn;
47 WORD32 i, j = 0, num_found_freeq = 0, is_first_polyn = 0;
48 FLOAT32 x_low, y_low, x_high, y_high, x_mid, y_mid, x_lin_interp;
49
50 p_sum_polyn = sum_polyn;
51 p_diff_polyn = diff_polyn;
52 *p_sum_polyn++ = 1.0f;
53 *p_diff_polyn++ = 1.0f;
54 sum_polyn[0] = 1.0f;
55 diff_polyn[0] = 1.0f;
56 p1_lpc = lpc + 1;
57 p2_lpc = lpc + ORDER;
58 for (i = 0; i <= ORDER_BY_2 - 1; i++) {
59 *p_sum_polyn = *p1_lpc + *p2_lpc - *(p_sum_polyn - 1);
60 p_sum_polyn++;
61 *p_diff_polyn = *p1_lpc++ - *p2_lpc-- + *(p_diff_polyn - 1);
62 p_diff_polyn++;
63 }
64 p_sum_polyn = sum_polyn;
65 x_low = iusace_chebyshev_polyn_grid[0];
66 y_low = iusace_lpc_eval_chebyshev_polyn(x_low, p_sum_polyn, ORDER_BY_2);
67
68 while ((num_found_freeq < ORDER) && (j < CHEBYSHEV_NUM_POINTS)) {
69 j++;
70 x_high = x_low;
71 y_high = y_low;
72 x_low = iusace_chebyshev_polyn_grid[j];
73 y_low = iusace_lpc_eval_chebyshev_polyn(x_low, p_sum_polyn, ORDER_BY_2);
74
75 if (y_low * y_high <= 0.0) /* if sign change new root exists */
76 {
77 j--;
78 for (i = 0; i < CHEBYSHEV_NUM_ITER; i++) {
79 x_mid = 0.5f * (x_low + x_high);
80 y_mid = iusace_lpc_eval_chebyshev_polyn(x_mid, p_sum_polyn, ORDER_BY_2);
81 if (y_low * y_mid <= 0.0) {
82 y_high = y_mid;
83 x_high = x_mid;
84 } else {
85 y_low = y_mid;
86 x_low = x_mid;
87 }
88 }
89
90 /* linear interpolation for evaluating the root */
91 x_lin_interp = x_low - y_low * (x_high - x_low) / (y_high - y_low);
92
93 lsp[num_found_freeq] = x_lin_interp;
94 num_found_freeq++;
95
96 is_first_polyn = 1 - is_first_polyn;
97 p_sum_polyn = is_first_polyn ? diff_polyn : sum_polyn;
98
99 x_low = x_lin_interp;
100 y_low = iusace_lpc_eval_chebyshev_polyn(x_low, p_sum_polyn, ORDER_BY_2);
101 }
102 }
103
104 /* Check if ORDER roots found */
105 /* if not use the LSPs from previous frame */
106 if (num_found_freeq < ORDER) {
107 for (i = 0; i < ORDER; i++) lsp[i] = prev_lsp[i];
108 }
109 }
110
iusace_compute_coeff_poly_f(FLOAT32 * lsp,FLOAT32 * poly1,FLOAT32 * poly2)111 static VOID iusace_compute_coeff_poly_f(FLOAT32 *lsp, FLOAT32 *poly1, FLOAT32 *poly2) {
112 FLOAT32 b1, b2;
113 FLOAT32 *ptr_lsp;
114 WORD32 i, j;
115
116 ptr_lsp = lsp;
117 poly1[0] = poly2[0] = 1.0f;
118
119 for (i = 1; i <= ORDER_BY_2; i++) {
120 b1 = -2.0f * (*ptr_lsp++);
121 b2 = -2.0f * (*ptr_lsp++);
122 poly1[i] = (b1 * poly1[i - 1]) + (2.0f * poly1[i - 2]);
123 poly2[i] = (b2 * poly2[i - 1]) + (2.0f * poly2[i - 2]);
124 for (j = i - 1; j > 0; j--) {
125 poly1[j] += (b1 * poly1[j - 1]) + poly1[j - 2];
126 poly2[j] += (b2 * poly2[j - 1]) + poly2[j - 2];
127 }
128 }
129 }
130
iusace_lsp_to_lp_conversion(FLOAT32 * lsp,FLOAT32 * lp_flt_coff_a)131 VOID iusace_lsp_to_lp_conversion(FLOAT32 *lsp, FLOAT32 *lp_flt_coff_a) {
132 WORD32 i;
133 FLOAT32 *ppoly_f1, *ppoly_f2;
134 FLOAT32 *plp_flt_coff_a_bott, *plp_flt_coff_a_top;
135 FLOAT32 poly1[ORDER_BY_2 + 2], poly2[ORDER_BY_2 + 2];
136
137 poly1[0] = 0.0f;
138 poly2[0] = 0.0f;
139
140 iusace_compute_coeff_poly_f(lsp, &poly1[1], &poly2[1]);
141
142 ppoly_f1 = poly1 + ORDER_BY_2 + 1;
143 ppoly_f2 = poly2 + ORDER_BY_2 + 1;
144
145 for (i = 0; i < ORDER_BY_2; i++) {
146 ppoly_f1[0] += ppoly_f1[-1];
147 ppoly_f2[0] -= ppoly_f2[-1];
148 ppoly_f1--;
149 ppoly_f2--;
150 }
151
152 plp_flt_coff_a_bott = lp_flt_coff_a;
153 *plp_flt_coff_a_bott++ = 1.0f;
154 plp_flt_coff_a_top = lp_flt_coff_a + ORDER;
155 ppoly_f1 = poly1 + 2;
156 ppoly_f2 = poly2 + 2;
157 for (i = 0; i < ORDER_BY_2; i++) {
158 *plp_flt_coff_a_bott++ = 0.5f * (*ppoly_f1 + *ppoly_f2);
159 *plp_flt_coff_a_top-- = 0.5f * (*ppoly_f1++ - *ppoly_f2++);
160 }
161 }
162
iusace_levinson_durbin_algo(FLOAT32 * auto_corr_input,FLOAT32 * lpc)163 VOID iusace_levinson_durbin_algo(FLOAT32 *auto_corr_input, FLOAT32 *lpc) {
164 WORD32 i, j;
165 FLOAT32 lpc_val, sum, sigma;
166 FLOAT32 reflection_coeffs[LEV_DUR_MAX_ORDER];
167
168 lpc[0] = 1.0f;
169
170 reflection_coeffs[0] = -auto_corr_input[1] / auto_corr_input[0];
171 lpc[1] = reflection_coeffs[0];
172 sigma = auto_corr_input[0] + auto_corr_input[1] * reflection_coeffs[0];
173
174 for (i = 2; i <= ORDER; i++) {
175 sum = 0.0f;
176 for (j = 0; j < i; j++) sum += auto_corr_input[i - j] * lpc[j];
177 reflection_coeffs[i - 1] = -sum / sigma;
178
179 sigma = sigma * (1.0f - reflection_coeffs[i - 1] * reflection_coeffs[i - 1]);
180
181 if (sigma <= 1.0E-09f) {
182 for (j = i; j <= ORDER; j++) {
183 reflection_coeffs[j - 1] = 0.0f;
184 lpc[j] = 0.0f;
185 }
186 break;
187 }
188
189 for (j = 1; j <= (i / 2); j++) {
190 lpc_val = lpc[j] + reflection_coeffs[i - 1] * lpc[i - j];
191 lpc[i - j] += reflection_coeffs[i - 1] * lpc[j];
192 lpc[j] = lpc_val;
193 }
194
195 lpc[i] = reflection_coeffs[i - 1];
196 }
197 }
198
iusace_get_weighted_lpc(FLOAT32 * lpc,FLOAT32 * weighted_lpc)199 VOID iusace_get_weighted_lpc(FLOAT32 *lpc, FLOAT32 *weighted_lpc) {
200 WORD32 i;
201 for (i = 0; i <= ORDER; i++) {
202 weighted_lpc[i] = iusace_gamma_table[i] * lpc[i];
203 }
204 }
205
iusace_lsp_2_lsf_conversion(FLOAT32 * lsp,FLOAT32 * lsf)206 VOID iusace_lsp_2_lsf_conversion(FLOAT32 *lsp, FLOAT32 *lsf) {
207 WORD32 i;
208 for (i = 0; i < ORDER; i++) {
209 lsf[i] = (FLOAT32)(acos(lsp[i]) * LSP_2_LSF_SCALE);
210 }
211 }
212
iusace_lsf_2_lsp_conversion(FLOAT32 * lsf,FLOAT32 * lsp)213 VOID iusace_lsf_2_lsp_conversion(FLOAT32 *lsf, FLOAT32 *lsp) {
214 WORD32 i;
215 for (i = 0; i < ORDER; i++) lsp[i] = (FLOAT32)cos((FLOAT64)lsf[i] * (FLOAT64)PI_BY_6400);
216 }
217