1 /* Copyright (c) 2007-2008 CSIRO 2 Copyright (c) 2007-2009 Xiph.Org Foundation 3 Copyright (c) 2008-2009 Gregory Maxwell 4 Written by Jean-Marc Valin and Gregory Maxwell */ 5 /* 6 Redistribution and use in source and binary forms, with or without 7 modification, are permitted provided that the following conditions 8 are met: 9 10 - Redistributions of source code must retain the above copyright 11 notice, this list of conditions and the following disclaimer. 12 13 - Redistributions in binary form must reproduce the above copyright 14 notice, this list of conditions and the following disclaimer in the 15 documentation and/or other materials provided with the distribution. 16 17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 18 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 19 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 20 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 21 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 22 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 23 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 24 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 25 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 26 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 27 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 28 */ 29 30 #ifdef HAVE_CONFIG_H 31 #include "config.h" 32 #endif 33 34 #include <math.h> 35 #include "bands.h" 36 #include "modes.h" 37 #include "vq.h" 38 #include "cwrs.h" 39 #include "stack_alloc.h" 40 #include "os_support.h" 41 #include "mathops.h" 42 #include "rate.h" 43 #include "quant_bands.h" 44 #include "pitch.h" 45 46 int hysteresis_decision(opus_val16 val, const opus_val16 *thresholds, const opus_val16 *hysteresis, int N, int prev) 47 { 48 int i; 49 for (i=0;i<N;i++) 50 { 51 if (val < thresholds[i]) 52 break; 53 } 54 if (i>prev && val < thresholds[prev]+hysteresis[prev]) 55 i=prev; 56 if (i<prev && val > thresholds[prev-1]-hysteresis[prev-1]) 57 i=prev; 58 return i; 59 } 60 61 opus_uint32 celt_lcg_rand(opus_uint32 seed) 62 { 63 return 1664525 * seed + 1013904223; 64 } 65 66 /* This is a cos() approximation designed to be bit-exact on any platform. Bit exactness 67 with this approximation is important because it has an impact on the bit allocation */ 68 static opus_int16 bitexact_cos(opus_int16 x) 69 { 70 opus_int32 tmp; 71 opus_int16 x2; 72 tmp = (4096+((opus_int32)(x)*(x)))>>13; 73 celt_assert(tmp<=32767); 74 x2 = tmp; 75 x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2))))); 76 celt_assert(x2<=32766); 77 return 1+x2; 78 } 79 80 static int bitexact_log2tan(int isin,int icos) 81 { 82 int lc; 83 int ls; 84 lc=EC_ILOG(icos); 85 ls=EC_ILOG(isin); 86 icos<<=15-lc; 87 isin<<=15-ls; 88 return (ls-lc)*(1<<11) 89 +FRAC_MUL16(isin, FRAC_MUL16(isin, -2597) + 7932) 90 -FRAC_MUL16(icos, FRAC_MUL16(icos, -2597) + 7932); 91 } 92 93 #ifdef FIXED_POINT 94 /* Compute the amplitude (sqrt energy) in each of the bands */ 95 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int M) 96 { 97 int i, c, N; 98 const opus_int16 *eBands = m->eBands; 99 N = M*m->shortMdctSize; 100 c=0; do { 101 for (i=0;i<end;i++) 102 { 103 int j; 104 opus_val32 maxval=0; 105 opus_val32 sum = 0; 106 107 j=M*eBands[i]; do { 108 maxval = MAX32(maxval, X[j+c*N]); 109 maxval = MAX32(maxval, -X[j+c*N]); 110 } while (++j<M*eBands[i+1]); 111 112 if (maxval > 0) 113 { 114 int shift = celt_ilog2(maxval)-10; 115 j=M*eBands[i]; do { 116 sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)), 117 EXTRACT16(VSHR32(X[j+c*N],shift))); 118 } while (++j<M*eBands[i+1]); 119 /* We're adding one here to ensure the normalized band isn't larger than unity norm */ 120 bandE[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift); 121 } else { 122 bandE[i+c*m->nbEBands] = EPSILON; 123 } 124 /*printf ("%f ", bandE[i+c*m->nbEBands]);*/ 125 } 126 } while (++c<C); 127 /*printf ("\n");*/ 128 } 129 130 /* Normalise each band such that the energy is one. */ 131 void normalise_bands(const CELTMode *m, const celt_sig * OPUS_RESTRICT freq, celt_norm * OPUS_RESTRICT X, const celt_ener *bandE, int end, int C, int M) 132 { 133 int i, c, N; 134 const opus_int16 *eBands = m->eBands; 135 N = M*m->shortMdctSize; 136 c=0; do { 137 i=0; do { 138 opus_val16 g; 139 int j,shift; 140 opus_val16 E; 141 shift = celt_zlog2(bandE[i+c*m->nbEBands])-13; 142 E = VSHR32(bandE[i+c*m->nbEBands], shift); 143 g = EXTRACT16(celt_rcp(SHL32(E,3))); 144 j=M*eBands[i]; do { 145 X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g); 146 } while (++j<M*eBands[i+1]); 147 } while (++i<end); 148 } while (++c<C); 149 } 150 151 #else /* FIXED_POINT */ 152 /* Compute the amplitude (sqrt energy) in each of the bands */ 153 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int M) 154 { 155 int i, c, N; 156 const opus_int16 *eBands = m->eBands; 157 N = M*m->shortMdctSize; 158 c=0; do { 159 for (i=0;i<end;i++) 160 { 161 int j; 162 opus_val32 sum = 1e-27f; 163 for (j=M*eBands[i];j<M*eBands[i+1];j++) 164 sum += X[j+c*N]*X[j+c*N]; 165 bandE[i+c*m->nbEBands] = celt_sqrt(sum); 166 /*printf ("%f ", bandE[i+c*m->nbEBands]);*/ 167 } 168 } while (++c<C); 169 /*printf ("\n");*/ 170 } 171 172 /* Normalise each band such that the energy is one. */ 173 void normalise_bands(const CELTMode *m, const celt_sig * OPUS_RESTRICT freq, celt_norm * OPUS_RESTRICT X, const celt_ener *bandE, int end, int C, int M) 174 { 175 int i, c, N; 176 const opus_int16 *eBands = m->eBands; 177 N = M*m->shortMdctSize; 178 c=0; do { 179 for (i=0;i<end;i++) 180 { 181 int j; 182 opus_val16 g = 1.f/(1e-27f+bandE[i+c*m->nbEBands]); 183 for (j=M*eBands[i];j<M*eBands[i+1];j++) 184 X[j+c*N] = freq[j+c*N]*g; 185 } 186 } while (++c<C); 187 } 188 189 #endif /* FIXED_POINT */ 190 191 /* De-normalise the energy to produce the synthesis from the unit-energy bands */ 192 void denormalise_bands(const CELTMode *m, const celt_norm * OPUS_RESTRICT X, 193 celt_sig * OPUS_RESTRICT freq, const opus_val16 *bandLogE, int start, int end, int C, int M) 194 { 195 int i, c, N; 196 const opus_int16 *eBands = m->eBands; 197 N = M*m->shortMdctSize; 198 celt_assert2(C<=2, "denormalise_bands() not implemented for >2 channels"); 199 c=0; do { 200 celt_sig * OPUS_RESTRICT f; 201 const celt_norm * OPUS_RESTRICT x; 202 f = freq+c*N; 203 x = X+c*N+M*eBands[start]; 204 for (i=0;i<M*eBands[start];i++) 205 *f++ = 0; 206 for (i=start;i<end;i++) 207 { 208 int j, band_end; 209 opus_val16 g; 210 opus_val16 lg; 211 #ifdef FIXED_POINT 212 int shift; 213 #endif 214 j=M*eBands[i]; 215 band_end = M*eBands[i+1]; 216 lg = ADD16(bandLogE[i+c*m->nbEBands], SHL16((opus_val16)eMeans[i],6)); 217 #ifdef FIXED_POINT 218 /* Handle the integer part of the log energy */ 219 shift = 16-(lg>>DB_SHIFT); 220 if (shift>31) 221 { 222 shift=0; 223 g=0; 224 } else { 225 /* Handle the fractional part. */ 226 g = celt_exp2_frac(lg&((1<<DB_SHIFT)-1)); 227 } 228 #else 229 g = celt_exp2(lg); 230 #endif 231 do { 232 *f++ = SHR32(MULT16_16(*x++, g), shift); 233 } while (++j<band_end); 234 } 235 celt_assert(start <= end); 236 for (i=M*eBands[end];i<N;i++) 237 *f++ = 0; 238 } while (++c<C); 239 } 240 241 /* This prevents energy collapse for transients with multiple short MDCTs */ 242 void anti_collapse(const CELTMode *m, celt_norm *X_, unsigned char *collapse_masks, int LM, int C, int size, 243 int start, int end, opus_val16 *logE, opus_val16 *prev1logE, 244 opus_val16 *prev2logE, int *pulses, opus_uint32 seed) 245 { 246 int c, i, j, k; 247 for (i=start;i<end;i++) 248 { 249 int N0; 250 opus_val16 thresh, sqrt_1; 251 int depth; 252 #ifdef FIXED_POINT 253 int shift; 254 opus_val32 thresh32; 255 #endif 256 257 N0 = m->eBands[i+1]-m->eBands[i]; 258 /* depth in 1/8 bits */ 259 depth = (1+pulses[i])/((m->eBands[i+1]-m->eBands[i])<<LM); 260 261 #ifdef FIXED_POINT 262 thresh32 = SHR32(celt_exp2(-SHL16(depth, 10-BITRES)),1); 263 thresh = MULT16_32_Q15(QCONST16(0.5f, 15), MIN32(32767,thresh32)); 264 { 265 opus_val32 t; 266 t = N0<<LM; 267 shift = celt_ilog2(t)>>1; 268 t = SHL32(t, (7-shift)<<1); 269 sqrt_1 = celt_rsqrt_norm(t); 270 } 271 #else 272 thresh = .5f*celt_exp2(-.125f*depth); 273 sqrt_1 = celt_rsqrt(N0<<LM); 274 #endif 275 276 c=0; do 277 { 278 celt_norm *X; 279 opus_val16 prev1; 280 opus_val16 prev2; 281 opus_val32 Ediff; 282 opus_val16 r; 283 int renormalize=0; 284 prev1 = prev1logE[c*m->nbEBands+i]; 285 prev2 = prev2logE[c*m->nbEBands+i]; 286 if (C==1) 287 { 288 prev1 = MAX16(prev1,prev1logE[m->nbEBands+i]); 289 prev2 = MAX16(prev2,prev2logE[m->nbEBands+i]); 290 } 291 Ediff = EXTEND32(logE[c*m->nbEBands+i])-EXTEND32(MIN16(prev1,prev2)); 292 Ediff = MAX32(0, Ediff); 293 294 #ifdef FIXED_POINT 295 if (Ediff < 16384) 296 { 297 opus_val32 r32 = SHR32(celt_exp2(-EXTRACT16(Ediff)),1); 298 r = 2*MIN16(16383,r32); 299 } else { 300 r = 0; 301 } 302 if (LM==3) 303 r = MULT16_16_Q14(23170, MIN32(23169, r)); 304 r = SHR16(MIN16(thresh, r),1); 305 r = SHR32(MULT16_16_Q15(sqrt_1, r),shift); 306 #else 307 /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because 308 short blocks don't have the same energy as long */ 309 r = 2.f*celt_exp2(-Ediff); 310 if (LM==3) 311 r *= 1.41421356f; 312 r = MIN16(thresh, r); 313 r = r*sqrt_1; 314 #endif 315 X = X_+c*size+(m->eBands[i]<<LM); 316 for (k=0;k<1<<LM;k++) 317 { 318 /* Detect collapse */ 319 if (!(collapse_masks[i*C+c]&1<<k)) 320 { 321 /* Fill with noise */ 322 for (j=0;j<N0;j++) 323 { 324 seed = celt_lcg_rand(seed); 325 X[(j<<LM)+k] = (seed&0x8000 ? r : -r); 326 } 327 renormalize = 1; 328 } 329 } 330 /* We just added some energy, so we need to renormalise */ 331 if (renormalize) 332 renormalise_vector(X, N0<<LM, Q15ONE); 333 } while (++c<C); 334 } 335 } 336 337 static void intensity_stereo(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bandE, int bandID, int N) 338 { 339 int i = bandID; 340 int j; 341 opus_val16 a1, a2; 342 opus_val16 left, right; 343 opus_val16 norm; 344 #ifdef FIXED_POINT 345 int shift = celt_zlog2(MAX32(bandE[i], bandE[i+m->nbEBands]))-13; 346 #endif 347 left = VSHR32(bandE[i],shift); 348 right = VSHR32(bandE[i+m->nbEBands],shift); 349 norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right)); 350 a1 = DIV32_16(SHL32(EXTEND32(left),14),norm); 351 a2 = DIV32_16(SHL32(EXTEND32(right),14),norm); 352 for (j=0;j<N;j++) 353 { 354 celt_norm r, l; 355 l = X[j]; 356 r = Y[j]; 357 X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r); 358 /* Side is not encoded, no need to calculate */ 359 } 360 } 361 362 static void stereo_split(celt_norm *X, celt_norm *Y, int N) 363 { 364 int j; 365 for (j=0;j<N;j++) 366 { 367 celt_norm r, l; 368 l = MULT16_16_Q15(QCONST16(.70710678f,15), X[j]); 369 r = MULT16_16_Q15(QCONST16(.70710678f,15), Y[j]); 370 X[j] = l+r; 371 Y[j] = r-l; 372 } 373 } 374 375 static void stereo_merge(celt_norm *X, celt_norm *Y, opus_val16 mid, int N) 376 { 377 int j; 378 opus_val32 xp=0, side=0; 379 opus_val32 El, Er; 380 opus_val16 mid2; 381 #ifdef FIXED_POINT 382 int kl, kr; 383 #endif 384 opus_val32 t, lgain, rgain; 385 386 /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */ 387 dual_inner_prod(Y, X, Y, N, &xp, &side); 388 /* Compensating for the mid normalization */ 389 xp = MULT16_32_Q15(mid, xp); 390 /* mid and side are in Q15, not Q14 like X and Y */ 391 mid2 = SHR32(mid, 1); 392 El = MULT16_16(mid2, mid2) + side - 2*xp; 393 Er = MULT16_16(mid2, mid2) + side + 2*xp; 394 if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28)) 395 { 396 for (j=0;j<N;j++) 397 Y[j] = X[j]; 398 return; 399 } 400 401 #ifdef FIXED_POINT 402 kl = celt_ilog2(El)>>1; 403 kr = celt_ilog2(Er)>>1; 404 #endif 405 t = VSHR32(El, (kl-7)<<1); 406 lgain = celt_rsqrt_norm(t); 407 t = VSHR32(Er, (kr-7)<<1); 408 rgain = celt_rsqrt_norm(t); 409 410 #ifdef FIXED_POINT 411 if (kl < 7) 412 kl = 7; 413 if (kr < 7) 414 kr = 7; 415 #endif 416 417 for (j=0;j<N;j++) 418 { 419 celt_norm r, l; 420 /* Apply mid scaling (side is already scaled) */ 421 l = MULT16_16_Q15(mid, X[j]); 422 r = Y[j]; 423 X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1)); 424 Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1)); 425 } 426 } 427 428 /* Decide whether we should spread the pulses in the current frame */ 429 int spreading_decision(const CELTMode *m, celt_norm *X, int *average, 430 int last_decision, int *hf_average, int *tapset_decision, int update_hf, 431 int end, int C, int M) 432 { 433 int i, c, N0; 434 int sum = 0, nbBands=0; 435 const opus_int16 * OPUS_RESTRICT eBands = m->eBands; 436 int decision; 437 int hf_sum=0; 438 439 celt_assert(end>0); 440 441 N0 = M*m->shortMdctSize; 442 443 if (M*(eBands[end]-eBands[end-1]) <= 8) 444 return SPREAD_NONE; 445 c=0; do { 446 for (i=0;i<end;i++) 447 { 448 int j, N, tmp=0; 449 int tcount[3] = {0,0,0}; 450 celt_norm * OPUS_RESTRICT x = X+M*eBands[i]+c*N0; 451 N = M*(eBands[i+1]-eBands[i]); 452 if (N<=8) 453 continue; 454 /* Compute rough CDF of |x[j]| */ 455 for (j=0;j<N;j++) 456 { 457 opus_val32 x2N; /* Q13 */ 458 459 x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N); 460 if (x2N < QCONST16(0.25f,13)) 461 tcount[0]++; 462 if (x2N < QCONST16(0.0625f,13)) 463 tcount[1]++; 464 if (x2N < QCONST16(0.015625f,13)) 465 tcount[2]++; 466 } 467 468 /* Only include four last bands (8 kHz and up) */ 469 if (i>m->nbEBands-4) 470 hf_sum += 32*(tcount[1]+tcount[0])/N; 471 tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N); 472 sum += tmp*256; 473 nbBands++; 474 } 475 } while (++c<C); 476 477 if (update_hf) 478 { 479 if (hf_sum) 480 hf_sum /= C*(4-m->nbEBands+end); 481 *hf_average = (*hf_average+hf_sum)>>1; 482 hf_sum = *hf_average; 483 if (*tapset_decision==2) 484 hf_sum += 4; 485 else if (*tapset_decision==0) 486 hf_sum -= 4; 487 if (hf_sum > 22) 488 *tapset_decision=2; 489 else if (hf_sum > 18) 490 *tapset_decision=1; 491 else 492 *tapset_decision=0; 493 } 494 /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/ 495 celt_assert(nbBands>0); /*M*(eBands[end]-eBands[end-1]) <= 8 assures this*/ 496 sum /= nbBands; 497 /* Recursive averaging */ 498 sum = (sum+*average)>>1; 499 *average = sum; 500 /* Hysteresis */ 501 sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2; 502 if (sum < 80) 503 { 504 decision = SPREAD_AGGRESSIVE; 505 } else if (sum < 256) 506 { 507 decision = SPREAD_NORMAL; 508 } else if (sum < 384) 509 { 510 decision = SPREAD_LIGHT; 511 } else { 512 decision = SPREAD_NONE; 513 } 514 #ifdef FUZZING 515 decision = rand()&0x3; 516 *tapset_decision=rand()%3; 517 #endif 518 return decision; 519 } 520 521 /* Indexing table for converting from natural Hadamard to ordery Hadamard 522 This is essentially a bit-reversed Gray, on top of which we've added 523 an inversion of the order because we want the DC at the end rather than 524 the beginning. The lines are for N=2, 4, 8, 16 */ 525 static const int ordery_table[] = { 526 1, 0, 527 3, 0, 2, 1, 528 7, 0, 4, 3, 6, 1, 5, 2, 529 15, 0, 8, 7, 12, 3, 11, 4, 14, 1, 9, 6, 13, 2, 10, 5, 530 }; 531 532 static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard) 533 { 534 int i,j; 535 VARDECL(celt_norm, tmp); 536 int N; 537 SAVE_STACK; 538 N = N0*stride; 539 ALLOC(tmp, N, celt_norm); 540 celt_assert(stride>0); 541 if (hadamard) 542 { 543 const int *ordery = ordery_table+stride-2; 544 for (i=0;i<stride;i++) 545 { 546 for (j=0;j<N0;j++) 547 tmp[ordery[i]*N0+j] = X[j*stride+i]; 548 } 549 } else { 550 for (i=0;i<stride;i++) 551 for (j=0;j<N0;j++) 552 tmp[i*N0+j] = X[j*stride+i]; 553 } 554 for (j=0;j<N;j++) 555 X[j] = tmp[j]; 556 RESTORE_STACK; 557 } 558 559 static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard) 560 { 561 int i,j; 562 VARDECL(celt_norm, tmp); 563 int N; 564 SAVE_STACK; 565 N = N0*stride; 566 ALLOC(tmp, N, celt_norm); 567 if (hadamard) 568 { 569 const int *ordery = ordery_table+stride-2; 570 for (i=0;i<stride;i++) 571 for (j=0;j<N0;j++) 572 tmp[j*stride+i] = X[ordery[i]*N0+j]; 573 } else { 574 for (i=0;i<stride;i++) 575 for (j=0;j<N0;j++) 576 tmp[j*stride+i] = X[i*N0+j]; 577 } 578 for (j=0;j<N;j++) 579 X[j] = tmp[j]; 580 RESTORE_STACK; 581 } 582 583 void haar1(celt_norm *X, int N0, int stride) 584 { 585 int i, j; 586 N0 >>= 1; 587 for (i=0;i<stride;i++) 588 for (j=0;j<N0;j++) 589 { 590 celt_norm tmp1, tmp2; 591 tmp1 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*2*j+i]); 592 tmp2 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*(2*j+1)+i]); 593 X[stride*2*j+i] = tmp1 + tmp2; 594 X[stride*(2*j+1)+i] = tmp1 - tmp2; 595 } 596 } 597 598 static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo) 599 { 600 static const opus_int16 exp2_table8[8] = 601 {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048}; 602 int qn, qb; 603 int N2 = 2*N-1; 604 if (stereo && N==2) 605 N2--; 606 /* The upper limit ensures that in a stereo split with itheta==16384, we'll 607 always have enough bits left over to code at least one pulse in the 608 side; otherwise it would collapse, since it doesn't get folded. */ 609 qb = IMIN(b-pulse_cap-(4<<BITRES), (b+N2*offset)/N2); 610 611 qb = IMIN(8<<BITRES, qb); 612 613 if (qb<(1<<BITRES>>1)) { 614 qn = 1; 615 } else { 616 qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES)); 617 qn = (qn+1)>>1<<1; 618 } 619 celt_assert(qn <= 256); 620 return qn; 621 } 622 623 struct band_ctx { 624 int encode; 625 const CELTMode *m; 626 int i; 627 int intensity; 628 int spread; 629 int tf_change; 630 ec_ctx *ec; 631 opus_int32 remaining_bits; 632 const celt_ener *bandE; 633 opus_uint32 seed; 634 }; 635 636 struct split_ctx { 637 int inv; 638 int imid; 639 int iside; 640 int delta; 641 int itheta; 642 int qalloc; 643 }; 644 645 static void compute_theta(struct band_ctx *ctx, struct split_ctx *sctx, 646 celt_norm *X, celt_norm *Y, int N, int *b, int B, int B0, 647 int LM, 648 int stereo, int *fill) 649 { 650 int qn; 651 int itheta=0; 652 int delta; 653 int imid, iside; 654 int qalloc; 655 int pulse_cap; 656 int offset; 657 opus_int32 tell; 658 int inv=0; 659 int encode; 660 const CELTMode *m; 661 int i; 662 int intensity; 663 ec_ctx *ec; 664 const celt_ener *bandE; 665 666 encode = ctx->encode; 667 m = ctx->m; 668 i = ctx->i; 669 intensity = ctx->intensity; 670 ec = ctx->ec; 671 bandE = ctx->bandE; 672 673 /* Decide on the resolution to give to the split parameter theta */ 674 pulse_cap = m->logN[i]+LM*(1<<BITRES); 675 offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET); 676 qn = compute_qn(N, *b, offset, pulse_cap, stereo); 677 if (stereo && i>=intensity) 678 qn = 1; 679 if (encode) 680 { 681 /* theta is the atan() of the ratio between the (normalized) 682 side and mid. With just that parameter, we can re-scale both 683 mid and side because we know that 1) they have unit norm and 684 2) they are orthogonal. */ 685 itheta = stereo_itheta(X, Y, stereo, N); 686 } 687 tell = ec_tell_frac(ec); 688 if (qn!=1) 689 { 690 if (encode) 691 itheta = (itheta*qn+8192)>>14; 692 693 /* Entropy coding of the angle. We use a uniform pdf for the 694 time split, a step for stereo, and a triangular one for the rest. */ 695 if (stereo && N>2) 696 { 697 int p0 = 3; 698 int x = itheta; 699 int x0 = qn/2; 700 int ft = p0*(x0+1) + x0; 701 /* Use a probability of p0 up to itheta=8192 and then use 1 after */ 702 if (encode) 703 { 704 ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft); 705 } else { 706 int fs; 707 fs=ec_decode(ec,ft); 708 if (fs<(x0+1)*p0) 709 x=fs/p0; 710 else 711 x=x0+1+(fs-(x0+1)*p0); 712 ec_dec_update(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft); 713 itheta = x; 714 } 715 } else if (B0>1 || stereo) { 716 /* Uniform pdf */ 717 if (encode) 718 ec_enc_uint(ec, itheta, qn+1); 719 else 720 itheta = ec_dec_uint(ec, qn+1); 721 } else { 722 int fs=1, ft; 723 ft = ((qn>>1)+1)*((qn>>1)+1); 724 if (encode) 725 { 726 int fl; 727 728 fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta; 729 fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 : 730 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1); 731 732 ec_encode(ec, fl, fl+fs, ft); 733 } else { 734 /* Triangular pdf */ 735 int fl=0; 736 int fm; 737 fm = ec_decode(ec, ft); 738 739 if (fm < ((qn>>1)*((qn>>1) + 1)>>1)) 740 { 741 itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1; 742 fs = itheta + 1; 743 fl = itheta*(itheta + 1)>>1; 744 } 745 else 746 { 747 itheta = (2*(qn + 1) 748 - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1; 749 fs = qn + 1 - itheta; 750 fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1); 751 } 752 753 ec_dec_update(ec, fl, fl+fs, ft); 754 } 755 } 756 itheta = (opus_int32)itheta*16384/qn; 757 if (encode && stereo) 758 { 759 if (itheta==0) 760 intensity_stereo(m, X, Y, bandE, i, N); 761 else 762 stereo_split(X, Y, N); 763 } 764 /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate. 765 Let's do that at higher complexity */ 766 } else if (stereo) { 767 if (encode) 768 { 769 inv = itheta > 8192; 770 if (inv) 771 { 772 int j; 773 for (j=0;j<N;j++) 774 Y[j] = -Y[j]; 775 } 776 intensity_stereo(m, X, Y, bandE, i, N); 777 } 778 if (*b>2<<BITRES && ctx->remaining_bits > 2<<BITRES) 779 { 780 if (encode) 781 ec_enc_bit_logp(ec, inv, 2); 782 else 783 inv = ec_dec_bit_logp(ec, 2); 784 } else 785 inv = 0; 786 itheta = 0; 787 } 788 qalloc = ec_tell_frac(ec) - tell; 789 *b -= qalloc; 790 791 if (itheta == 0) 792 { 793 imid = 32767; 794 iside = 0; 795 *fill &= (1<<B)-1; 796 delta = -16384; 797 } else if (itheta == 16384) 798 { 799 imid = 0; 800 iside = 32767; 801 *fill &= ((1<<B)-1)<<B; 802 delta = 16384; 803 } else { 804 imid = bitexact_cos((opus_int16)itheta); 805 iside = bitexact_cos((opus_int16)(16384-itheta)); 806 /* This is the mid vs side allocation that minimizes squared error 807 in that band. */ 808 delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid)); 809 } 810 811 sctx->inv = inv; 812 sctx->imid = imid; 813 sctx->iside = iside; 814 sctx->delta = delta; 815 sctx->itheta = itheta; 816 sctx->qalloc = qalloc; 817 } 818 static unsigned quant_band_n1(struct band_ctx *ctx, celt_norm *X, celt_norm *Y, int b, 819 celt_norm *lowband_out) 820 { 821 #ifdef RESYNTH 822 int resynth = 1; 823 #else 824 int resynth = !ctx->encode; 825 #endif 826 int c; 827 int stereo; 828 celt_norm *x = X; 829 int encode; 830 ec_ctx *ec; 831 832 encode = ctx->encode; 833 ec = ctx->ec; 834 835 stereo = Y != NULL; 836 c=0; do { 837 int sign=0; 838 if (ctx->remaining_bits>=1<<BITRES) 839 { 840 if (encode) 841 { 842 sign = x[0]<0; 843 ec_enc_bits(ec, sign, 1); 844 } else { 845 sign = ec_dec_bits(ec, 1); 846 } 847 ctx->remaining_bits -= 1<<BITRES; 848 b-=1<<BITRES; 849 } 850 if (resynth) 851 x[0] = sign ? -NORM_SCALING : NORM_SCALING; 852 x = Y; 853 } while (++c<1+stereo); 854 if (lowband_out) 855 lowband_out[0] = SHR16(X[0],4); 856 return 1; 857 } 858 859 /* This function is responsible for encoding and decoding a mono partition. 860 It can split the band in two and transmit the energy difference with 861 the two half-bands. It can be called recursively so bands can end up being 862 split in 8 parts. */ 863 static unsigned quant_partition(struct band_ctx *ctx, celt_norm *X, 864 int N, int b, int B, celt_norm *lowband, 865 int LM, 866 opus_val16 gain, int fill) 867 { 868 const unsigned char *cache; 869 int q; 870 int curr_bits; 871 int imid=0, iside=0; 872 int N_B=N; 873 int B0=B; 874 opus_val16 mid=0, side=0; 875 unsigned cm=0; 876 #ifdef RESYNTH 877 int resynth = 1; 878 #else 879 int resynth = !ctx->encode; 880 #endif 881 celt_norm *Y=NULL; 882 int encode; 883 const CELTMode *m; 884 int i; 885 int spread; 886 ec_ctx *ec; 887 888 encode = ctx->encode; 889 m = ctx->m; 890 i = ctx->i; 891 spread = ctx->spread; 892 ec = ctx->ec; 893 894 N_B /= B; 895 896 /* If we need 1.5 more bit than we can produce, split the band in two. */ 897 cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i]; 898 if (LM != -1 && b > cache[cache[0]]+12 && N>2) 899 { 900 int mbits, sbits, delta; 901 int itheta; 902 int qalloc; 903 struct split_ctx sctx; 904 celt_norm *next_lowband2=NULL; 905 opus_int32 rebalance; 906 907 N >>= 1; 908 Y = X+N; 909 LM -= 1; 910 if (B==1) 911 fill = (fill&1)|(fill<<1); 912 B = (B+1)>>1; 913 914 compute_theta(ctx, &sctx, X, Y, N, &b, B, B0, 915 LM, 0, &fill); 916 imid = sctx.imid; 917 iside = sctx.iside; 918 delta = sctx.delta; 919 itheta = sctx.itheta; 920 qalloc = sctx.qalloc; 921 #ifdef FIXED_POINT 922 mid = imid; 923 side = iside; 924 #else 925 mid = (1.f/32768)*imid; 926 side = (1.f/32768)*iside; 927 #endif 928 929 /* Give more bits to low-energy MDCTs than they would otherwise deserve */ 930 if (B0>1 && (itheta&0x3fff)) 931 { 932 if (itheta > 8192) 933 /* Rough approximation for pre-echo masking */ 934 delta -= delta>>(4-LM); 935 else 936 /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */ 937 delta = IMIN(0, delta + (N<<BITRES>>(5-LM))); 938 } 939 mbits = IMAX(0, IMIN(b, (b-delta)/2)); 940 sbits = b-mbits; 941 ctx->remaining_bits -= qalloc; 942 943 if (lowband) 944 next_lowband2 = lowband+N; /* >32-bit split case */ 945 946 rebalance = ctx->remaining_bits; 947 if (mbits >= sbits) 948 { 949 cm = quant_partition(ctx, X, N, mbits, B, 950 lowband, LM, 951 MULT16_16_P15(gain,mid), fill); 952 rebalance = mbits - (rebalance-ctx->remaining_bits); 953 if (rebalance > 3<<BITRES && itheta!=0) 954 sbits += rebalance - (3<<BITRES); 955 cm |= quant_partition(ctx, Y, N, sbits, B, 956 next_lowband2, LM, 957 MULT16_16_P15(gain,side), fill>>B)<<(B0>>1); 958 } else { 959 cm = quant_partition(ctx, Y, N, sbits, B, 960 next_lowband2, LM, 961 MULT16_16_P15(gain,side), fill>>B)<<(B0>>1); 962 rebalance = sbits - (rebalance-ctx->remaining_bits); 963 if (rebalance > 3<<BITRES && itheta!=16384) 964 mbits += rebalance - (3<<BITRES); 965 cm |= quant_partition(ctx, X, N, mbits, B, 966 lowband, LM, 967 MULT16_16_P15(gain,mid), fill); 968 } 969 } else { 970 /* This is the basic no-split case */ 971 q = bits2pulses(m, i, LM, b); 972 curr_bits = pulses2bits(m, i, LM, q); 973 ctx->remaining_bits -= curr_bits; 974 975 /* Ensures we can never bust the budget */ 976 while (ctx->remaining_bits < 0 && q > 0) 977 { 978 ctx->remaining_bits += curr_bits; 979 q--; 980 curr_bits = pulses2bits(m, i, LM, q); 981 ctx->remaining_bits -= curr_bits; 982 } 983 984 if (q!=0) 985 { 986 int K = get_pulses(q); 987 988 /* Finally do the actual quantization */ 989 if (encode) 990 { 991 cm = alg_quant(X, N, K, spread, B, ec 992 #ifdef RESYNTH 993 , gain 994 #endif 995 ); 996 } else { 997 cm = alg_unquant(X, N, K, spread, B, ec, gain); 998 } 999 } else { 1000 /* If there's no pulse, fill the band anyway */ 1001 int j; 1002 if (resynth) 1003 { 1004 unsigned cm_mask; 1005 /* B can be as large as 16, so this shift might overflow an int on a 1006 16-bit platform; use a long to get defined behavior.*/ 1007 cm_mask = (unsigned)(1UL<<B)-1; 1008 fill &= cm_mask; 1009 if (!fill) 1010 { 1011 for (j=0;j<N;j++) 1012 X[j] = 0; 1013 } else { 1014 if (lowband == NULL) 1015 { 1016 /* Noise */ 1017 for (j=0;j<N;j++) 1018 { 1019 ctx->seed = celt_lcg_rand(ctx->seed); 1020 X[j] = (celt_norm)((opus_int32)ctx->seed>>20); 1021 } 1022 cm = cm_mask; 1023 } else { 1024 /* Folded spectrum */ 1025 for (j=0;j<N;j++) 1026 { 1027 opus_val16 tmp; 1028 ctx->seed = celt_lcg_rand(ctx->seed); 1029 /* About 48 dB below the "normal" folding level */ 1030 tmp = QCONST16(1.0f/256, 10); 1031 tmp = (ctx->seed)&0x8000 ? tmp : -tmp; 1032 X[j] = lowband[j]+tmp; 1033 } 1034 cm = fill; 1035 } 1036 renormalise_vector(X, N, gain); 1037 } 1038 } 1039 } 1040 } 1041 1042 return cm; 1043 } 1044 1045 1046 /* This function is responsible for encoding and decoding a band for the mono case. */ 1047 static unsigned quant_band(struct band_ctx *ctx, celt_norm *X, 1048 int N, int b, int B, celt_norm *lowband, 1049 int LM, celt_norm *lowband_out, 1050 opus_val16 gain, celt_norm *lowband_scratch, int fill) 1051 { 1052 int N0=N; 1053 int N_B=N; 1054 int N_B0; 1055 int B0=B; 1056 int time_divide=0; 1057 int recombine=0; 1058 int longBlocks; 1059 unsigned cm=0; 1060 #ifdef RESYNTH 1061 int resynth = 1; 1062 #else 1063 int resynth = !ctx->encode; 1064 #endif 1065 int k; 1066 int encode; 1067 int tf_change; 1068 1069 encode = ctx->encode; 1070 tf_change = ctx->tf_change; 1071 1072 longBlocks = B0==1; 1073 1074 N_B /= B; 1075 N_B0 = N_B; 1076 1077 /* Special case for one sample */ 1078 if (N==1) 1079 { 1080 return quant_band_n1(ctx, X, NULL, b, lowband_out); 1081 } 1082 1083 if (tf_change>0) 1084 recombine = tf_change; 1085 /* Band recombining to increase frequency resolution */ 1086 1087 if (lowband_scratch && lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1)) 1088 { 1089 int j; 1090 for (j=0;j<N;j++) 1091 lowband_scratch[j] = lowband[j]; 1092 lowband = lowband_scratch; 1093 } 1094 1095 for (k=0;k<recombine;k++) 1096 { 1097 static const unsigned char bit_interleave_table[16]={ 1098 0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3 1099 }; 1100 if (encode) 1101 haar1(X, N>>k, 1<<k); 1102 if (lowband) 1103 haar1(lowband, N>>k, 1<<k); 1104 fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2; 1105 } 1106 B>>=recombine; 1107 N_B<<=recombine; 1108 1109 /* Increasing the time resolution */ 1110 while ((N_B&1) == 0 && tf_change<0) 1111 { 1112 if (encode) 1113 haar1(X, N_B, B); 1114 if (lowband) 1115 haar1(lowband, N_B, B); 1116 fill |= fill<<B; 1117 B <<= 1; 1118 N_B >>= 1; 1119 time_divide++; 1120 tf_change++; 1121 } 1122 B0=B; 1123 N_B0 = N_B; 1124 1125 /* Reorganize the samples in time order instead of frequency order */ 1126 if (B0>1) 1127 { 1128 if (encode) 1129 deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks); 1130 if (lowband) 1131 deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks); 1132 } 1133 1134 cm = quant_partition(ctx, X, N, b, B, lowband, 1135 LM, gain, fill); 1136 1137 /* This code is used by the decoder and by the resynthesis-enabled encoder */ 1138 if (resynth) 1139 { 1140 /* Undo the sample reorganization going from time order to frequency order */ 1141 if (B0>1) 1142 interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks); 1143 1144 /* Undo time-freq changes that we did earlier */ 1145 N_B = N_B0; 1146 B = B0; 1147 for (k=0;k<time_divide;k++) 1148 { 1149 B >>= 1; 1150 N_B <<= 1; 1151 cm |= cm>>B; 1152 haar1(X, N_B, B); 1153 } 1154 1155 for (k=0;k<recombine;k++) 1156 { 1157 static const unsigned char bit_deinterleave_table[16]={ 1158 0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F, 1159 0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF 1160 }; 1161 cm = bit_deinterleave_table[cm]; 1162 haar1(X, N0>>k, 1<<k); 1163 } 1164 B<<=recombine; 1165 1166 /* Scale output for later folding */ 1167 if (lowband_out) 1168 { 1169 int j; 1170 opus_val16 n; 1171 n = celt_sqrt(SHL32(EXTEND32(N0),22)); 1172 for (j=0;j<N0;j++) 1173 lowband_out[j] = MULT16_16_Q15(n,X[j]); 1174 } 1175 cm &= (1<<B)-1; 1176 } 1177 return cm; 1178 } 1179 1180 1181 /* This function is responsible for encoding and decoding a band for the stereo case. */ 1182 static unsigned quant_band_stereo(struct band_ctx *ctx, celt_norm *X, celt_norm *Y, 1183 int N, int b, int B, celt_norm *lowband, 1184 int LM, celt_norm *lowband_out, 1185 celt_norm *lowband_scratch, int fill) 1186 { 1187 int imid=0, iside=0; 1188 int inv = 0; 1189 opus_val16 mid=0, side=0; 1190 unsigned cm=0; 1191 #ifdef RESYNTH 1192 int resynth = 1; 1193 #else 1194 int resynth = !ctx->encode; 1195 #endif 1196 int mbits, sbits, delta; 1197 int itheta; 1198 int qalloc; 1199 struct split_ctx sctx; 1200 int orig_fill; 1201 int encode; 1202 ec_ctx *ec; 1203 1204 encode = ctx->encode; 1205 ec = ctx->ec; 1206 1207 /* Special case for one sample */ 1208 if (N==1) 1209 { 1210 return quant_band_n1(ctx, X, Y, b, lowband_out); 1211 } 1212 1213 orig_fill = fill; 1214 1215 compute_theta(ctx, &sctx, X, Y, N, &b, B, B, 1216 LM, 1, &fill); 1217 inv = sctx.inv; 1218 imid = sctx.imid; 1219 iside = sctx.iside; 1220 delta = sctx.delta; 1221 itheta = sctx.itheta; 1222 qalloc = sctx.qalloc; 1223 #ifdef FIXED_POINT 1224 mid = imid; 1225 side = iside; 1226 #else 1227 mid = (1.f/32768)*imid; 1228 side = (1.f/32768)*iside; 1229 #endif 1230 1231 /* This is a special case for N=2 that only works for stereo and takes 1232 advantage of the fact that mid and side are orthogonal to encode 1233 the side with just one bit. */ 1234 if (N==2) 1235 { 1236 int c; 1237 int sign=0; 1238 celt_norm *x2, *y2; 1239 mbits = b; 1240 sbits = 0; 1241 /* Only need one bit for the side. */ 1242 if (itheta != 0 && itheta != 16384) 1243 sbits = 1<<BITRES; 1244 mbits -= sbits; 1245 c = itheta > 8192; 1246 ctx->remaining_bits -= qalloc+sbits; 1247 1248 x2 = c ? Y : X; 1249 y2 = c ? X : Y; 1250 if (sbits) 1251 { 1252 if (encode) 1253 { 1254 /* Here we only need to encode a sign for the side. */ 1255 sign = x2[0]*y2[1] - x2[1]*y2[0] < 0; 1256 ec_enc_bits(ec, sign, 1); 1257 } else { 1258 sign = ec_dec_bits(ec, 1); 1259 } 1260 } 1261 sign = 1-2*sign; 1262 /* We use orig_fill here because we want to fold the side, but if 1263 itheta==16384, we'll have cleared the low bits of fill. */ 1264 cm = quant_band(ctx, x2, N, mbits, B, lowband, 1265 LM, lowband_out, Q15ONE, lowband_scratch, orig_fill); 1266 /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse), 1267 and there's no need to worry about mixing with the other channel. */ 1268 y2[0] = -sign*x2[1]; 1269 y2[1] = sign*x2[0]; 1270 if (resynth) 1271 { 1272 celt_norm tmp; 1273 X[0] = MULT16_16_Q15(mid, X[0]); 1274 X[1] = MULT16_16_Q15(mid, X[1]); 1275 Y[0] = MULT16_16_Q15(side, Y[0]); 1276 Y[1] = MULT16_16_Q15(side, Y[1]); 1277 tmp = X[0]; 1278 X[0] = SUB16(tmp,Y[0]); 1279 Y[0] = ADD16(tmp,Y[0]); 1280 tmp = X[1]; 1281 X[1] = SUB16(tmp,Y[1]); 1282 Y[1] = ADD16(tmp,Y[1]); 1283 } 1284 } else { 1285 /* "Normal" split code */ 1286 opus_int32 rebalance; 1287 1288 mbits = IMAX(0, IMIN(b, (b-delta)/2)); 1289 sbits = b-mbits; 1290 ctx->remaining_bits -= qalloc; 1291 1292 rebalance = ctx->remaining_bits; 1293 if (mbits >= sbits) 1294 { 1295 /* In stereo mode, we do not apply a scaling to the mid because we need the normalized 1296 mid for folding later. */ 1297 cm = quant_band(ctx, X, N, mbits, B, 1298 lowband, LM, lowband_out, 1299 Q15ONE, lowband_scratch, fill); 1300 rebalance = mbits - (rebalance-ctx->remaining_bits); 1301 if (rebalance > 3<<BITRES && itheta!=0) 1302 sbits += rebalance - (3<<BITRES); 1303 1304 /* For a stereo split, the high bits of fill are always zero, so no 1305 folding will be done to the side. */ 1306 cm |= quant_band(ctx, Y, N, sbits, B, 1307 NULL, LM, NULL, 1308 side, NULL, fill>>B); 1309 } else { 1310 /* For a stereo split, the high bits of fill are always zero, so no 1311 folding will be done to the side. */ 1312 cm = quant_band(ctx, Y, N, sbits, B, 1313 NULL, LM, NULL, 1314 side, NULL, fill>>B); 1315 rebalance = sbits - (rebalance-ctx->remaining_bits); 1316 if (rebalance > 3<<BITRES && itheta!=16384) 1317 mbits += rebalance - (3<<BITRES); 1318 /* In stereo mode, we do not apply a scaling to the mid because we need the normalized 1319 mid for folding later. */ 1320 cm |= quant_band(ctx, X, N, mbits, B, 1321 lowband, LM, lowband_out, 1322 Q15ONE, lowband_scratch, fill); 1323 } 1324 } 1325 1326 1327 /* This code is used by the decoder and by the resynthesis-enabled encoder */ 1328 if (resynth) 1329 { 1330 if (N!=2) 1331 stereo_merge(X, Y, mid, N); 1332 if (inv) 1333 { 1334 int j; 1335 for (j=0;j<N;j++) 1336 Y[j] = -Y[j]; 1337 } 1338 } 1339 return cm; 1340 } 1341 1342 1343 void quant_all_bands(int encode, const CELTMode *m, int start, int end, 1344 celt_norm *X_, celt_norm *Y_, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses, 1345 int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res, 1346 opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed) 1347 { 1348 int i; 1349 opus_int32 remaining_bits; 1350 const opus_int16 * OPUS_RESTRICT eBands = m->eBands; 1351 celt_norm * OPUS_RESTRICT norm, * OPUS_RESTRICT norm2; 1352 VARDECL(celt_norm, _norm); 1353 celt_norm *lowband_scratch; 1354 int B; 1355 int M; 1356 int lowband_offset; 1357 int update_lowband = 1; 1358 int C = Y_ != NULL ? 2 : 1; 1359 int norm_offset; 1360 #ifdef RESYNTH 1361 int resynth = 1; 1362 #else 1363 int resynth = !encode; 1364 #endif 1365 struct band_ctx ctx; 1366 SAVE_STACK; 1367 1368 M = 1<<LM; 1369 B = shortBlocks ? M : 1; 1370 norm_offset = M*eBands[start]; 1371 /* No need to allocate norm for the last band because we don't need an 1372 output in that band. */ 1373 ALLOC(_norm, C*(M*eBands[m->nbEBands-1]-norm_offset), celt_norm); 1374 norm = _norm; 1375 norm2 = norm + M*eBands[m->nbEBands-1]-norm_offset; 1376 /* We can use the last band as scratch space because we don't need that 1377 scratch space for the last band. */ 1378 lowband_scratch = X_+M*eBands[m->nbEBands-1]; 1379 1380 lowband_offset = 0; 1381 ctx.bandE = bandE; 1382 ctx.ec = ec; 1383 ctx.encode = encode; 1384 ctx.intensity = intensity; 1385 ctx.m = m; 1386 ctx.seed = *seed; 1387 ctx.spread = spread; 1388 for (i=start;i<end;i++) 1389 { 1390 opus_int32 tell; 1391 int b; 1392 int N; 1393 opus_int32 curr_balance; 1394 int effective_lowband=-1; 1395 celt_norm * OPUS_RESTRICT X, * OPUS_RESTRICT Y; 1396 int tf_change=0; 1397 unsigned x_cm; 1398 unsigned y_cm; 1399 int last; 1400 1401 ctx.i = i; 1402 last = (i==end-1); 1403 1404 X = X_+M*eBands[i]; 1405 if (Y_!=NULL) 1406 Y = Y_+M*eBands[i]; 1407 else 1408 Y = NULL; 1409 N = M*eBands[i+1]-M*eBands[i]; 1410 tell = ec_tell_frac(ec); 1411 1412 /* Compute how many bits we want to allocate to this band */ 1413 if (i != start) 1414 balance -= tell; 1415 remaining_bits = total_bits-tell-1; 1416 ctx.remaining_bits = remaining_bits; 1417 if (i <= codedBands-1) 1418 { 1419 curr_balance = balance / IMIN(3, codedBands-i); 1420 b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance))); 1421 } else { 1422 b = 0; 1423 } 1424 1425 if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0)) 1426 lowband_offset = i; 1427 1428 tf_change = tf_res[i]; 1429 ctx.tf_change = tf_change; 1430 if (i>=m->effEBands) 1431 { 1432 X=norm; 1433 if (Y_!=NULL) 1434 Y = norm; 1435 lowband_scratch = NULL; 1436 } 1437 if (i==end-1) 1438 lowband_scratch = NULL; 1439 1440 /* Get a conservative estimate of the collapse_mask's for the bands we're 1441 going to be folding from. */ 1442 if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0)) 1443 { 1444 int fold_start; 1445 int fold_end; 1446 int fold_i; 1447 /* This ensures we never repeat spectral content within one band */ 1448 effective_lowband = IMAX(0, M*eBands[lowband_offset]-norm_offset-N); 1449 fold_start = lowband_offset; 1450 while(M*eBands[--fold_start] > effective_lowband+norm_offset); 1451 fold_end = lowband_offset-1; 1452 while(M*eBands[++fold_end] < effective_lowband+norm_offset+N); 1453 x_cm = y_cm = 0; 1454 fold_i = fold_start; do { 1455 x_cm |= collapse_masks[fold_i*C+0]; 1456 y_cm |= collapse_masks[fold_i*C+C-1]; 1457 } while (++fold_i<fold_end); 1458 } 1459 /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost 1460 always) be non-zero. */ 1461 else 1462 x_cm = y_cm = (1<<B)-1; 1463 1464 if (dual_stereo && i==intensity) 1465 { 1466 int j; 1467 1468 /* Switch off dual stereo to do intensity. */ 1469 dual_stereo = 0; 1470 if (resynth) 1471 for (j=0;j<M*eBands[i]-norm_offset;j++) 1472 norm[j] = HALF32(norm[j]+norm2[j]); 1473 } 1474 if (dual_stereo) 1475 { 1476 x_cm = quant_band(&ctx, X, N, b/2, B, 1477 effective_lowband != -1 ? norm+effective_lowband : NULL, LM, 1478 last?NULL:norm+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, x_cm); 1479 y_cm = quant_band(&ctx, Y, N, b/2, B, 1480 effective_lowband != -1 ? norm2+effective_lowband : NULL, LM, 1481 last?NULL:norm2+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, y_cm); 1482 } else { 1483 if (Y!=NULL) 1484 { 1485 x_cm = quant_band_stereo(&ctx, X, Y, N, b, B, 1486 effective_lowband != -1 ? norm+effective_lowband : NULL, LM, 1487 last?NULL:norm+M*eBands[i]-norm_offset, lowband_scratch, x_cm|y_cm); 1488 } else { 1489 x_cm = quant_band(&ctx, X, N, b, B, 1490 effective_lowband != -1 ? norm+effective_lowband : NULL, LM, 1491 last?NULL:norm+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, x_cm|y_cm); 1492 } 1493 y_cm = x_cm; 1494 } 1495 collapse_masks[i*C+0] = (unsigned char)x_cm; 1496 collapse_masks[i*C+C-1] = (unsigned char)y_cm; 1497 balance += pulses[i] + tell; 1498 1499 /* Update the folding position only as long as we have 1 bit/sample depth. */ 1500 update_lowband = b>(N<<BITRES); 1501 } 1502 *seed = ctx.seed; 1503 1504 RESTORE_STACK; 1505 } 1506 1507
