1 /******************************************************************************
2 *
3 * Copyright 2022 Google LLC
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
19 #include "tns.h"
20 #include "tables.h"
21
22
23 /* ----------------------------------------------------------------------------
24 * Filter Coefficients
25 * -------------------------------------------------------------------------- */
26
27 /**
28 * Resolve LPC Weighting indication according bitrate
29 * dt, nbytes Duration and size of the frame
30 * return True when LPC Weighting enabled
31 */
resolve_lpc_weighting(enum lc3_dt dt,int nbytes)32 static bool resolve_lpc_weighting(enum lc3_dt dt, int nbytes)
33 {
34 return nbytes * 8 < 120 * (int)(1 + dt);
35 }
36
37 /**
38 * Return dot product of 2 vectors
39 * a, b, n The 2 vectors of size `n`
40 * return sum( a[i] * b[i] ), i = [0..n-1]
41 */
dot(const float * a,const float * b,int n)42 LC3_HOT static inline float dot(const float *a, const float *b, int n)
43 {
44 float v = 0;
45
46 while (n--)
47 v += *(a++) * *(b++);
48
49 return v;
50 }
51
52 /**
53 * LPC Coefficients
54 * dt, bw Duration and bandwidth of the frame
55 * maxorder Maximum order of filter
56 * x Spectral coefficients
57 * gain, a Output the prediction gains and LPC coefficients
58 */
compute_lpc_coeffs(enum lc3_dt dt,enum lc3_bandwidth bw,int maxorder,const float * x,float * gain,float (* a)[9])59 LC3_HOT static void compute_lpc_coeffs(
60 enum lc3_dt dt, enum lc3_bandwidth bw, int maxorder,
61 const float *x, float *gain, float (*a)[9])
62 {
63
64 #if LC3_PLUS
65
66 static const int sub_2m5_nb[] = { 3, 10, 20 };
67 static const int sub_2m5_wb[] = { 3, 20, 40 };
68 static const int sub_2m5_sswb[] = { 3, 30, 60 };
69 static const int sub_2m5_swb[] = { 3, 40, 80 };
70 static const int sub_2m5_fb[] = { 3, 51, 100 };
71
72 static const int sub_5m_nb[] = { 6, 23, 40 };
73 static const int sub_5m_wb[] = { 6, 43, 80 };
74 static const int sub_5m_sswb[] = { 6, 63, 120 };
75 static const int sub_5m_swb[] = { 6, 43, 80, 120, 160 };
76 static const int sub_5m_fb[] = { 6, 53, 100, 150, 200 };
77
78 #endif /* LC3_PLUS */
79
80 static const int sub_7m5_nb[] = { 9, 26, 43, 60 };
81 static const int sub_7m5_wb[] = { 9, 46, 83, 120 };
82 static const int sub_7m5_sswb[] = { 9, 66, 123, 180 };
83 static const int sub_7m5_swb[] = { 9, 46, 82, 120, 159, 200, 240 };
84 static const int sub_7m5_fb[] = { 9, 56, 103, 150, 200, 250, 300 };
85
86 static const int sub_10m_nb[] = { 12, 34, 57, 80 };
87 static const int sub_10m_wb[] = { 12, 61, 110, 160 };
88 static const int sub_10m_sswb[] = { 12, 88, 164, 240 };
89 static const int sub_10m_swb[] = { 12, 61, 110, 160, 213, 266, 320 };
90 static const int sub_10m_fb[] = { 12, 74, 137, 200, 266, 333, 400 };
91
92 static const float lag_window[] = {
93 1.00000000e+00, 9.98028026e-01, 9.92135406e-01, 9.82391584e-01,
94 9.68910791e-01, 9.51849807e-01, 9.31404933e-01, 9.07808230e-01,
95 8.81323137e-01
96 };
97
98 const int *sub = (const int * const [LC3_NUM_DT][LC3_NUM_BANDWIDTH]){
99
100 #if LC3_PLUS
101
102 [LC3_DT_2M5] = {
103 sub_2m5_nb, sub_2m5_wb, sub_2m5_sswb, sub_2m5_swb,
104 sub_2m5_fb, sub_2m5_fb, sub_2m5_fb },
105
106 [LC3_DT_5M] = {
107 sub_5m_nb , sub_5m_wb , sub_5m_sswb , sub_5m_swb ,
108 sub_5m_fb , sub_5m_fb , sub_5m_fb },
109
110 #endif /* LC3_PLUS */
111
112 [LC3_DT_7M5] = {
113 sub_7m5_nb, sub_7m5_wb, sub_7m5_sswb, sub_7m5_swb,
114 sub_7m5_fb },
115
116 [LC3_DT_10M] = {
117 sub_10m_nb, sub_10m_wb, sub_10m_sswb, sub_10m_swb,
118 sub_10m_fb, sub_10m_fb, sub_10m_fb },
119
120 }[dt][bw];
121
122 /* --- Normalized autocorrelation --- */
123
124 int nfilters = 1 + (dt >= LC3_DT_5M && bw >= LC3_BANDWIDTH_SWB);
125 int nsubdivisions = 2 + (dt >= LC3_DT_7M5);
126
127 const float *xs, *xe = x + *sub;
128 float r[2][9];
129
130 for (int f = 0; f < nfilters; f++) {
131 float c[9][3] = { 0 };
132
133 for (int s = 0; s < nsubdivisions; s++) {
134 xs = xe, xe = x + *(++sub);
135
136 for (int k = 0; k <= maxorder; k++)
137 c[k][s] = dot(xs, xs + k, (xe - xs) - k);
138 }
139
140 r[f][0] = nsubdivisions;
141 if (nsubdivisions == 2) {
142 float e0 = c[0][0], e1 = c[0][1];
143 for (int k = 1; k <= maxorder; k++)
144 r[f][k] = e0 == 0 || e1 == 0 ? 0 :
145 (c[k][0]/e0 + c[k][1]/e1) * lag_window[k];
146
147 } else {
148 float e0 = c[0][0], e1 = c[0][1], e2 = c[0][2];
149 for (int k = 1; k <= maxorder; k++)
150 r[f][k] = e0 == 0 || e1 == 0 || e2 == 0 ? 0 :
151 (c[k][0]/e0 + c[k][1]/e1 + c[k][2]/e2) * lag_window[k];
152 }
153 }
154
155 /* --- Levinson-Durbin recursion --- */
156
157 for (int f = 0; f < nfilters; f++) {
158 float *a0 = a[f], a1[9];
159 float err = r[f][0], rc;
160
161 gain[f] = err;
162
163 a0[0] = 1;
164 for (int k = 1; k <= maxorder; ) {
165
166 rc = -r[f][k];
167 for (int i = 1; i < k; i++)
168 rc -= a0[i] * r[f][k-i];
169
170 rc /= err;
171 err *= 1 - rc * rc;
172
173 for (int i = 1; i < k; i++)
174 a1[i] = a0[i] + rc * a0[k-i];
175 a1[k++] = rc;
176
177 rc = -r[f][k];
178 for (int i = 1; i < k; i++)
179 rc -= a1[i] * r[f][k-i];
180
181 rc /= err;
182 err *= 1 - rc * rc;
183
184 for (int i = 1; i < k; i++)
185 a0[i] = a1[i] + rc * a1[k-i];
186 a0[k++] = rc;
187 }
188
189 gain[f] /= err;
190 }
191 }
192
193 /**
194 * LPC Weighting
195 * gain, a Prediction gain and LPC coefficients, weighted as output
196 */
lpc_weighting(float pred_gain,float * a)197 LC3_HOT static void lpc_weighting(float pred_gain, float *a)
198 {
199 float gamma = 1.f - (1.f - 0.85f) * (2.f - pred_gain) / (2.f - 1.5f);
200 float g = 1.f;
201
202 for (int i = 1; i < 9; i++)
203 a[i] *= (g *= gamma);
204 }
205
206 /**
207 * LPC reflection
208 * a, maxorder LPC coefficients, and maximum order (4 or 8)
209 * rc Output refelection coefficients
210 */
lpc_reflection(const float * a,int maxorder,float * rc)211 LC3_HOT static void lpc_reflection(
212 const float *a, int maxorder, float *rc)
213 {
214 float e, b[2][7], *b0, *b1;
215
216 rc[maxorder-1] = a[maxorder];
217 e = 1 - rc[maxorder-1] * rc[maxorder-1];
218
219 b1 = b[1];
220 for (int i = 0; i < maxorder-1; i++)
221 b1[i] = (a[1+i] - rc[maxorder-1] * a[(maxorder-1)-i]) / e;
222
223 for (int k = maxorder-2; k > 0; k--) {
224 b0 = b1, b1 = b[k & 1];
225
226 rc[k] = b0[k];
227 e = 1 - rc[k] * rc[k];
228
229 for (int i = 0; i < k; i++)
230 b1[i] = (b0[i] - rc[k] * b0[k-1-i]) / e;
231 }
232
233 rc[0] = b1[0];
234 }
235
236 /**
237 * Quantization of RC coefficients
238 * rc, maxorder Refelection coefficients, and maximum order (4 or 8)
239 * order Return order of coefficients
240 * rc_i Return quantized coefficients
241 */
quantize_rc(const float * rc,int maxorder,int * order,int * rc_q)242 static void quantize_rc(const float *rc, int maxorder, int *order, int *rc_q)
243 {
244 /* Quantization table, sin(delta * (i + 0.5)), delta = Pi / 17,
245 * rounded to fixed point Q15 value (LC3-Plus HR requirements). */
246
247 static float q_thr[] = {
248 0x0bcfp-15, 0x2307p-15, 0x390ep-15, 0x4d23p-15,
249 0x5e98p-15, 0x6cd4p-15, 0x775bp-15, 0x7dd2p-15,
250 };
251
252 *order = maxorder;
253
254 for (int i = 0; i < maxorder; i++) {
255 float rc_m = fabsf(rc[i]);
256
257 rc_q[i] = 4 * (rc_m >= q_thr[4]);
258 for (int j = 0; j < 4 && rc_m >= q_thr[rc_q[i]]; j++, rc_q[i]++);
259
260 if (rc[i] < 0)
261 rc_q[i] = -rc_q[i];
262
263 *order = rc_q[i] != 0 ? maxorder : *order - 1;
264 }
265 }
266
267 /**
268 * Unquantization of RC coefficients
269 * rc_q, order Quantized coefficients, and order
270 * rc Return refelection coefficients
271 */
unquantize_rc(const int * rc_q,int order,float rc[8])272 static void unquantize_rc(const int *rc_q, int order, float rc[8])
273 {
274 /* Quantization table, sin(delta * i), delta = Pi / 17,
275 * rounded to fixed point Q15 value (LC3-Plus HR requirements). */
276
277 static float q_inv[] = {
278 0x0000p-15, 0x1785p-15, 0x2e3dp-15, 0x4362p-15,
279 0x563cp-15, 0x6625p-15, 0x7295p-15, 0x7b1dp-15, 0x7f74p-15,
280 };
281
282 int i;
283
284 for (i = 0; i < order; i++) {
285 float rc_m = q_inv[LC3_ABS(rc_q[i])];
286 rc[i] = rc_q[i] < 0 ? -rc_m : rc_m;
287 }
288 }
289
290
291 /* ----------------------------------------------------------------------------
292 * Filtering
293 * -------------------------------------------------------------------------- */
294
295 /**
296 * Forward filtering
297 * dt, bw Duration and bandwidth of the frame
298 * rc_order, rc Order of coefficients, and coefficients
299 * x Spectral coefficients, filtered as output
300 */
forward_filtering(enum lc3_dt dt,enum lc3_bandwidth bw,const int rc_order[2],float (* const rc)[8],float * x)301 LC3_HOT static void forward_filtering(
302 enum lc3_dt dt, enum lc3_bandwidth bw,
303 const int rc_order[2], float (* const rc)[8], float *x)
304 {
305 int nfilters = 1 + (dt >= LC3_DT_5M && bw >= LC3_BANDWIDTH_SWB);
306 int nf = lc3_ne(dt, (enum lc3_srate)LC3_MIN(bw, LC3_BANDWIDTH_FB))
307 >> (nfilters - 1);
308 int i0, ie = 3*(1 + dt);
309
310 float s[8] = { 0 };
311
312 for (int f = 0; f < nfilters; f++) {
313
314 i0 = ie;
315 ie = nf * (1 + f);
316
317 if (!rc_order[f])
318 continue;
319
320 for (int i = i0; i < ie; i++) {
321 float xi = x[i];
322 float s0, s1 = xi;
323
324 for (int k = 0; k < rc_order[f]; k++) {
325 s0 = s[k];
326 s[k] = s1;
327
328 s1 = rc[f][k] * xi + s0;
329 xi += rc[f][k] * s0;
330 }
331
332 x[i] = xi;
333 }
334 }
335 }
336
337 /**
338 * Inverse filtering
339 * dt, bw Duration and bandwidth of the frame
340 * rc_order, rc Order of coefficients, and unquantized coefficients
341 * x Spectral coefficients, filtered as output
342 */
inverse_filtering(enum lc3_dt dt,enum lc3_bandwidth bw,const int rc_order[2],float (* const rc)[8],float * x)343 LC3_HOT static void inverse_filtering(
344 enum lc3_dt dt, enum lc3_bandwidth bw,
345 const int rc_order[2], float (* const rc)[8], float *x)
346 {
347 int nfilters = 1 + (dt >= LC3_DT_5M && bw >= LC3_BANDWIDTH_SWB);
348 int nf = lc3_ne(dt, (enum lc3_srate)LC3_MIN(bw, LC3_BANDWIDTH_FB))
349 >> (nfilters - 1);
350 int i0, ie = 3*(1 + dt);
351
352 float s[8] = { 0 };
353
354 for (int f = 0; f < nfilters; f++) {
355
356 i0 = ie;
357 ie = nf * (1 + f);
358
359 if (!rc_order[f])
360 continue;
361
362 for (int i = i0; i < ie; i++) {
363 float xi = x[i];
364
365 xi -= s[7] * rc[f][7];
366 for (int k = 6; k >= 0; k--) {
367 xi -= s[k] * rc[f][k];
368 s[k+1] = s[k] + rc[f][k] * xi;
369 }
370 s[0] = xi;
371 x[i] = xi;
372 }
373
374 for (int k = 7; k >= rc_order[f]; k--)
375 s[k] = 0;
376 }
377 }
378
379
380 /* ----------------------------------------------------------------------------
381 * Interface
382 * -------------------------------------------------------------------------- */
383
384 /**
385 * TNS analysis
386 */
lc3_tns_analyze(enum lc3_dt dt,enum lc3_bandwidth bw,bool nn_flag,int nbytes,struct lc3_tns_data * data,float * x)387 void lc3_tns_analyze(enum lc3_dt dt, enum lc3_bandwidth bw,
388 bool nn_flag, int nbytes, struct lc3_tns_data *data, float *x)
389 {
390 /* Processing steps :
391 * - Determine the LPC (Linear Predictive Coding) Coefficients
392 * - Check is the filtering is disabled
393 * - The coefficients are weighted on low bitrates and predicition gain
394 * - Convert to reflection coefficients and quantize
395 * - Finally filter the spectral coefficients */
396
397 float pred_gain[2], a[2][9];
398 float rc[2][8];
399
400 data->lpc_weighting = resolve_lpc_weighting(dt, nbytes);
401 data->nfilters = 1 + (dt >= LC3_DT_5M && bw >= LC3_BANDWIDTH_SWB);
402 int maxorder = dt <= LC3_DT_5M ? 4 : 8;
403
404 compute_lpc_coeffs(dt, bw, maxorder, x, pred_gain, a);
405
406 for (int f = 0; f < data->nfilters; f++) {
407
408 data->rc_order[f] = 0;
409 if (nn_flag || pred_gain[f] <= 1.5f)
410 continue;
411
412 if (data->lpc_weighting && pred_gain[f] < 2.f)
413 lpc_weighting(pred_gain[f], a[f]);
414
415 lpc_reflection(a[f], maxorder, rc[f]);
416
417 quantize_rc(rc[f], maxorder, &data->rc_order[f], data->rc[f]);
418 unquantize_rc(data->rc[f], data->rc_order[f], rc[f]);
419 }
420
421 forward_filtering(dt, bw, data->rc_order, rc, x);
422 }
423
424 /**
425 * TNS synthesis
426 */
lc3_tns_synthesize(enum lc3_dt dt,enum lc3_bandwidth bw,const struct lc3_tns_data * data,float * x)427 void lc3_tns_synthesize(enum lc3_dt dt, enum lc3_bandwidth bw,
428 const struct lc3_tns_data *data, float *x)
429 {
430 float rc[2][8] = { 0 };
431
432 for (int f = 0; f < data->nfilters; f++)
433 if (data->rc_order[f])
434 unquantize_rc(data->rc[f], data->rc_order[f], rc[f]);
435
436 inverse_filtering(dt, bw, data->rc_order, rc, x);
437 }
438
439 /**
440 * Bit consumption of bitstream data
441 */
lc3_tns_get_nbits(const struct lc3_tns_data * data)442 int lc3_tns_get_nbits(const struct lc3_tns_data *data)
443 {
444 int nbits = 0;
445
446 for (int f = 0; f < data->nfilters; f++) {
447
448 int nbits_2048 = 2048;
449 int rc_order = data->rc_order[f];
450
451 nbits_2048 += rc_order > 0 ? lc3_tns_order_bits
452 [data->lpc_weighting][rc_order-1] : 0;
453
454 for (int i = 0; i < rc_order; i++)
455 nbits_2048 += lc3_tns_coeffs_bits[i][8 + data->rc[f][i]];
456
457 nbits += (nbits_2048 + (1 << 11) - 1) >> 11;
458 }
459
460 return nbits;
461 }
462
463 /**
464 * Put bitstream data
465 */
lc3_tns_put_data(lc3_bits_t * bits,const struct lc3_tns_data * data)466 void lc3_tns_put_data(lc3_bits_t *bits, const struct lc3_tns_data *data)
467 {
468 for (int f = 0; f < data->nfilters; f++) {
469 int rc_order = data->rc_order[f];
470
471 lc3_put_bits(bits, rc_order > 0, 1);
472 if (rc_order <= 0)
473 continue;
474
475 lc3_put_symbol(bits,
476 lc3_tns_order_models + data->lpc_weighting, rc_order-1);
477
478 for (int i = 0; i < rc_order; i++)
479 lc3_put_symbol(bits,
480 lc3_tns_coeffs_models + i, 8 + data->rc[f][i]);
481 }
482 }
483
484 /**
485 * Get bitstream data
486 */
lc3_tns_get_data(lc3_bits_t * bits,enum lc3_dt dt,enum lc3_bandwidth bw,int nbytes,lc3_tns_data_t * data)487 int lc3_tns_get_data(lc3_bits_t *bits,
488 enum lc3_dt dt, enum lc3_bandwidth bw, int nbytes, lc3_tns_data_t *data)
489 {
490 data->nfilters = 1 + (dt >= LC3_DT_5M && bw >= LC3_BANDWIDTH_SWB);
491 data->lpc_weighting = resolve_lpc_weighting(dt, nbytes);
492
493 for (int f = 0; f < data->nfilters; f++) {
494
495 data->rc_order[f] = lc3_get_bit(bits);
496 if (!data->rc_order[f])
497 continue;
498
499 data->rc_order[f] += lc3_get_symbol(bits,
500 lc3_tns_order_models + data->lpc_weighting);
501 if (dt <= LC3_DT_5M && data->rc_order[f] > 4)
502 return -1;
503
504 for (int i = 0; i < data->rc_order[f]; i++)
505 data->rc[f][i] = (int)lc3_get_symbol(bits,
506 lc3_tns_coeffs_models + i) - 8;
507 }
508
509 return 0;
510 }
511
