bands.c (61074B)
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 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_sig_assert(tmp<=32767); 74 x2 = tmp; 75 x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2))))); 76 celt_sig_assert(x2<=32766); 77 return 1+x2; 78 } 79 80 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 LM, int arch) 96 { 97 int i, c, N; 98 const opus_int16 *eBands = m->eBands; 99 (void)arch; 100 N = m->shortMdctSize<<LM; 101 c=0; do { 102 for (i=0;i<end;i++) 103 { 104 int j; 105 opus_val32 maxval=0; 106 opus_val32 sum = 0; 107 108 maxval = celt_maxabs32(&X[c*N+(eBands[i]<<LM)], (eBands[i+1]-eBands[i])<<LM); 109 if (maxval > 0) 110 { 111 int shift = IMAX(0, 30 - celt_ilog2(maxval+(maxval>>14)+1) - ((((m->logN[i]+7)>>BITRES)+LM+1)>>1)); 112 j=eBands[i]<<LM; do { 113 opus_val32 x = SHL32(X[j+c*N],shift); 114 sum = ADD32(sum, MULT32_32_Q31(x, x)); 115 } while (++j<eBands[i+1]<<LM); 116 bandE[i+c*m->nbEBands] = MAX32(maxval, PSHR32(celt_sqrt32(SHR32(sum,1)), shift)); 117 } else { 118 bandE[i+c*m->nbEBands] = EPSILON; 119 } 120 } 121 } while (++c<C); 122 } 123 124 /* Normalise each band such that the energy is one. */ 125 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) 126 { 127 int i, c, N; 128 const opus_int16 *eBands = m->eBands; 129 N = M*m->shortMdctSize; 130 c=0; do { 131 i=0; do { 132 int j,shift; 133 opus_val32 E; 134 opus_val32 g; 135 E = bandE[i+c*m->nbEBands]; 136 /* For very low energies, we need this to make sure not to prevent energy rounding from 137 blowing up the normalized signal. */ 138 if (E < 10) E += EPSILON; 139 shift = 30-celt_zlog2(E); 140 E = SHL32(E, shift); 141 g = celt_rcp_norm32(E); 142 j=M*eBands[i]; do { 143 X[j+c*N] = PSHR32(MULT32_32_Q31(g, SHL32(freq[j+c*N], shift)), 30-NORM_SHIFT); 144 } while (++j<M*eBands[i+1]); 145 } while (++i<end); 146 } while (++c<C); 147 } 148 149 #else /* FIXED_POINT */ 150 /* Compute the amplitude (sqrt energy) in each of the bands */ 151 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int LM, int arch) 152 { 153 int i, c, N; 154 const opus_int16 *eBands = m->eBands; 155 N = m->shortMdctSize<<LM; 156 c=0; do { 157 for (i=0;i<end;i++) 158 { 159 opus_val32 sum; 160 sum = 1e-27f + celt_inner_prod(&X[c*N+(eBands[i]<<LM)], &X[c*N+(eBands[i]<<LM)], (eBands[i+1]-eBands[i])<<LM, arch); 161 bandE[i+c*m->nbEBands] = celt_sqrt(sum); 162 /*printf ("%f ", bandE[i+c*m->nbEBands]);*/ 163 } 164 } while (++c<C); 165 /*printf ("\n");*/ 166 } 167 168 /* Normalise each band such that the energy is one. */ 169 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) 170 { 171 int i, c, N; 172 const opus_int16 *eBands = m->eBands; 173 N = M*m->shortMdctSize; 174 c=0; do { 175 for (i=0;i<end;i++) 176 { 177 int j; 178 opus_val16 g = 1.f/(1e-27f+bandE[i+c*m->nbEBands]); 179 for (j=M*eBands[i];j<M*eBands[i+1];j++) 180 X[j+c*N] = freq[j+c*N]*g; 181 } 182 } while (++c<C); 183 } 184 185 #endif /* FIXED_POINT */ 186 187 /* De-normalise the energy to produce the synthesis from the unit-energy bands */ 188 void denormalise_bands(const CELTMode *m, const celt_norm * OPUS_RESTRICT X, 189 celt_sig * OPUS_RESTRICT freq, const celt_glog *bandLogE, int start, 190 int end, int M, int downsample, int silence) 191 { 192 int i, N; 193 int bound; 194 celt_sig * OPUS_RESTRICT f; 195 const celt_norm * OPUS_RESTRICT x; 196 const opus_int16 *eBands = m->eBands; 197 N = M*m->shortMdctSize; 198 bound = M*eBands[end]; 199 if (downsample!=1) 200 bound = IMIN(bound, N/downsample); 201 if (silence) 202 { 203 bound = 0; 204 start = end = 0; 205 } 206 f = freq; 207 x = X+M*eBands[start]; 208 if (start != 0) 209 { 210 for (i=0;i<M*eBands[start];i++) 211 *f++ = 0; 212 } else { 213 f += M*eBands[start]; 214 } 215 for (i=start;i<end;i++) 216 { 217 int j, band_end; 218 opus_val32 g; 219 celt_glog lg; 220 #ifdef FIXED_POINT 221 int shift; 222 #endif 223 j=M*eBands[i]; 224 band_end = M*eBands[i+1]; 225 lg = ADD32(bandLogE[i], SHL32((opus_val32)eMeans[i],DB_SHIFT-4)); 226 #ifndef FIXED_POINT 227 g = celt_exp2_db(MIN32(32.f, lg)); 228 #else 229 /* Handle the integer part of the log energy */ 230 shift = 17-(lg>>DB_SHIFT); 231 if (shift>=31) 232 { 233 shift=0; 234 g=0; 235 } else { 236 /* Handle the fractional part. */ 237 g = SHL32(celt_exp2_db_frac((lg&((1<<DB_SHIFT)-1))), 2); 238 } 239 /* Handle extreme gains with negative shift. */ 240 if (shift<0) 241 { 242 /* To avoid overflow, we're 243 capping the gain here, which is equivalent to a cap of 18 on lg. 244 This shouldn't trigger unless the bitstream is already corrupted. */ 245 g = 2147483647; 246 shift = 0; 247 } 248 #endif 249 do { 250 *f++ = PSHR32(MULT32_32_Q31(SHL32(*x, 30-NORM_SHIFT), g), shift); 251 x++; 252 } while (++j<band_end); 253 } 254 celt_assert(start <= end); 255 OPUS_CLEAR(&freq[bound], N-bound); 256 } 257 258 /* This prevents energy collapse for transients with multiple short MDCTs */ 259 void anti_collapse(const CELTMode *m, celt_norm *X_, unsigned char *collapse_masks, int LM, int C, int size, 260 int start, int end, const celt_glog *logE, const celt_glog *prev1logE, 261 const celt_glog *prev2logE, const int *pulses, opus_uint32 seed, int encode, int arch) 262 { 263 int c, i, j, k; 264 for (i=start;i<end;i++) 265 { 266 int N0; 267 opus_val16 thresh, sqrt_1; 268 int depth; 269 #ifdef FIXED_POINT 270 int shift; 271 opus_val32 thresh32; 272 #endif 273 274 N0 = m->eBands[i+1]-m->eBands[i]; 275 /* depth in 1/8 bits */ 276 celt_sig_assert(pulses[i]>=0); 277 depth = celt_udiv(1+pulses[i], (m->eBands[i+1]-m->eBands[i]))>>LM; 278 279 #ifdef FIXED_POINT 280 thresh32 = SHR32(celt_exp2(-SHL16(depth, 10-BITRES)),1); 281 thresh = MULT16_32_Q15(QCONST16(0.5f, 15), MIN32(32767,thresh32)); 282 { 283 opus_val32 t; 284 t = N0<<LM; 285 shift = celt_ilog2(t)>>1; 286 t = SHL32(t, (7-shift)<<1); 287 sqrt_1 = celt_rsqrt_norm(t); 288 } 289 #else 290 thresh = .5f*celt_exp2(-.125f*depth); 291 sqrt_1 = celt_rsqrt(N0<<LM); 292 #endif 293 294 c=0; do 295 { 296 celt_norm *X; 297 celt_glog prev1; 298 celt_glog prev2; 299 opus_val32 Ediff; 300 celt_norm r; 301 int renormalize=0; 302 prev1 = prev1logE[c*m->nbEBands+i]; 303 prev2 = prev2logE[c*m->nbEBands+i]; 304 if (!encode && C==1) 305 { 306 prev1 = MAXG(prev1,prev1logE[m->nbEBands+i]); 307 prev2 = MAXG(prev2,prev2logE[m->nbEBands+i]); 308 } 309 Ediff = logE[c*m->nbEBands+i]-MING(prev1,prev2); 310 Ediff = MAX32(0, Ediff); 311 312 #ifdef FIXED_POINT 313 if (Ediff < GCONST(16.f)) 314 { 315 opus_val32 r32 = SHR32(celt_exp2_db(-Ediff),1); 316 r = 2*MIN16(16383,r32); 317 } else { 318 r = 0; 319 } 320 if (LM==3) 321 r = MULT16_16_Q14(23170, MIN32(23169, r)); 322 r = SHR16(MIN16(thresh, r),1); 323 r = VSHR32(MULT16_16_Q15(sqrt_1, r),shift+14-NORM_SHIFT); 324 #else 325 /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because 326 short blocks don't have the same energy as long */ 327 r = 2.f*celt_exp2_db(-Ediff); 328 if (LM==3) 329 r *= 1.41421356f; 330 r = MIN16(thresh, r); 331 r = r*sqrt_1; 332 #endif 333 X = X_+c*size+(m->eBands[i]<<LM); 334 for (k=0;k<1<<LM;k++) 335 { 336 /* Detect collapse */ 337 if (!(collapse_masks[i*C+c]&1<<k)) 338 { 339 /* Fill with noise */ 340 for (j=0;j<N0;j++) 341 { 342 seed = celt_lcg_rand(seed); 343 X[(j<<LM)+k] = (seed&0x8000 ? r : -r); 344 } 345 renormalize = 1; 346 } 347 } 348 /* We just added some energy, so we need to renormalise */ 349 if (renormalize) 350 renormalise_vector(X, N0<<LM, Q31ONE, arch); 351 } while (++c<C); 352 } 353 } 354 355 /* Compute the weights to use for optimizing normalized distortion across 356 channels. We use the amplitude to weight square distortion, which means 357 that we use the square root of the value we would have been using if we 358 wanted to minimize the MSE in the non-normalized domain. This roughly 359 corresponds to some quick-and-dirty perceptual experiments I ran to 360 measure inter-aural masking (there doesn't seem to be any published data 361 on the topic). */ 362 static void compute_channel_weights(celt_ener Ex, celt_ener Ey, opus_val16 w[2]) 363 { 364 celt_ener minE; 365 #ifdef FIXED_POINT 366 int shift; 367 #endif 368 minE = MIN32(Ex, Ey); 369 /* Adjustment to make the weights a bit more conservative. */ 370 Ex = ADD32(Ex, minE/3); 371 Ey = ADD32(Ey, minE/3); 372 #ifdef FIXED_POINT 373 shift = celt_ilog2(EPSILON+MAX32(Ex, Ey))-14; 374 #endif 375 w[0] = VSHR32(Ex, shift); 376 w[1] = VSHR32(Ey, shift); 377 } 378 379 static void intensity_stereo(const CELTMode *m, celt_norm * OPUS_RESTRICT X, const celt_norm * OPUS_RESTRICT Y, const celt_ener *bandE, int bandID, int N) 380 { 381 int i = bandID; 382 int j; 383 opus_val16 a1, a2; 384 opus_val16 left, right; 385 opus_val16 norm; 386 #ifdef FIXED_POINT 387 int shift = celt_zlog2(MAX32(bandE[i], bandE[i+m->nbEBands]))-13; 388 #endif 389 left = VSHR32(bandE[i],shift); 390 right = VSHR32(bandE[i+m->nbEBands],shift); 391 norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right)); 392 #ifdef FIXED_POINT 393 left = MIN32(left, norm-1); 394 right = MIN32(right, norm-1); 395 #endif 396 a1 = DIV32_16(SHL32(EXTEND32(left),15),norm); 397 a2 = DIV32_16(SHL32(EXTEND32(right),15),norm); 398 for (j=0;j<N;j++) 399 { 400 X[j] = ADD32(MULT16_32_Q15(a1, X[j]), MULT16_32_Q15(a2, Y[j])); 401 /* Side is not encoded, no need to calculate */ 402 } 403 } 404 405 static void stereo_split(celt_norm * OPUS_RESTRICT X, celt_norm * OPUS_RESTRICT Y, int N) 406 { 407 int j; 408 for (j=0;j<N;j++) 409 { 410 opus_val32 r, l; 411 l = MULT32_32_Q31(QCONST32(.70710678f,31), X[j]); 412 r = MULT32_32_Q31(QCONST32(.70710678f,31), Y[j]); 413 X[j] = ADD32(l, r); 414 Y[j] = SUB32(r, l); 415 } 416 } 417 418 static void stereo_merge(celt_norm * OPUS_RESTRICT X, celt_norm * OPUS_RESTRICT Y, opus_val32 mid, int N, int arch) 419 { 420 int j; 421 opus_val32 xp=0, side=0; 422 opus_val32 El, Er; 423 #ifdef FIXED_POINT 424 int kl, kr; 425 #endif 426 opus_val32 t, lgain, rgain; 427 428 /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */ 429 xp = celt_inner_prod_norm_shift(Y, X, N, arch); 430 side = celt_inner_prod_norm_shift(Y, Y, N, arch); 431 /* Compensating for the mid normalization */ 432 xp = MULT32_32_Q31(mid, xp); 433 /* mid and side are in Q15, not Q14 like X and Y */ 434 El = SHR32(MULT32_32_Q31(mid, mid),3) + side - 2*xp; 435 Er = SHR32(MULT32_32_Q31(mid, mid),3) + side + 2*xp; 436 if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28)) 437 { 438 OPUS_COPY(Y, X, N); 439 return; 440 } 441 442 #ifdef FIXED_POINT 443 kl = celt_ilog2(El)>>1; 444 kr = celt_ilog2(Er)>>1; 445 #endif 446 t = VSHR32(El, (kl<<1)-29); 447 lgain = celt_rsqrt_norm32(t); 448 t = VSHR32(Er, (kr<<1)-29); 449 rgain = celt_rsqrt_norm32(t); 450 451 #ifdef FIXED_POINT 452 if (kl < 7) 453 kl = 7; 454 if (kr < 7) 455 kr = 7; 456 #endif 457 458 for (j=0;j<N;j++) 459 { 460 celt_norm r, l; 461 /* Apply mid scaling (side is already scaled) */ 462 l = MULT32_32_Q31(mid, X[j]); 463 r = Y[j]; 464 X[j] = VSHR32(MULT32_32_Q31(lgain, SUB32(l,r)), kl-15); 465 Y[j] = VSHR32(MULT32_32_Q31(rgain, ADD32(l,r)), kr-15); 466 } 467 } 468 469 /* Decide whether we should spread the pulses in the current frame */ 470 int spreading_decision(const CELTMode *m, const celt_norm *X, int *average, 471 int last_decision, int *hf_average, int *tapset_decision, int update_hf, 472 int end, int C, int M, const int *spread_weight) 473 { 474 int i, c, N0; 475 int sum = 0, nbBands=0; 476 const opus_int16 * OPUS_RESTRICT eBands = m->eBands; 477 int decision; 478 int hf_sum=0; 479 480 celt_assert(end>0); 481 482 N0 = M*m->shortMdctSize; 483 484 if (M*(eBands[end]-eBands[end-1]) <= 8) 485 return SPREAD_NONE; 486 c=0; do { 487 for (i=0;i<end;i++) 488 { 489 int j, N, tmp=0; 490 int tcount[3] = {0,0,0}; 491 const celt_norm * OPUS_RESTRICT x = X+M*eBands[i]+c*N0; 492 N = M*(eBands[i+1]-eBands[i]); 493 if (N<=8) 494 continue; 495 /* Compute rough CDF of |x[j]| */ 496 for (j=0;j<N;j++) 497 { 498 opus_val32 x2N; /* Q13 */ 499 500 x2N = MULT16_16(MULT16_16_Q15(SHR32(x[j], NORM_SHIFT-14), SHR32(x[j], NORM_SHIFT-14)), N); 501 if (x2N < QCONST16(0.25f,13)) 502 tcount[0]++; 503 if (x2N < QCONST16(0.0625f,13)) 504 tcount[1]++; 505 if (x2N < QCONST16(0.015625f,13)) 506 tcount[2]++; 507 } 508 509 /* Only include four last bands (8 kHz and up) */ 510 if (i>m->nbEBands-4) 511 hf_sum += celt_udiv(32*(tcount[1]+tcount[0]), N); 512 tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N); 513 sum += tmp*spread_weight[i]; 514 nbBands+=spread_weight[i]; 515 } 516 } while (++c<C); 517 518 if (update_hf) 519 { 520 if (hf_sum) 521 hf_sum = celt_udiv(hf_sum, C*(4-m->nbEBands+end)); 522 *hf_average = (*hf_average+hf_sum)>>1; 523 hf_sum = *hf_average; 524 if (*tapset_decision==2) 525 hf_sum += 4; 526 else if (*tapset_decision==0) 527 hf_sum -= 4; 528 if (hf_sum > 22) 529 *tapset_decision=2; 530 else if (hf_sum > 18) 531 *tapset_decision=1; 532 else 533 *tapset_decision=0; 534 } 535 /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/ 536 celt_assert(nbBands>0); /* end has to be non-zero */ 537 celt_assert(sum>=0); 538 sum = celt_udiv((opus_int32)sum<<8, nbBands); 539 /* Recursive averaging */ 540 sum = (sum+*average)>>1; 541 *average = sum; 542 /* Hysteresis */ 543 sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2; 544 if (sum < 80) 545 { 546 decision = SPREAD_AGGRESSIVE; 547 } else if (sum < 256) 548 { 549 decision = SPREAD_NORMAL; 550 } else if (sum < 384) 551 { 552 decision = SPREAD_LIGHT; 553 } else { 554 decision = SPREAD_NONE; 555 } 556 #ifdef FUZZING 557 decision = rand()&0x3; 558 *tapset_decision=rand()%3; 559 #endif 560 return decision; 561 } 562 563 /* Indexing table for converting from natural Hadamard to ordery Hadamard 564 This is essentially a bit-reversed Gray, on top of which we've added 565 an inversion of the order because we want the DC at the end rather than 566 the beginning. The lines are for N=2, 4, 8, 16 */ 567 static const int ordery_table[] = { 568 1, 0, 569 3, 0, 2, 1, 570 7, 0, 4, 3, 6, 1, 5, 2, 571 15, 0, 8, 7, 12, 3, 11, 4, 14, 1, 9, 6, 13, 2, 10, 5, 572 }; 573 574 static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard) 575 { 576 int i,j; 577 VARDECL(celt_norm, tmp); 578 int N; 579 SAVE_STACK; 580 N = N0*stride; 581 ALLOC(tmp, N, celt_norm); 582 celt_assert(stride>0); 583 if (hadamard) 584 { 585 const int *ordery = ordery_table+stride-2; 586 for (i=0;i<stride;i++) 587 { 588 for (j=0;j<N0;j++) 589 tmp[ordery[i]*N0+j] = X[j*stride+i]; 590 } 591 } else { 592 for (i=0;i<stride;i++) 593 for (j=0;j<N0;j++) 594 tmp[i*N0+j] = X[j*stride+i]; 595 } 596 OPUS_COPY(X, tmp, N); 597 RESTORE_STACK; 598 } 599 600 static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard) 601 { 602 int i,j; 603 VARDECL(celt_norm, tmp); 604 int N; 605 SAVE_STACK; 606 N = N0*stride; 607 ALLOC(tmp, N, celt_norm); 608 if (hadamard) 609 { 610 const int *ordery = ordery_table+stride-2; 611 for (i=0;i<stride;i++) 612 for (j=0;j<N0;j++) 613 tmp[j*stride+i] = X[ordery[i]*N0+j]; 614 } else { 615 for (i=0;i<stride;i++) 616 for (j=0;j<N0;j++) 617 tmp[j*stride+i] = X[i*N0+j]; 618 } 619 OPUS_COPY(X, tmp, N); 620 RESTORE_STACK; 621 } 622 623 void haar1(celt_norm *X, int N0, int stride) 624 { 625 int i, j; 626 N0 >>= 1; 627 for (i=0;i<stride;i++) 628 for (j=0;j<N0;j++) 629 { 630 opus_val32 tmp1, tmp2; 631 tmp1 = MULT32_32_Q31(QCONST32(.70710678f,31), X[stride*2*j+i]); 632 tmp2 = MULT32_32_Q31(QCONST32(.70710678f,31), X[stride*(2*j+1)+i]); 633 X[stride*2*j+i] = ADD32(tmp1, tmp2); 634 X[stride*(2*j+1)+i] = SUB32(tmp1, tmp2); 635 } 636 } 637 638 static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo) 639 { 640 static const opus_int16 exp2_table8[8] = 641 {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048}; 642 int qn, qb; 643 int N2 = 2*N-1; 644 if (stereo && N==2) 645 N2--; 646 /* The upper limit ensures that in a stereo split with itheta==16384, we'll 647 always have enough bits left over to code at least one pulse in the 648 side; otherwise it would collapse, since it doesn't get folded. */ 649 qb = celt_sudiv(b+N2*offset, N2); 650 qb = IMIN(b-pulse_cap-(4<<BITRES), qb); 651 652 qb = IMIN(8<<BITRES, qb); 653 654 if (qb<(1<<BITRES>>1)) { 655 qn = 1; 656 } else { 657 qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES)); 658 qn = (qn+1)>>1<<1; 659 } 660 celt_assert(qn <= 256); 661 return qn; 662 } 663 664 struct band_ctx { 665 int encode; 666 int resynth; 667 const CELTMode *m; 668 int i; 669 int intensity; 670 int spread; 671 int tf_change; 672 ec_ctx *ec; 673 opus_int32 remaining_bits; 674 const celt_ener *bandE; 675 opus_uint32 seed; 676 int arch; 677 int theta_round; 678 int disable_inv; 679 int avoid_split_noise; 680 #ifdef ENABLE_QEXT 681 ec_ctx *ext_ec; 682 int extra_bits; 683 opus_int32 ext_total_bits; 684 int extra_bands; 685 #endif 686 }; 687 688 struct split_ctx { 689 int inv; 690 int imid; 691 int iside; 692 int delta; 693 int itheta; 694 #ifdef ENABLE_QEXT 695 int itheta_q30; 696 #endif 697 int qalloc; 698 }; 699 700 static void compute_theta(struct band_ctx *ctx, struct split_ctx *sctx, 701 celt_norm *X, celt_norm *Y, int N, int *b, int B, int B0, 702 int LM, 703 int stereo, int *fill ARG_QEXT(int *ext_b)) 704 { 705 int qn; 706 int itheta=0; 707 int itheta_q30=0; 708 int delta; 709 int imid, iside; 710 int qalloc; 711 int pulse_cap; 712 int offset; 713 opus_int32 tell; 714 int inv=0; 715 int encode; 716 const CELTMode *m; 717 int i; 718 int intensity; 719 ec_ctx *ec; 720 const celt_ener *bandE; 721 722 encode = ctx->encode; 723 m = ctx->m; 724 i = ctx->i; 725 intensity = ctx->intensity; 726 ec = ctx->ec; 727 bandE = ctx->bandE; 728 729 /* Decide on the resolution to give to the split parameter theta */ 730 pulse_cap = m->logN[i]+LM*(1<<BITRES); 731 offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET); 732 qn = compute_qn(N, *b, offset, pulse_cap, stereo); 733 if (stereo && i>=intensity) 734 qn = 1; 735 if (encode) 736 { 737 /* theta is the atan() of the ratio between the (normalized) 738 side and mid. With just that parameter, we can re-scale both 739 mid and side because we know that 1) they have unit norm and 740 2) they are orthogonal. */ 741 itheta_q30 = stereo_itheta(X, Y, stereo, N, ctx->arch); 742 itheta = itheta_q30>>16; 743 } 744 tell = ec_tell_frac(ec); 745 if (qn!=1) 746 { 747 if (encode) 748 { 749 if (!stereo || ctx->theta_round == 0) 750 { 751 itheta = (itheta*(opus_int32)qn+8192)>>14; 752 if (!stereo && ctx->avoid_split_noise && itheta > 0 && itheta < qn) 753 { 754 /* Check if the selected value of theta will cause the bit allocation 755 to inject noise on one side. If so, make sure the energy of that side 756 is zero. */ 757 int unquantized = celt_udiv((opus_int32)itheta*16384, qn); 758 imid = bitexact_cos((opus_int16)unquantized); 759 iside = bitexact_cos((opus_int16)(16384-unquantized)); 760 delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid)); 761 if (delta > *b) 762 itheta = qn; 763 else if (delta < -*b) 764 itheta = 0; 765 } 766 } else { 767 int down; 768 /* Bias quantization towards itheta=0 and itheta=16384. */ 769 int bias = itheta > 8192 ? 32767/qn : -32767/qn; 770 down = IMIN(qn-1, IMAX(0, (itheta*(opus_int32)qn + bias)>>14)); 771 if (ctx->theta_round < 0) 772 itheta = down; 773 else 774 itheta = down+1; 775 } 776 } 777 /* Entropy coding of the angle. We use a uniform pdf for the 778 time split, a step for stereo, and a triangular one for the rest. */ 779 if (stereo && N>2) 780 { 781 int p0 = 3; 782 int x = itheta; 783 int x0 = qn/2; 784 int ft = p0*(x0+1) + x0; 785 /* Use a probability of p0 up to itheta=8192 and then use 1 after */ 786 if (encode) 787 { 788 ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft); 789 } else { 790 int fs; 791 fs=ec_decode(ec,ft); 792 if (fs<(x0+1)*p0) 793 x=fs/p0; 794 else 795 x=x0+1+(fs-(x0+1)*p0); 796 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); 797 itheta = x; 798 } 799 } else if (B0>1 || stereo) { 800 /* Uniform pdf */ 801 if (encode) 802 ec_enc_uint(ec, itheta, qn+1); 803 else 804 itheta = ec_dec_uint(ec, qn+1); 805 } else { 806 int fs=1, ft; 807 ft = ((qn>>1)+1)*((qn>>1)+1); 808 if (encode) 809 { 810 int fl; 811 812 fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta; 813 fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 : 814 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1); 815 816 ec_encode(ec, fl, fl+fs, ft); 817 } else { 818 /* Triangular pdf */ 819 int fl=0; 820 int fm; 821 fm = ec_decode(ec, ft); 822 823 if (fm < ((qn>>1)*((qn>>1) + 1)>>1)) 824 { 825 itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1; 826 fs = itheta + 1; 827 fl = itheta*(itheta + 1)>>1; 828 } 829 else 830 { 831 itheta = (2*(qn + 1) 832 - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1; 833 fs = qn + 1 - itheta; 834 fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1); 835 } 836 837 ec_dec_update(ec, fl, fl+fs, ft); 838 } 839 } 840 celt_assert(itheta>=0); 841 itheta = celt_udiv((opus_int32)itheta*16384, qn); 842 #ifdef ENABLE_QEXT 843 *ext_b = IMIN(*ext_b, ctx->ext_total_bits - (opus_int32)ec_tell_frac(ctx->ext_ec)); 844 if (*ext_b >= 2*N<<BITRES && ctx->ext_total_bits-ec_tell_frac(ctx->ext_ec)-1 > 2<<BITRES) { 845 int extra_bits; 846 int ext_tell = ec_tell_frac(ctx->ext_ec); 847 extra_bits = IMIN(12, IMAX(2, celt_sudiv(*ext_b, (2*N-1)<<BITRES))); 848 if (encode) { 849 itheta_q30 = itheta_q30 - (itheta<<16); 850 itheta_q30 = (itheta_q30*(opus_int64)qn*((1<<extra_bits)-1)+(1<<29))>>30; 851 itheta_q30 += (1<<(extra_bits-1))-1; 852 itheta_q30 = IMAX(0, IMIN((1<<extra_bits)-2, itheta_q30)); 853 ec_enc_uint(ctx->ext_ec, itheta_q30, (1<<extra_bits)-1); 854 } else { 855 itheta_q30 = ec_dec_uint(ctx->ext_ec, (1<<extra_bits)-1); 856 } 857 itheta_q30 -= (1<<(extra_bits-1))-1; 858 itheta_q30 = (itheta<<16) + itheta_q30*(opus_int64)(1<<30)/(qn*((1<<extra_bits)-1)); 859 /* Hard bounds on itheta (can only trigger on corrupted bitstreams). */ 860 itheta_q30 = IMAX(0, IMIN(itheta_q30, 1073741824)); 861 *ext_b -= ec_tell_frac(ctx->ext_ec) - ext_tell; 862 } else { 863 itheta_q30 = (opus_int32)itheta<<16; 864 } 865 #endif 866 if (encode && stereo) 867 { 868 if (itheta==0) 869 intensity_stereo(m, X, Y, bandE, i, N); 870 else 871 stereo_split(X, Y, N); 872 } 873 /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate. 874 Let's do that at higher complexity */ 875 } else if (stereo) { 876 if (encode) 877 { 878 inv = itheta > 8192 && !ctx->disable_inv; 879 if (inv) 880 { 881 int j; 882 for (j=0;j<N;j++) 883 Y[j] = -Y[j]; 884 } 885 intensity_stereo(m, X, Y, bandE, i, N); 886 } 887 if (*b>2<<BITRES && ctx->remaining_bits > 2<<BITRES) 888 { 889 if (encode) 890 ec_enc_bit_logp(ec, inv, 2); 891 else 892 inv = ec_dec_bit_logp(ec, 2); 893 } else 894 inv = 0; 895 /* inv flag override to avoid problems with downmixing. */ 896 if (ctx->disable_inv) 897 inv = 0; 898 itheta = 0; 899 itheta_q30 = 0; 900 } 901 qalloc = ec_tell_frac(ec) - tell; 902 *b -= qalloc; 903 904 if (itheta == 0) 905 { 906 imid = 32767; 907 iside = 0; 908 *fill &= (1<<B)-1; 909 delta = -16384; 910 } else if (itheta == 16384) 911 { 912 imid = 0; 913 iside = 32767; 914 *fill &= ((1<<B)-1)<<B; 915 delta = 16384; 916 } else { 917 imid = bitexact_cos((opus_int16)itheta); 918 iside = bitexact_cos((opus_int16)(16384-itheta)); 919 /* This is the mid vs side allocation that minimizes squared error 920 in that band. */ 921 delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid)); 922 } 923 924 sctx->inv = inv; 925 sctx->imid = imid; 926 sctx->iside = iside; 927 sctx->delta = delta; 928 sctx->itheta = itheta; 929 #ifdef ENABLE_QEXT 930 sctx->itheta_q30 = itheta_q30; 931 #endif 932 sctx->qalloc = qalloc; 933 } 934 static unsigned quant_band_n1(struct band_ctx *ctx, celt_norm *X, celt_norm *Y, 935 celt_norm *lowband_out) 936 { 937 int c; 938 int stereo; 939 celt_norm *x = X; 940 int encode; 941 ec_ctx *ec; 942 943 encode = ctx->encode; 944 ec = ctx->ec; 945 946 stereo = Y != NULL; 947 c=0; do { 948 int sign=0; 949 if (ctx->remaining_bits>=1<<BITRES) 950 { 951 if (encode) 952 { 953 sign = x[0]<0; 954 ec_enc_bits(ec, sign, 1); 955 } else { 956 sign = ec_dec_bits(ec, 1); 957 } 958 ctx->remaining_bits -= 1<<BITRES; 959 } 960 if (ctx->resynth) 961 x[0] = sign ? -NORM_SCALING : NORM_SCALING; 962 x = Y; 963 } while (++c<1+stereo); 964 if (lowband_out) 965 lowband_out[0] = SHR32(X[0],4); 966 return 1; 967 } 968 969 /* This function is responsible for encoding and decoding a mono partition. 970 It can split the band in two and transmit the energy difference with 971 the two half-bands. It can be called recursively so bands can end up being 972 split in 8 parts. */ 973 static unsigned quant_partition(struct band_ctx *ctx, celt_norm *X, 974 int N, int b, int B, celt_norm *lowband, 975 int LM, 976 opus_val32 gain, int fill 977 ARG_QEXT(int ext_b)) 978 { 979 const unsigned char *cache; 980 int q; 981 int curr_bits; 982 int imid=0, iside=0; 983 int B0=B; 984 opus_val32 mid=0, side=0; 985 unsigned cm=0; 986 celt_norm *Y=NULL; 987 int encode; 988 const CELTMode *m; 989 int i; 990 int spread; 991 ec_ctx *ec; 992 993 encode = ctx->encode; 994 m = ctx->m; 995 i = ctx->i; 996 spread = ctx->spread; 997 ec = ctx->ec; 998 999 /* If we need 1.5 more bit than we can produce, split the band in two. */ 1000 cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i]; 1001 if (LM != -1 && b > cache[cache[0]]+12 && N>2) 1002 { 1003 int mbits, sbits, delta; 1004 int itheta; 1005 int qalloc; 1006 struct split_ctx sctx; 1007 celt_norm *next_lowband2=NULL; 1008 opus_int32 rebalance; 1009 1010 N >>= 1; 1011 Y = X+N; 1012 LM -= 1; 1013 if (B==1) 1014 fill = (fill&1)|(fill<<1); 1015 B = (B+1)>>1; 1016 1017 compute_theta(ctx, &sctx, X, Y, N, &b, B, B0, LM, 0, &fill ARG_QEXT(&ext_b)); 1018 imid = sctx.imid; 1019 iside = sctx.iside; 1020 delta = sctx.delta; 1021 itheta = sctx.itheta; 1022 qalloc = sctx.qalloc; 1023 #ifdef FIXED_POINT 1024 # ifdef ENABLE_QEXT 1025 (void)imid; 1026 (void)iside; 1027 mid = celt_cos_norm32(sctx.itheta_q30); 1028 side = celt_cos_norm32((1<<30)-sctx.itheta_q30); 1029 # else 1030 mid = SHL32(EXTEND32(imid), 16); 1031 side = SHL32(EXTEND32(iside), 16); 1032 # endif 1033 #else 1034 # ifdef ENABLE_QEXT 1035 (void)imid; 1036 (void)iside; 1037 mid = celt_cos_norm2(sctx.itheta_q30*(1.f/(1<<30))); 1038 side = celt_cos_norm2(1.f-sctx.itheta_q30*(1.f/(1<<30))); 1039 # else 1040 mid = (1.f/32768)*imid; 1041 side = (1.f/32768)*iside; 1042 # endif 1043 #endif 1044 1045 /* Give more bits to low-energy MDCTs than they would otherwise deserve */ 1046 if (B0>1 && (itheta&0x3fff)) 1047 { 1048 if (itheta > 8192) 1049 /* Rough approximation for pre-echo masking */ 1050 delta -= delta>>(4-LM); 1051 else 1052 /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */ 1053 delta = IMIN(0, delta + (N<<BITRES>>(5-LM))); 1054 } 1055 mbits = IMAX(0, IMIN(b, (b-delta)/2)); 1056 sbits = b-mbits; 1057 ctx->remaining_bits -= qalloc; 1058 1059 if (lowband) 1060 next_lowband2 = lowband+N; /* >32-bit split case */ 1061 1062 rebalance = ctx->remaining_bits; 1063 if (mbits >= sbits) 1064 { 1065 cm = quant_partition(ctx, X, N, mbits, B, lowband, LM, 1066 MULT32_32_Q31(gain,mid), fill ARG_QEXT(ext_b/2)); 1067 rebalance = mbits - (rebalance-ctx->remaining_bits); 1068 if (rebalance > 3<<BITRES && itheta!=0) 1069 sbits += rebalance - (3<<BITRES); 1070 cm |= quant_partition(ctx, Y, N, sbits, B, next_lowband2, LM, 1071 MULT32_32_Q31(gain,side), fill>>B ARG_QEXT(ext_b/2))<<(B0>>1); 1072 } else { 1073 cm = quant_partition(ctx, Y, N, sbits, B, next_lowband2, LM, 1074 MULT32_32_Q31(gain,side), fill>>B ARG_QEXT(ext_b/2))<<(B0>>1); 1075 rebalance = sbits - (rebalance-ctx->remaining_bits); 1076 if (rebalance > 3<<BITRES && itheta!=16384) 1077 mbits += rebalance - (3<<BITRES); 1078 cm |= quant_partition(ctx, X, N, mbits, B, lowband, LM, 1079 MULT32_32_Q31(gain,mid), fill ARG_QEXT(ext_b/2)); 1080 } 1081 } else { 1082 #ifdef ENABLE_QEXT 1083 int extra_bits; 1084 int ext_remaining_bits; 1085 extra_bits = ext_b/(N-1)>>BITRES; 1086 ext_remaining_bits = ctx->ext_total_bits-(opus_int32)ec_tell_frac(ctx->ext_ec); 1087 if (ext_remaining_bits < ((extra_bits+1)*(N-1)+N)<<BITRES) { 1088 extra_bits = (ext_remaining_bits-(N<<BITRES))/(N-1)>>BITRES; 1089 extra_bits = IMAX(extra_bits-1, 0); 1090 } 1091 extra_bits = IMIN(12, extra_bits); 1092 #endif 1093 /* This is the basic no-split case */ 1094 q = bits2pulses(m, i, LM, b); 1095 curr_bits = pulses2bits(m, i, LM, q); 1096 ctx->remaining_bits -= curr_bits; 1097 1098 /* Ensures we can never bust the budget */ 1099 while (ctx->remaining_bits < 0 && q > 0) 1100 { 1101 ctx->remaining_bits += curr_bits; 1102 q--; 1103 curr_bits = pulses2bits(m, i, LM, q); 1104 ctx->remaining_bits -= curr_bits; 1105 } 1106 1107 if (q!=0) 1108 { 1109 int K = get_pulses(q); 1110 1111 /* Finally do the actual quantization */ 1112 if (encode) 1113 { 1114 cm = alg_quant(X, N, K, spread, B, ec, gain, ctx->resynth 1115 ARG_QEXT(ctx->ext_ec) ARG_QEXT(extra_bits), 1116 ctx->arch); 1117 } else { 1118 cm = alg_unquant(X, N, K, spread, B, ec, gain 1119 ARG_QEXT(ctx->ext_ec) ARG_QEXT(extra_bits)); 1120 } 1121 #ifdef ENABLE_QEXT 1122 } else if (ext_b > 2*N<<BITRES) 1123 { 1124 extra_bits = ext_b/(N-1)>>BITRES; 1125 ext_remaining_bits = ctx->ext_total_bits-ec_tell_frac(ctx->ext_ec); 1126 if (ext_remaining_bits < ((extra_bits+1)*(N-1)+N)<<BITRES) { 1127 extra_bits = (ext_remaining_bits-(N<<BITRES))/(N-1)>>BITRES; 1128 extra_bits = IMAX(extra_bits-1, 0); 1129 } 1130 extra_bits = IMIN(14, extra_bits); 1131 if (encode) cm = cubic_quant(X, N, extra_bits, B, ctx->ext_ec, gain, ctx->resynth); 1132 else cm = cubic_unquant(X, N, extra_bits, B, ctx->ext_ec, gain); 1133 #endif 1134 } else { 1135 /* If there's no pulse, fill the band anyway */ 1136 int j; 1137 if (ctx->resynth) 1138 { 1139 unsigned cm_mask; 1140 /* B can be as large as 16, so this shift might overflow an int on a 1141 16-bit platform; use a long to get defined behavior.*/ 1142 cm_mask = (unsigned)(1UL<<B)-1; 1143 fill &= cm_mask; 1144 if (!fill) 1145 { 1146 OPUS_CLEAR(X, N); 1147 } else { 1148 if (lowband == NULL) 1149 { 1150 /* Noise */ 1151 for (j=0;j<N;j++) 1152 { 1153 ctx->seed = celt_lcg_rand(ctx->seed); 1154 X[j] = SHL32((celt_norm)((opus_int32)ctx->seed>>20), NORM_SHIFT-14); 1155 } 1156 cm = cm_mask; 1157 } else { 1158 /* Folded spectrum */ 1159 for (j=0;j<N;j++) 1160 { 1161 opus_val16 tmp; 1162 ctx->seed = celt_lcg_rand(ctx->seed); 1163 /* About 48 dB below the "normal" folding level */ 1164 tmp = QCONST16(1.0f/256, NORM_SHIFT-4); 1165 tmp = (ctx->seed)&0x8000 ? tmp : -tmp; 1166 X[j] = lowband[j]+tmp; 1167 } 1168 cm = fill; 1169 } 1170 renormalise_vector(X, N, gain, ctx->arch); 1171 } 1172 } 1173 } 1174 } 1175 1176 return cm; 1177 } 1178 1179 #ifdef ENABLE_QEXT 1180 static unsigned cubic_quant_partition(struct band_ctx *ctx, celt_norm *X, int N, int b, int B, ec_ctx *ec, int LM, opus_val32 gain, int resynth, int encode) 1181 { 1182 celt_assert(LM>=0); 1183 ctx->remaining_bits = ctx->ec->storage*8*8 - ec_tell_frac(ctx->ec); 1184 b = IMIN(b, ctx->remaining_bits); 1185 /* As long as we have at least two bits of depth, split all the way to LM=0 (not -1 like PVQ). */ 1186 if (LM==0 || b<=2*N<<BITRES) { 1187 int res, ret; 1188 b = IMIN(b + ((N-1)<<BITRES)/2, ctx->remaining_bits); 1189 /* Resolution left after taking into account coding the cube face. */ 1190 res = (b-(1<<BITRES)-ctx->m->logN[ctx->i]-(LM<<BITRES)-1)/(N-1)>>BITRES; 1191 res = IMIN(14, IMAX(0, res)); 1192 if (encode) ret = cubic_quant(X, N, res, B, ec, gain, resynth); 1193 else ret = cubic_unquant(X, N, res, B, ec, gain); 1194 ctx->remaining_bits = ctx->ec->storage*8*8 - ec_tell_frac(ctx->ec); 1195 return ret; 1196 } else { 1197 celt_norm *Y; 1198 opus_int32 itheta_q30; 1199 opus_val32 g1, g2; 1200 opus_int32 theta_res; 1201 opus_int32 qtheta; 1202 int delta; 1203 int b1, b2; 1204 int cm; 1205 int N0; 1206 N0 = N; 1207 N >>= 1; 1208 Y = X+N; 1209 LM -= 1; 1210 B = (B+1)>>1; 1211 theta_res = IMIN(16, (b>>BITRES)/(N0-1) + 1); 1212 if (encode) { 1213 itheta_q30 = stereo_itheta(X, Y, 0, N, ctx->arch); 1214 qtheta = (itheta_q30+(1<<(29-theta_res)))>>(30-theta_res); 1215 ec_enc_uint(ec, qtheta, (1<<theta_res)+1); 1216 } else { 1217 qtheta = ec_dec_uint(ec, (1<<theta_res)+1); 1218 } 1219 itheta_q30 = qtheta<<(30-theta_res); 1220 b -= theta_res<<BITRES; 1221 delta = (N0-1) * 23 * ((itheta_q30>>16)-8192) >> (17-BITRES); 1222 1223 #ifdef FIXED_POINT 1224 g1 = celt_cos_norm32(itheta_q30); 1225 g2 = celt_cos_norm32((1<<30)-itheta_q30); 1226 #else 1227 g1 = celt_cos_norm2(itheta_q30*(1.f/(1<<30))); 1228 g2 = celt_cos_norm2(1.f-itheta_q30*(1.f/(1<<30))); 1229 #endif 1230 if (itheta_q30 == 0) { 1231 b1=b; 1232 b2=0; 1233 } else if (itheta_q30==1073741824) { 1234 b1=0; 1235 b2=b; 1236 } else { 1237 b1 = IMIN(b, IMAX(0, (b-delta)/2)); 1238 b2 = b-b1; 1239 } 1240 cm = cubic_quant_partition(ctx, X, N, b1, B, ec, LM, MULT32_32_Q31(gain, g1), resynth, encode); 1241 cm |= cubic_quant_partition(ctx, Y, N, b2, B, ec, LM, MULT32_32_Q31(gain, g2), resynth, encode); 1242 return cm; 1243 } 1244 } 1245 #endif 1246 1247 /* This function is responsible for encoding and decoding a band for the mono case. */ 1248 static unsigned quant_band(struct band_ctx *ctx, celt_norm *X, 1249 int N, int b, int B, celt_norm *lowband, 1250 int LM, celt_norm *lowband_out, 1251 opus_val32 gain, celt_norm *lowband_scratch, int fill 1252 ARG_QEXT(int ext_b)) 1253 { 1254 int N0=N; 1255 int N_B=N; 1256 int N_B0; 1257 int B0=B; 1258 int time_divide=0; 1259 int recombine=0; 1260 int longBlocks; 1261 unsigned cm=0; 1262 int k; 1263 int encode; 1264 int tf_change; 1265 1266 encode = ctx->encode; 1267 tf_change = ctx->tf_change; 1268 1269 longBlocks = B0==1; 1270 1271 N_B = celt_udiv(N_B, B); 1272 1273 /* Special case for one sample */ 1274 if (N==1) 1275 { 1276 return quant_band_n1(ctx, X, NULL, lowband_out); 1277 } 1278 1279 if (tf_change>0) 1280 recombine = tf_change; 1281 /* Band recombining to increase frequency resolution */ 1282 1283 if (lowband_scratch && lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1)) 1284 { 1285 OPUS_COPY(lowband_scratch, lowband, N); 1286 lowband = lowband_scratch; 1287 } 1288 1289 for (k=0;k<recombine;k++) 1290 { 1291 static const unsigned char bit_interleave_table[16]={ 1292 0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3 1293 }; 1294 if (encode) 1295 haar1(X, N>>k, 1<<k); 1296 if (lowband) 1297 haar1(lowband, N>>k, 1<<k); 1298 fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2; 1299 } 1300 B>>=recombine; 1301 N_B<<=recombine; 1302 1303 /* Increasing the time resolution */ 1304 while ((N_B&1) == 0 && tf_change<0) 1305 { 1306 if (encode) 1307 haar1(X, N_B, B); 1308 if (lowband) 1309 haar1(lowband, N_B, B); 1310 fill |= fill<<B; 1311 B <<= 1; 1312 N_B >>= 1; 1313 time_divide++; 1314 tf_change++; 1315 } 1316 B0=B; 1317 N_B0 = N_B; 1318 1319 /* Reorganize the samples in time order instead of frequency order */ 1320 if (B0>1) 1321 { 1322 if (encode) 1323 deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks); 1324 if (lowband) 1325 deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks); 1326 } 1327 1328 #ifdef ENABLE_QEXT 1329 if (ctx->extra_bands && b > (3*N<<BITRES)+(ctx->m->logN[ctx->i]+8+8*LM)) { 1330 cm = cubic_quant_partition(ctx, X, N, b, B, ctx->ec, LM, gain, ctx->resynth, encode); 1331 } else 1332 #endif 1333 { 1334 cm = quant_partition(ctx, X, N, b, B, lowband, LM, gain, fill ARG_QEXT(ext_b)); 1335 } 1336 1337 /* This code is used by the decoder and by the resynthesis-enabled encoder */ 1338 if (ctx->resynth) 1339 { 1340 /* Undo the sample reorganization going from time order to frequency order */ 1341 if (B0>1) 1342 interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks); 1343 1344 /* Undo time-freq changes that we did earlier */ 1345 N_B = N_B0; 1346 B = B0; 1347 for (k=0;k<time_divide;k++) 1348 { 1349 B >>= 1; 1350 N_B <<= 1; 1351 cm |= cm>>B; 1352 haar1(X, N_B, B); 1353 } 1354 1355 for (k=0;k<recombine;k++) 1356 { 1357 static const unsigned char bit_deinterleave_table[16]={ 1358 0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F, 1359 0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF 1360 }; 1361 cm = bit_deinterleave_table[cm]; 1362 haar1(X, N0>>k, 1<<k); 1363 } 1364 B<<=recombine; 1365 1366 /* Scale output for later folding */ 1367 if (lowband_out) 1368 { 1369 int j; 1370 opus_val16 n; 1371 n = celt_sqrt(SHL32(EXTEND32(N0),22)); 1372 for (j=0;j<N0;j++) 1373 lowband_out[j] = MULT16_32_Q15(n,X[j]); 1374 } 1375 cm &= (1<<B)-1; 1376 } 1377 return cm; 1378 } 1379 1380 #ifdef FIXED_POINT 1381 #define MIN_STEREO_ENERGY 2 1382 #else 1383 #define MIN_STEREO_ENERGY 1e-10f 1384 #endif 1385 1386 /* This function is responsible for encoding and decoding a band for the stereo case. */ 1387 static unsigned quant_band_stereo(struct band_ctx *ctx, celt_norm *X, celt_norm *Y, 1388 int N, int b, int B, celt_norm *lowband, 1389 int LM, celt_norm *lowband_out, 1390 celt_norm *lowband_scratch, int fill 1391 ARG_QEXT(int ext_b) ARG_QEXT(const int *cap)) 1392 { 1393 int imid=0, iside=0; 1394 int inv = 0; 1395 opus_val32 mid=0, side=0; 1396 unsigned cm=0; 1397 int mbits, sbits, delta; 1398 int itheta; 1399 int qalloc; 1400 struct split_ctx sctx; 1401 int orig_fill; 1402 int encode; 1403 ec_ctx *ec; 1404 1405 encode = ctx->encode; 1406 ec = ctx->ec; 1407 1408 /* Special case for one sample */ 1409 if (N==1) 1410 { 1411 return quant_band_n1(ctx, X, Y, lowband_out); 1412 } 1413 1414 orig_fill = fill; 1415 1416 if (encode) { 1417 if (ctx->bandE[ctx->i] < MIN_STEREO_ENERGY || ctx->bandE[ctx->m->nbEBands+ctx->i] < MIN_STEREO_ENERGY) { 1418 if (ctx->bandE[ctx->i] > ctx->bandE[ctx->m->nbEBands+ctx->i]) OPUS_COPY(Y, X, N); 1419 else OPUS_COPY(X, Y, N); 1420 } 1421 } 1422 compute_theta(ctx, &sctx, X, Y, N, &b, B, B, LM, 1, &fill ARG_QEXT(&ext_b)); 1423 inv = sctx.inv; 1424 imid = sctx.imid; 1425 iside = sctx.iside; 1426 delta = sctx.delta; 1427 itheta = sctx.itheta; 1428 qalloc = sctx.qalloc; 1429 #ifdef FIXED_POINT 1430 # ifdef ENABLE_QEXT 1431 (void)imid; 1432 (void)iside; 1433 mid = celt_cos_norm32(sctx.itheta_q30); 1434 side = celt_cos_norm32((1<<30)-sctx.itheta_q30); 1435 # else 1436 mid = SHL32(EXTEND32(imid), 16); 1437 side = SHL32(EXTEND32(iside), 16); 1438 # endif 1439 #else 1440 # ifdef ENABLE_QEXT 1441 (void)imid; 1442 (void)iside; 1443 mid = celt_cos_norm2(sctx.itheta_q30*(1.f/(1<<30))); 1444 side = celt_cos_norm2(1.f-sctx.itheta_q30*(1.f/(1<<30))); 1445 # else 1446 mid = (1.f/32768)*imid; 1447 side = (1.f/32768)*iside; 1448 # endif 1449 #endif 1450 1451 /* This is a special case for N=2 that only works for stereo and takes 1452 advantage of the fact that mid and side are orthogonal to encode 1453 the side with just one bit. */ 1454 if (N==2) 1455 { 1456 int c; 1457 int sign=0; 1458 celt_norm *x2, *y2; 1459 mbits = b; 1460 sbits = 0; 1461 /* Only need one bit for the side. */ 1462 if (itheta != 0 && itheta != 16384) 1463 sbits = 1<<BITRES; 1464 mbits -= sbits; 1465 c = itheta > 8192; 1466 ctx->remaining_bits -= qalloc+sbits; 1467 1468 x2 = c ? Y : X; 1469 y2 = c ? X : Y; 1470 if (sbits) 1471 { 1472 if (encode) 1473 { 1474 /* Here we only need to encode a sign for the side. */ 1475 /* FIXME: Need to increase fixed-point precision? */ 1476 sign = MULT32_32_Q31(x2[0],y2[1]) - MULT32_32_Q31(x2[1],y2[0]) < 0; 1477 ec_enc_bits(ec, sign, 1); 1478 } else { 1479 sign = ec_dec_bits(ec, 1); 1480 } 1481 } 1482 sign = 1-2*sign; 1483 /* We use orig_fill here because we want to fold the side, but if 1484 itheta==16384, we'll have cleared the low bits of fill. */ 1485 cm = quant_band(ctx, x2, N, mbits, B, lowband, LM, lowband_out, Q31ONE, 1486 lowband_scratch, orig_fill ARG_QEXT(ext_b)); 1487 /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse), 1488 and there's no need to worry about mixing with the other channel. */ 1489 y2[0] = -sign*x2[1]; 1490 y2[1] = sign*x2[0]; 1491 if (ctx->resynth) 1492 { 1493 celt_norm tmp; 1494 X[0] = MULT32_32_Q31(mid, X[0]); 1495 X[1] = MULT32_32_Q31(mid, X[1]); 1496 Y[0] = MULT32_32_Q31(side, Y[0]); 1497 Y[1] = MULT32_32_Q31(side, Y[1]); 1498 tmp = X[0]; 1499 X[0] = SUB32(tmp,Y[0]); 1500 Y[0] = ADD32(tmp,Y[0]); 1501 tmp = X[1]; 1502 X[1] = SUB32(tmp,Y[1]); 1503 Y[1] = ADD32(tmp,Y[1]); 1504 } 1505 } else { 1506 /* "Normal" split code */ 1507 opus_int32 rebalance; 1508 1509 mbits = IMAX(0, IMIN(b, (b-delta)/2)); 1510 sbits = b-mbits; 1511 ctx->remaining_bits -= qalloc; 1512 1513 rebalance = ctx->remaining_bits; 1514 if (mbits >= sbits) 1515 { 1516 #ifdef ENABLE_QEXT 1517 int qext_extra = 0; 1518 /* Reallocate any mid bits that cannot be used to extra mid bits. */ 1519 if (cap != NULL && ext_b != 0) qext_extra = IMAX(0, IMIN(ext_b/2, mbits - cap[ctx->i]/2)); 1520 #endif 1521 /* In stereo mode, we do not apply a scaling to the mid because we need the normalized 1522 mid for folding later. */ 1523 cm = quant_band(ctx, X, N, mbits, B, lowband, LM, lowband_out, Q31ONE, 1524 lowband_scratch, fill ARG_QEXT(ext_b/2+qext_extra)); 1525 rebalance = mbits - (rebalance-ctx->remaining_bits); 1526 if (rebalance > 3<<BITRES && itheta!=0) 1527 sbits += rebalance - (3<<BITRES); 1528 #ifdef ENABLE_QEXT 1529 /* Guard against overflowing the EC with the angle if the cubic quant used too many bits for the mid. */ 1530 if (ctx->extra_bands) sbits = IMIN(sbits, ctx->remaining_bits); 1531 #endif 1532 /* For a stereo split, the high bits of fill are always zero, so no 1533 folding will be done to the side. */ 1534 cm |= quant_band(ctx, Y, N, sbits, B, NULL, LM, NULL, side, NULL, fill>>B ARG_QEXT(ext_b/2-qext_extra)); 1535 } else { 1536 #ifdef ENABLE_QEXT 1537 int qext_extra = 0; 1538 /* Reallocate any side bits that cannot be used to extra side bits. */ 1539 if (cap != NULL && ext_b != 0) qext_extra = IMAX(0, IMIN(ext_b/2, sbits - cap[ctx->i]/2)); 1540 #endif 1541 /* For a stereo split, the high bits of fill are always zero, so no 1542 folding will be done to the side. */ 1543 cm = quant_band(ctx, Y, N, sbits, B, NULL, LM, NULL, side, NULL, fill>>B ARG_QEXT(ext_b/2+qext_extra)); 1544 rebalance = sbits - (rebalance-ctx->remaining_bits); 1545 if (rebalance > 3<<BITRES && itheta!=16384) 1546 mbits += rebalance - (3<<BITRES); 1547 #ifdef ENABLE_QEXT 1548 /* Guard against overflowing the EC with the angle if the cubic quant used too many bits for the side. */ 1549 if (ctx->extra_bands) mbits = IMIN(mbits, ctx->remaining_bits); 1550 #endif 1551 /* In stereo mode, we do not apply a scaling to the mid because we need the normalized 1552 mid for folding later. */ 1553 cm |= quant_band(ctx, X, N, mbits, B, lowband, LM, lowband_out, Q31ONE, 1554 lowband_scratch, fill ARG_QEXT(ext_b/2-qext_extra)); 1555 } 1556 } 1557 1558 1559 /* This code is used by the decoder and by the resynthesis-enabled encoder */ 1560 if (ctx->resynth) 1561 { 1562 if (N!=2) 1563 stereo_merge(X, Y, mid, N, ctx->arch); 1564 if (inv) 1565 { 1566 int j; 1567 for (j=0;j<N;j++) 1568 Y[j] = -Y[j]; 1569 } 1570 } 1571 return cm; 1572 } 1573 1574 #ifndef DISABLE_UPDATE_DRAFT 1575 static void special_hybrid_folding(const CELTMode *m, celt_norm *norm, celt_norm *norm2, int start, int M, int dual_stereo) 1576 { 1577 int n1, n2; 1578 const opus_int16 * OPUS_RESTRICT eBands = m->eBands; 1579 n1 = M*(eBands[start+1]-eBands[start]); 1580 n2 = M*(eBands[start+2]-eBands[start+1]); 1581 /* Duplicate enough of the first band folding data to be able to fold the second band. 1582 Copies no data for CELT-only mode. */ 1583 OPUS_COPY(&norm[n1], &norm[2*n1 - n2], n2-n1); 1584 if (dual_stereo) 1585 OPUS_COPY(&norm2[n1], &norm2[2*n1 - n2], n2-n1); 1586 } 1587 #endif 1588 1589 void quant_all_bands(int encode, const CELTMode *m, int start, int end, 1590 celt_norm *X_, celt_norm *Y_, unsigned char *collapse_masks, 1591 const celt_ener *bandE, int *pulses, int shortBlocks, int spread, 1592 int dual_stereo, int intensity, int *tf_res, opus_int32 total_bits, 1593 opus_int32 balance, ec_ctx *ec, int LM, int codedBands, 1594 opus_uint32 *seed, int complexity, int arch, int disable_inv 1595 ARG_QEXT(ec_ctx *ext_ec) ARG_QEXT(int *extra_pulses) 1596 ARG_QEXT(opus_int32 ext_total_bits) ARG_QEXT(const int *cap)) 1597 { 1598 int i; 1599 opus_int32 remaining_bits; 1600 const opus_int16 * OPUS_RESTRICT eBands = m->eBands; 1601 celt_norm * OPUS_RESTRICT norm, * OPUS_RESTRICT norm2; 1602 VARDECL(celt_norm, _norm); 1603 VARDECL(celt_norm, _lowband_scratch); 1604 VARDECL(celt_norm, X_save); 1605 VARDECL(celt_norm, Y_save); 1606 VARDECL(celt_norm, X_save2); 1607 VARDECL(celt_norm, Y_save2); 1608 VARDECL(celt_norm, norm_save2); 1609 VARDECL(unsigned char, bytes_save); 1610 int resynth_alloc; 1611 celt_norm *lowband_scratch; 1612 int B; 1613 int M; 1614 int lowband_offset; 1615 int update_lowband = 1; 1616 int C = Y_ != NULL ? 2 : 1; 1617 int norm_offset; 1618 int theta_rdo = encode && Y_!=NULL && !dual_stereo && complexity>=8; 1619 #ifdef RESYNTH 1620 int resynth = 1; 1621 #else 1622 int resynth = !encode || theta_rdo; 1623 #endif 1624 struct band_ctx ctx; 1625 #ifdef ENABLE_QEXT 1626 int ext_b; 1627 opus_int32 ext_balance=0; 1628 opus_int32 ext_tell=0; 1629 VARDECL(unsigned char, ext_bytes_save); 1630 #endif 1631 SAVE_STACK; 1632 1633 M = 1<<LM; 1634 B = shortBlocks ? M : 1; 1635 norm_offset = M*eBands[start]; 1636 /* No need to allocate norm for the last band because we don't need an 1637 output in that band. */ 1638 ALLOC(_norm, C*(M*eBands[m->nbEBands-1]-norm_offset), celt_norm); 1639 norm = _norm; 1640 norm2 = norm + M*eBands[m->nbEBands-1]-norm_offset; 1641 1642 /* For decoding, we can use the last band as scratch space because we don't need that 1643 scratch space for the last band and we don't care about the data there until we're 1644 decoding the last band. */ 1645 if (encode && resynth) 1646 resynth_alloc = M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]); 1647 else 1648 resynth_alloc = ALLOC_NONE; 1649 ALLOC(_lowband_scratch, resynth_alloc, celt_norm); 1650 if (encode && resynth) 1651 lowband_scratch = _lowband_scratch; 1652 else 1653 lowband_scratch = X_+M*eBands[m->effEBands-1]; 1654 ALLOC(X_save, resynth_alloc, celt_norm); 1655 ALLOC(Y_save, resynth_alloc, celt_norm); 1656 ALLOC(X_save2, resynth_alloc, celt_norm); 1657 ALLOC(Y_save2, resynth_alloc, celt_norm); 1658 ALLOC(norm_save2, resynth_alloc, celt_norm); 1659 1660 lowband_offset = 0; 1661 ctx.bandE = bandE; 1662 ctx.ec = ec; 1663 ctx.encode = encode; 1664 ctx.intensity = intensity; 1665 ctx.m = m; 1666 ctx.seed = *seed; 1667 ctx.spread = spread; 1668 ctx.arch = arch; 1669 ctx.disable_inv = disable_inv; 1670 ctx.resynth = resynth; 1671 ctx.theta_round = 0; 1672 #ifdef ENABLE_QEXT 1673 ctx.ext_ec = ext_ec; 1674 ctx.ext_total_bits = ext_total_bits; 1675 ctx.extra_bands = end == NB_QEXT_BANDS || end == 2; 1676 if (ctx.extra_bands) theta_rdo = 0; 1677 ALLOC(ext_bytes_save, theta_rdo ? QEXT_PACKET_SIZE_CAP : ALLOC_NONE, unsigned char); 1678 #endif 1679 ALLOC(bytes_save, theta_rdo ? 1275 : ALLOC_NONE, unsigned char); 1680 1681 /* Avoid injecting noise in the first band on transients. */ 1682 ctx.avoid_split_noise = B > 1; 1683 for (i=start;i<end;i++) 1684 { 1685 opus_int32 tell; 1686 int b; 1687 int N; 1688 opus_int32 curr_balance; 1689 int effective_lowband=-1; 1690 celt_norm * OPUS_RESTRICT X, * OPUS_RESTRICT Y; 1691 int tf_change=0; 1692 unsigned x_cm; 1693 unsigned y_cm; 1694 int last; 1695 1696 ctx.i = i; 1697 last = (i==end-1); 1698 1699 X = X_+M*eBands[i]; 1700 if (Y_!=NULL) 1701 Y = Y_+M*eBands[i]; 1702 else 1703 Y = NULL; 1704 N = M*eBands[i+1]-M*eBands[i]; 1705 celt_assert(N > 0); 1706 tell = ec_tell_frac(ec); 1707 1708 /* Compute how many bits we want to allocate to this band */ 1709 if (i != start) 1710 balance -= tell; 1711 remaining_bits = total_bits-tell-1; 1712 ctx.remaining_bits = remaining_bits; 1713 #ifdef ENABLE_QEXT 1714 if (i != start) { 1715 ext_balance += extra_pulses[i-1] + ext_tell; 1716 } 1717 ext_tell = ec_tell_frac(ext_ec); 1718 ctx.extra_bits = extra_pulses[i]; 1719 if (i != start) 1720 ext_balance -= ext_tell; 1721 if (i <= codedBands-1) 1722 { 1723 opus_int32 ext_curr_balance = celt_sudiv(ext_balance, IMIN(3, codedBands-i)); 1724 ext_b = IMAX(0, IMIN(16383, IMIN(ext_total_bits-ext_tell,extra_pulses[i]+ext_curr_balance))); 1725 } else { 1726 ext_b = 0; 1727 } 1728 #endif 1729 if (i <= codedBands-1) 1730 { 1731 curr_balance = celt_sudiv(balance, IMIN(3, codedBands-i)); 1732 b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance))); 1733 } else { 1734 b = 0; 1735 } 1736 1737 #ifndef DISABLE_UPDATE_DRAFT 1738 if (resynth && (M*eBands[i]-N >= M*eBands[start] || i==start+1) && (update_lowband || lowband_offset==0)) 1739 lowband_offset = i; 1740 if (i == start+1) 1741 special_hybrid_folding(m, norm, norm2, start, M, dual_stereo); 1742 #else 1743 if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0)) 1744 lowband_offset = i; 1745 #endif 1746 1747 tf_change = tf_res[i]; 1748 ctx.tf_change = tf_change; 1749 if (i>=m->effEBands) 1750 { 1751 X=norm; 1752 if (Y_!=NULL) 1753 Y = norm; 1754 lowband_scratch = NULL; 1755 } 1756 if (last && !theta_rdo) 1757 lowband_scratch = NULL; 1758 1759 /* Get a conservative estimate of the collapse_mask's for the bands we're 1760 going to be folding from. */ 1761 if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0)) 1762 { 1763 int fold_start; 1764 int fold_end; 1765 int fold_i; 1766 /* This ensures we never repeat spectral content within one band */ 1767 effective_lowband = IMAX(0, M*eBands[lowband_offset]-norm_offset-N); 1768 fold_start = lowband_offset; 1769 while(M*eBands[--fold_start] > effective_lowband+norm_offset); 1770 fold_end = lowband_offset-1; 1771 #ifndef DISABLE_UPDATE_DRAFT 1772 while(++fold_end < i && M*eBands[fold_end] < effective_lowband+norm_offset+N); 1773 #else 1774 while(M*eBands[++fold_end] < effective_lowband+norm_offset+N); 1775 #endif 1776 x_cm = y_cm = 0; 1777 fold_i = fold_start; do { 1778 x_cm |= collapse_masks[fold_i*C+0]; 1779 y_cm |= collapse_masks[fold_i*C+C-1]; 1780 } while (++fold_i<fold_end); 1781 } 1782 /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost 1783 always) be non-zero. */ 1784 else 1785 x_cm = y_cm = (1<<B)-1; 1786 1787 if (dual_stereo && i==intensity) 1788 { 1789 int j; 1790 1791 /* Switch off dual stereo to do intensity. */ 1792 dual_stereo = 0; 1793 if (resynth) 1794 for (j=0;j<M*eBands[i]-norm_offset;j++) 1795 norm[j] = HALF32(norm[j]+norm2[j]); 1796 } 1797 if (dual_stereo) 1798 { 1799 x_cm = quant_band(&ctx, X, N, b/2, B, 1800 effective_lowband != -1 ? norm+effective_lowband : NULL, LM, 1801 last?NULL:norm+M*eBands[i]-norm_offset, Q31ONE, lowband_scratch, x_cm ARG_QEXT(ext_b/2)); 1802 y_cm = quant_band(&ctx, Y, N, b/2, B, 1803 effective_lowband != -1 ? norm2+effective_lowband : NULL, LM, 1804 last?NULL:norm2+M*eBands[i]-norm_offset, Q31ONE, lowband_scratch, y_cm ARG_QEXT(ext_b/2)); 1805 } else { 1806 if (Y!=NULL) 1807 { 1808 if (theta_rdo && i < intensity) 1809 { 1810 ec_ctx ec_save, ec_save2; 1811 struct band_ctx ctx_save, ctx_save2; 1812 opus_val32 dist0, dist1; 1813 unsigned cm, cm2; 1814 int nstart_bytes, nend_bytes, save_bytes; 1815 unsigned char *bytes_buf; 1816 #ifdef ENABLE_QEXT 1817 ec_ctx ext_ec_save, ext_ec_save2; 1818 unsigned char *ext_bytes_buf; 1819 int ext_nstart_bytes, ext_nend_bytes, ext_save_bytes; 1820 #endif 1821 opus_val16 w[2]; 1822 compute_channel_weights(bandE[i], bandE[i+m->nbEBands], w); 1823 /* Make a copy. */ 1824 cm = x_cm|y_cm; 1825 ec_save = *ec; 1826 #ifdef ENABLE_QEXT 1827 ext_ec_save = *ext_ec; 1828 #endif 1829 ctx_save = ctx; 1830 OPUS_COPY(X_save, X, N); 1831 OPUS_COPY(Y_save, Y, N); 1832 /* Encode and round down. */ 1833 ctx.theta_round = -1; 1834 x_cm = quant_band_stereo(&ctx, X, Y, N, b, B, 1835 effective_lowband != -1 ? norm+effective_lowband : NULL, LM, 1836 last?NULL:norm+M*eBands[i]-norm_offset, lowband_scratch, cm ARG_QEXT(ext_b) ARG_QEXT(cap)); 1837 dist0 = MULT16_32_Q15(w[0], celt_inner_prod_norm_shift(X_save, X, N, arch)) + MULT16_32_Q15(w[1], celt_inner_prod_norm_shift(Y_save, Y, N, arch)); 1838 1839 /* Save first result. */ 1840 cm2 = x_cm; 1841 ec_save2 = *ec; 1842 #ifdef ENABLE_QEXT 1843 ext_ec_save2 = *ext_ec; 1844 #endif 1845 ctx_save2 = ctx; 1846 OPUS_COPY(X_save2, X, N); 1847 OPUS_COPY(Y_save2, Y, N); 1848 if (!last) 1849 OPUS_COPY(norm_save2, norm+M*eBands[i]-norm_offset, N); 1850 nstart_bytes = ec_save.offs; 1851 nend_bytes = ec_save.storage; 1852 bytes_buf = ec_save.buf+nstart_bytes; 1853 save_bytes = nend_bytes-nstart_bytes; 1854 OPUS_COPY(bytes_save, bytes_buf, save_bytes); 1855 #ifdef ENABLE_QEXT 1856 ext_nstart_bytes = ext_ec_save.offs; 1857 ext_nend_bytes = ext_ec_save.storage; 1858 ext_bytes_buf = ext_ec_save.buf!=NULL ? ext_ec_save.buf+ext_nstart_bytes : NULL; 1859 ext_save_bytes = ext_nend_bytes-ext_nstart_bytes; 1860 if (ext_save_bytes) OPUS_COPY(ext_bytes_save, ext_bytes_buf, ext_save_bytes); 1861 #endif 1862 /* Restore */ 1863 *ec = ec_save; 1864 #ifdef ENABLE_QEXT 1865 *ext_ec = ext_ec_save; 1866 #endif 1867 ctx = ctx_save; 1868 OPUS_COPY(X, X_save, N); 1869 OPUS_COPY(Y, Y_save, N); 1870 #ifndef DISABLE_UPDATE_DRAFT 1871 if (i == start+1) 1872 special_hybrid_folding(m, norm, norm2, start, M, dual_stereo); 1873 #endif 1874 /* Encode and round up. */ 1875 ctx.theta_round = 1; 1876 x_cm = quant_band_stereo(&ctx, X, Y, N, b, B, 1877 effective_lowband != -1 ? norm+effective_lowband : NULL, LM, 1878 last?NULL:norm+M*eBands[i]-norm_offset, lowband_scratch, cm ARG_QEXT(ext_b) ARG_QEXT(cap)); 1879 dist1 = MULT16_32_Q15(w[0], celt_inner_prod_norm_shift(X_save, X, N, arch)) + MULT16_32_Q15(w[1], celt_inner_prod_norm_shift(Y_save, Y, N, arch)); 1880 if (dist0 >= dist1) { 1881 x_cm = cm2; 1882 *ec = ec_save2; 1883 #ifdef ENABLE_QEXT 1884 *ext_ec = ext_ec_save2; 1885 #endif 1886 ctx = ctx_save2; 1887 OPUS_COPY(X, X_save2, N); 1888 OPUS_COPY(Y, Y_save2, N); 1889 if (!last) 1890 OPUS_COPY(norm+M*eBands[i]-norm_offset, norm_save2, N); 1891 OPUS_COPY(bytes_buf, bytes_save, save_bytes); 1892 #ifdef ENABLE_QEXT 1893 if (ext_save_bytes) OPUS_COPY(ext_bytes_buf, ext_bytes_save, ext_save_bytes); 1894 #endif 1895 } 1896 } else { 1897 ctx.theta_round = 0; 1898 x_cm = quant_band_stereo(&ctx, X, Y, N, b, B, 1899 effective_lowband != -1 ? norm+effective_lowband : NULL, LM, 1900 last?NULL:norm+M*eBands[i]-norm_offset, lowband_scratch, x_cm|y_cm ARG_QEXT(ext_b) ARG_QEXT(cap)); 1901 } 1902 } else { 1903 x_cm = quant_band(&ctx, X, N, b, B, 1904 effective_lowband != -1 ? norm+effective_lowband : NULL, LM, 1905 last?NULL:norm+M*eBands[i]-norm_offset, Q31ONE, lowband_scratch, x_cm|y_cm ARG_QEXT(ext_b)); 1906 } 1907 y_cm = x_cm; 1908 } 1909 collapse_masks[i*C+0] = (unsigned char)x_cm; 1910 collapse_masks[i*C+C-1] = (unsigned char)y_cm; 1911 balance += pulses[i] + tell; 1912 1913 /* Update the folding position only as long as we have 1 bit/sample depth. */ 1914 update_lowband = b>(N<<BITRES); 1915 /* We only need to avoid noise on a split for the first band. After that, we 1916 have folding. */ 1917 ctx.avoid_split_noise = 0; 1918 } 1919 *seed = ctx.seed; 1920 1921 RESTORE_STACK; 1922 }