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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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