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 < (dt == LC3_DT_7M5 ? 360/8 : 480/8);
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 * x Spectral coefficients
56 * gain, a Output the prediction gains and LPC coefficients
57 */
compute_lpc_coeffs(enum lc3_dt dt,enum lc3_bandwidth bw,const float * x,float * gain,float (* a)[9])58 LC3_HOT static void compute_lpc_coeffs(
59 enum lc3_dt dt, enum lc3_bandwidth bw,
60 const float *x, float *gain, float (*a)[9])
61 {
62 static const int sub_7m5_nb[] = { 9, 26, 43, 60 };
63 static const int sub_7m5_wb[] = { 9, 46, 83, 120 };
64 static const int sub_7m5_sswb[] = { 9, 66, 123, 180 };
65 static const int sub_7m5_swb[] = { 9, 46, 82, 120, 159, 200, 240 };
66 static const int sub_7m5_fb[] = { 9, 56, 103, 150, 200, 250, 300 };
67
68 static const int sub_10m_nb[] = { 12, 34, 57, 80 };
69 static const int sub_10m_wb[] = { 12, 61, 110, 160 };
70 static const int sub_10m_sswb[] = { 12, 88, 164, 240 };
71 static const int sub_10m_swb[] = { 12, 61, 110, 160, 213, 266, 320 };
72 static const int sub_10m_fb[] = { 12, 74, 137, 200, 266, 333, 400 };
73
74 /* --- Normalized autocorrelation --- */
75
76 static const float lag_window[] = {
77 1.00000000e+00, 9.98028026e-01, 9.92135406e-01, 9.82391584e-01,
78 9.68910791e-01, 9.51849807e-01, 9.31404933e-01, 9.07808230e-01,
79 8.81323137e-01
80 };
81
82 const int *sub = (const int * const [LC3_NUM_DT][LC3_NUM_SRATE]){
83 { sub_7m5_nb, sub_7m5_wb, sub_7m5_sswb, sub_7m5_swb, sub_7m5_fb },
84 { sub_10m_nb, sub_10m_wb, sub_10m_sswb, sub_10m_swb, sub_10m_fb },
85 }[dt][bw];
86
87 int nfilters = 1 + (bw >= LC3_BANDWIDTH_SWB);
88
89 const float *xs, *xe = x + *sub;
90 float r[2][9];
91
92 for (int f = 0; f < nfilters; f++) {
93 float c[9][3];
94
95 for (int s = 0; s < 3; s++) {
96 xs = xe, xe = x + *(++sub);
97
98 for (int k = 0; k < 9; k++)
99 c[k][s] = dot(xs, xs + k, (xe - xs) - k);
100 }
101
102 float e0 = c[0][0], e1 = c[0][1], e2 = c[0][2];
103
104 r[f][0] = 3;
105 for (int k = 1; k < 9; k++)
106 r[f][k] = e0 == 0 || e1 == 0 || e2 == 0 ? 0 :
107 (c[k][0]/e0 + c[k][1]/e1 + c[k][2]/e2) * lag_window[k];
108 }
109
110 /* --- Levinson-Durbin recursion --- */
111
112 for (int f = 0; f < nfilters; f++) {
113 float *a0 = a[f], a1[9];
114 float err = r[f][0], rc;
115
116 gain[f] = err;
117
118 a0[0] = 1;
119 for (int k = 1; k < 9; ) {
120
121 rc = -r[f][k];
122 for (int i = 1; i < k; i++)
123 rc -= a0[i] * r[f][k-i];
124
125 rc /= err;
126 err *= 1 - rc * rc;
127
128 for (int i = 1; i < k; i++)
129 a1[i] = a0[i] + rc * a0[k-i];
130 a1[k++] = rc;
131
132 rc = -r[f][k];
133 for (int i = 1; i < k; i++)
134 rc -= a1[i] * r[f][k-i];
135
136 rc /= err;
137 err *= 1 - rc * rc;
138
139 for (int i = 1; i < k; i++)
140 a0[i] = a1[i] + rc * a1[k-i];
141 a0[k++] = rc;
142 }
143
144 gain[f] /= err;
145 }
146 }
147
148 /**
149 * LPC Weighting
150 * gain, a Prediction gain and LPC coefficients, weighted as output
151 */
lpc_weighting(float pred_gain,float * a)152 LC3_HOT static void lpc_weighting(float pred_gain, float *a)
153 {
154 float gamma = 1.f - (1.f - 0.85f) * (2.f - pred_gain) / (2.f - 1.5f);
155 float g = 1.f;
156
157 for (int i = 1; i < 9; i++)
158 a[i] *= (g *= gamma);
159 }
160
161 /**
162 * LPC reflection
163 * a LPC coefficients
164 * rc Output refelection coefficients
165 */
lpc_reflection(const float * a,float * rc)166 LC3_HOT static void lpc_reflection(const float *a, float *rc)
167 {
168 float e, b[2][7], *b0, *b1;
169
170 rc[7] = a[1+7];
171 e = 1 - rc[7] * rc[7];
172
173 b1 = b[1];
174 for (int i = 0; i < 7; i++)
175 b1[i] = (a[1+i] - rc[7] * a[7-i]) / e;
176
177 for (int k = 6; k > 0; k--) {
178 b0 = b1, b1 = b[k & 1];
179
180 rc[k] = b0[k];
181 e = 1 - rc[k] * rc[k];
182
183 for (int i = 0; i < k; i++)
184 b1[i] = (b0[i] - rc[k] * b0[k-1-i]) / e;
185 }
186
187 rc[0] = b1[0];
188 }
189
190 /**
191 * Quantization of RC coefficients
192 * rc Refelection coefficients
193 * rc_order Return order of coefficients
194 * rc_i Return quantized coefficients
195 */
quantize_rc(const float * rc,int * rc_order,int * rc_q)196 static void quantize_rc(const float *rc, int *rc_order, int *rc_q)
197 {
198 /* Quantization table, sin(delta * (i + 0.5)), delta = Pi / 17 */
199
200 static float q_thr[] = {
201 9.22683595e-02, 2.73662990e-01, 4.45738356e-01, 6.02634636e-01,
202 7.39008917e-01, 8.50217136e-01, 9.32472229e-01, 9.82973100e-01
203 };
204
205 *rc_order = 8;
206
207 for (int i = 0; i < 8; i++) {
208 float rc_m = fabsf(rc[i]);
209
210 rc_q[i] = 4 * (rc_m >= q_thr[4]);
211 for (int j = 0; j < 4 && rc_m >= q_thr[rc_q[i]]; j++, rc_q[i]++);
212
213 if (rc[i] < 0)
214 rc_q[i] = -rc_q[i];
215
216 *rc_order = rc_q[i] != 0 ? 8 : *rc_order - 1;
217 }
218 }
219
220 /**
221 * Unquantization of RC coefficients
222 * rc_q Quantized coefficients
223 * rc_order Order of coefficients
224 * rc Return refelection coefficients
225 */
unquantize_rc(const int * rc_q,int rc_order,float rc[8])226 static void unquantize_rc(const int *rc_q, int rc_order, float rc[8])
227 {
228 /* Quantization table, sin(delta * i), delta = Pi / 17 */
229
230 static float q_inv[] = {
231 0.00000000e+00, 1.83749517e-01, 3.61241664e-01, 5.26432173e-01,
232 6.73695641e-01, 7.98017215e-01, 8.95163302e-01, 9.61825645e-01,
233 9.95734176e-01
234 };
235
236 int i;
237
238 for (i = 0; i < rc_order; i++) {
239 float rc_m = q_inv[LC3_ABS(rc_q[i])];
240 rc[i] = rc_q[i] < 0 ? -rc_m : rc_m;
241 }
242 }
243
244
245 /* ----------------------------------------------------------------------------
246 * Filtering
247 * -------------------------------------------------------------------------- */
248
249 /**
250 * Forward filtering
251 * dt, bw Duration and bandwidth of the frame
252 * rc_order, rc Order of coefficients, and coefficients
253 * x Spectral coefficients, filtered as output
254 */
forward_filtering(enum lc3_dt dt,enum lc3_bandwidth bw,const int rc_order[2],float (* const rc)[8],float * x)255 LC3_HOT static void forward_filtering(
256 enum lc3_dt dt, enum lc3_bandwidth bw,
257 const int rc_order[2], float (* const rc)[8], float *x)
258 {
259 int nfilters = 1 + (bw >= LC3_BANDWIDTH_SWB);
260 int nf = LC3_NE(dt, bw) >> (nfilters - 1);
261 int i0, ie = 3*(3 + dt);
262
263 float s[8] = { 0 };
264
265 for (int f = 0; f < nfilters; f++) {
266
267 i0 = ie;
268 ie = nf * (1 + f);
269
270 if (!rc_order[f])
271 continue;
272
273 for (int i = i0; i < ie; i++) {
274 float xi = x[i];
275 float s0, s1 = xi;
276
277 for (int k = 0; k < rc_order[f]; k++) {
278 s0 = s[k];
279 s[k] = s1;
280
281 s1 = rc[f][k] * xi + s0;
282 xi += rc[f][k] * s0;
283 }
284
285 x[i] = xi;
286 }
287 }
288 }
289
290 /**
291 * Inverse filtering
292 * dt, bw Duration and bandwidth of the frame
293 * rc_order, rc Order of coefficients, and unquantized coefficients
294 * x Spectral coefficients, filtered as output
295 */
inverse_filtering(enum lc3_dt dt,enum lc3_bandwidth bw,const int rc_order[2],float (* const rc)[8],float * x)296 LC3_HOT static void inverse_filtering(
297 enum lc3_dt dt, enum lc3_bandwidth bw,
298 const int rc_order[2], float (* const rc)[8], float *x)
299 {
300 int nfilters = 1 + (bw >= LC3_BANDWIDTH_SWB);
301 int nf = LC3_NE(dt, bw) >> (nfilters - 1);
302 int i0, ie = 3*(3 + dt);
303
304 float s[8] = { 0 };
305
306 for (int f = 0; f < nfilters; f++) {
307
308 i0 = ie;
309 ie = nf * (1 + f);
310
311 if (!rc_order[f])
312 continue;
313
314 for (int i = i0; i < ie; i++) {
315 float xi = x[i];
316
317 xi -= s[7] * rc[f][7];
318 for (int k = 6; k >= 0; k--) {
319 xi -= s[k] * rc[f][k];
320 s[k+1] = s[k] + rc[f][k] * xi;
321 }
322 s[0] = xi;
323 x[i] = xi;
324 }
325
326 for (int k = 7; k >= rc_order[f]; k--)
327 s[k] = 0;
328 }
329 }
330
331
332 /* ----------------------------------------------------------------------------
333 * Interface
334 * -------------------------------------------------------------------------- */
335
336 /**
337 * TNS analysis
338 */
lc3_tns_analyze(enum lc3_dt dt,enum lc3_bandwidth bw,bool nn_flag,int nbytes,struct lc3_tns_data * data,float * x)339 void lc3_tns_analyze(enum lc3_dt dt, enum lc3_bandwidth bw,
340 bool nn_flag, int nbytes, struct lc3_tns_data *data, float *x)
341 {
342 /* Processing steps :
343 * - Determine the LPC (Linear Predictive Coding) Coefficients
344 * - Check is the filtering is disabled
345 * - The coefficients are weighted on low bitrates and predicition gain
346 * - Convert to reflection coefficients and quantize
347 * - Finally filter the spectral coefficients */
348
349 float pred_gain[2], a[2][9];
350 float rc[2][8];
351
352 data->nfilters = 1 + (bw >= LC3_BANDWIDTH_SWB);
353 data->lpc_weighting = resolve_lpc_weighting(dt, nbytes);
354
355 compute_lpc_coeffs(dt, bw, x, pred_gain, a);
356
357 for (int f = 0; f < data->nfilters; f++) {
358
359 data->rc_order[f] = 0;
360 if (nn_flag || pred_gain[f] <= 1.5f)
361 continue;
362
363 if (data->lpc_weighting && pred_gain[f] < 2.f)
364 lpc_weighting(pred_gain[f], a[f]);
365
366 lpc_reflection(a[f], rc[f]);
367
368 quantize_rc(rc[f], &data->rc_order[f], data->rc[f]);
369 unquantize_rc(data->rc[f], data->rc_order[f], rc[f]);
370 }
371
372 forward_filtering(dt, bw, data->rc_order, rc, x);
373 }
374
375 /**
376 * TNS synthesis
377 */
lc3_tns_synthesize(enum lc3_dt dt,enum lc3_bandwidth bw,const struct lc3_tns_data * data,float * x)378 void lc3_tns_synthesize(enum lc3_dt dt, enum lc3_bandwidth bw,
379 const struct lc3_tns_data *data, float *x)
380 {
381 float rc[2][8] = { 0 };
382
383 for (int f = 0; f < data->nfilters; f++)
384 if (data->rc_order[f])
385 unquantize_rc(data->rc[f], data->rc_order[f], rc[f]);
386
387 inverse_filtering(dt, bw, data->rc_order, rc, x);
388 }
389
390 /**
391 * Bit consumption of bitstream data
392 */
lc3_tns_get_nbits(const struct lc3_tns_data * data)393 int lc3_tns_get_nbits(const struct lc3_tns_data *data)
394 {
395 int nbits = 0;
396
397 for (int f = 0; f < data->nfilters; f++) {
398
399 int nbits_2048 = 2048;
400 int rc_order = data->rc_order[f];
401
402 nbits_2048 += rc_order > 0 ? lc3_tns_order_bits
403 [data->lpc_weighting][rc_order-1] : 0;
404
405 for (int i = 0; i < rc_order; i++)
406 nbits_2048 += lc3_tns_coeffs_bits[i][8 + data->rc[f][i]];
407
408 nbits += (nbits_2048 + (1 << 11) - 1) >> 11;
409 }
410
411 return nbits;
412 }
413
414 /**
415 * Put bitstream data
416 */
lc3_tns_put_data(lc3_bits_t * bits,const struct lc3_tns_data * data)417 void lc3_tns_put_data(lc3_bits_t *bits, const struct lc3_tns_data *data)
418 {
419 for (int f = 0; f < data->nfilters; f++) {
420 int rc_order = data->rc_order[f];
421
422 lc3_put_bits(bits, rc_order > 0, 1);
423 if (rc_order <= 0)
424 continue;
425
426 lc3_put_symbol(bits,
427 lc3_tns_order_models + data->lpc_weighting, rc_order-1);
428
429 for (int i = 0; i < rc_order; i++)
430 lc3_put_symbol(bits,
431 lc3_tns_coeffs_models + i, 8 + data->rc[f][i]);
432 }
433 }
434
435 /**
436 * Get bitstream data
437 */
lc3_tns_get_data(lc3_bits_t * bits,enum lc3_dt dt,enum lc3_bandwidth bw,int nbytes,lc3_tns_data_t * data)438 void lc3_tns_get_data(lc3_bits_t *bits,
439 enum lc3_dt dt, enum lc3_bandwidth bw, int nbytes, lc3_tns_data_t *data)
440 {
441 data->nfilters = 1 + (bw >= LC3_BANDWIDTH_SWB);
442 data->lpc_weighting = resolve_lpc_weighting(dt, nbytes);
443
444 for (int f = 0; f < data->nfilters; f++) {
445
446 data->rc_order[f] = lc3_get_bit(bits);
447 if (!data->rc_order[f])
448 continue;
449
450 data->rc_order[f] += lc3_get_symbol(bits,
451 lc3_tns_order_models + data->lpc_weighting);
452
453 for (int i = 0; i < data->rc_order[f]; i++)
454 data->rc[f][i] = (int)lc3_get_symbol(bits,
455 lc3_tns_coeffs_models + i) - 8;
456 }
457 }
458
