1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libopus/celt/bands.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,1518 @@
1.4 +/* Copyright (c) 2007-2008 CSIRO
1.5 + Copyright (c) 2007-2009 Xiph.Org Foundation
1.6 + Copyright (c) 2008-2009 Gregory Maxwell
1.7 + Written by Jean-Marc Valin and Gregory Maxwell */
1.8 +/*
1.9 + Redistribution and use in source and binary forms, with or without
1.10 + modification, are permitted provided that the following conditions
1.11 + are met:
1.12 +
1.13 + - Redistributions of source code must retain the above copyright
1.14 + notice, this list of conditions and the following disclaimer.
1.15 +
1.16 + - Redistributions in binary form must reproduce the above copyright
1.17 + notice, this list of conditions and the following disclaimer in the
1.18 + documentation and/or other materials provided with the distribution.
1.19 +
1.20 + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
1.21 + ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
1.22 + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
1.23 + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
1.24 + OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
1.25 + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
1.26 + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
1.27 + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
1.28 + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
1.29 + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
1.30 + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
1.31 +*/
1.32 +
1.33 +#ifdef HAVE_CONFIG_H
1.34 +#include "config.h"
1.35 +#endif
1.36 +
1.37 +#include <math.h>
1.38 +#include "bands.h"
1.39 +#include "modes.h"
1.40 +#include "vq.h"
1.41 +#include "cwrs.h"
1.42 +#include "stack_alloc.h"
1.43 +#include "os_support.h"
1.44 +#include "mathops.h"
1.45 +#include "rate.h"
1.46 +#include "quant_bands.h"
1.47 +#include "pitch.h"
1.48 +
1.49 +int hysteresis_decision(opus_val16 val, const opus_val16 *thresholds, const opus_val16 *hysteresis, int N, int prev)
1.50 +{
1.51 + int i;
1.52 + for (i=0;i<N;i++)
1.53 + {
1.54 + if (val < thresholds[i])
1.55 + break;
1.56 + }
1.57 + if (i>prev && val < thresholds[prev]+hysteresis[prev])
1.58 + i=prev;
1.59 + if (i<prev && val > thresholds[prev-1]-hysteresis[prev-1])
1.60 + i=prev;
1.61 + return i;
1.62 +}
1.63 +
1.64 +opus_uint32 celt_lcg_rand(opus_uint32 seed)
1.65 +{
1.66 + return 1664525 * seed + 1013904223;
1.67 +}
1.68 +
1.69 +/* This is a cos() approximation designed to be bit-exact on any platform. Bit exactness
1.70 + with this approximation is important because it has an impact on the bit allocation */
1.71 +static opus_int16 bitexact_cos(opus_int16 x)
1.72 +{
1.73 + opus_int32 tmp;
1.74 + opus_int16 x2;
1.75 + tmp = (4096+((opus_int32)(x)*(x)))>>13;
1.76 + celt_assert(tmp<=32767);
1.77 + x2 = tmp;
1.78 + x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));
1.79 + celt_assert(x2<=32766);
1.80 + return 1+x2;
1.81 +}
1.82 +
1.83 +static int bitexact_log2tan(int isin,int icos)
1.84 +{
1.85 + int lc;
1.86 + int ls;
1.87 + lc=EC_ILOG(icos);
1.88 + ls=EC_ILOG(isin);
1.89 + icos<<=15-lc;
1.90 + isin<<=15-ls;
1.91 + return (ls-lc)*(1<<11)
1.92 + +FRAC_MUL16(isin, FRAC_MUL16(isin, -2597) + 7932)
1.93 + -FRAC_MUL16(icos, FRAC_MUL16(icos, -2597) + 7932);
1.94 +}
1.95 +
1.96 +#ifdef FIXED_POINT
1.97 +/* Compute the amplitude (sqrt energy) in each of the bands */
1.98 +void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int M)
1.99 +{
1.100 + int i, c, N;
1.101 + const opus_int16 *eBands = m->eBands;
1.102 + N = M*m->shortMdctSize;
1.103 + c=0; do {
1.104 + for (i=0;i<end;i++)
1.105 + {
1.106 + int j;
1.107 + opus_val32 maxval=0;
1.108 + opus_val32 sum = 0;
1.109 +
1.110 + j=M*eBands[i]; do {
1.111 + maxval = MAX32(maxval, X[j+c*N]);
1.112 + maxval = MAX32(maxval, -X[j+c*N]);
1.113 + } while (++j<M*eBands[i+1]);
1.114 +
1.115 + if (maxval > 0)
1.116 + {
1.117 + int shift = celt_ilog2(maxval)-10;
1.118 + j=M*eBands[i]; do {
1.119 + sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
1.120 + EXTRACT16(VSHR32(X[j+c*N],shift)));
1.121 + } while (++j<M*eBands[i+1]);
1.122 + /* We're adding one here to ensure the normalized band isn't larger than unity norm */
1.123 + bandE[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
1.124 + } else {
1.125 + bandE[i+c*m->nbEBands] = EPSILON;
1.126 + }
1.127 + /*printf ("%f ", bandE[i+c*m->nbEBands]);*/
1.128 + }
1.129 + } while (++c<C);
1.130 + /*printf ("\n");*/
1.131 +}
1.132 +
1.133 +/* Normalise each band such that the energy is one. */
1.134 +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)
1.135 +{
1.136 + int i, c, N;
1.137 + const opus_int16 *eBands = m->eBands;
1.138 + N = M*m->shortMdctSize;
1.139 + c=0; do {
1.140 + i=0; do {
1.141 + opus_val16 g;
1.142 + int j,shift;
1.143 + opus_val16 E;
1.144 + shift = celt_zlog2(bandE[i+c*m->nbEBands])-13;
1.145 + E = VSHR32(bandE[i+c*m->nbEBands], shift);
1.146 + g = EXTRACT16(celt_rcp(SHL32(E,3)));
1.147 + j=M*eBands[i]; do {
1.148 + X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
1.149 + } while (++j<M*eBands[i+1]);
1.150 + } while (++i<end);
1.151 + } while (++c<C);
1.152 +}
1.153 +
1.154 +#else /* FIXED_POINT */
1.155 +/* Compute the amplitude (sqrt energy) in each of the bands */
1.156 +void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int M)
1.157 +{
1.158 + int i, c, N;
1.159 + const opus_int16 *eBands = m->eBands;
1.160 + N = M*m->shortMdctSize;
1.161 + c=0; do {
1.162 + for (i=0;i<end;i++)
1.163 + {
1.164 + int j;
1.165 + opus_val32 sum = 1e-27f;
1.166 + for (j=M*eBands[i];j<M*eBands[i+1];j++)
1.167 + sum += X[j+c*N]*X[j+c*N];
1.168 + bandE[i+c*m->nbEBands] = celt_sqrt(sum);
1.169 + /*printf ("%f ", bandE[i+c*m->nbEBands]);*/
1.170 + }
1.171 + } while (++c<C);
1.172 + /*printf ("\n");*/
1.173 +}
1.174 +
1.175 +/* Normalise each band such that the energy is one. */
1.176 +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)
1.177 +{
1.178 + int i, c, N;
1.179 + const opus_int16 *eBands = m->eBands;
1.180 + N = M*m->shortMdctSize;
1.181 + c=0; do {
1.182 + for (i=0;i<end;i++)
1.183 + {
1.184 + int j;
1.185 + opus_val16 g = 1.f/(1e-27f+bandE[i+c*m->nbEBands]);
1.186 + for (j=M*eBands[i];j<M*eBands[i+1];j++)
1.187 + X[j+c*N] = freq[j+c*N]*g;
1.188 + }
1.189 + } while (++c<C);
1.190 +}
1.191 +
1.192 +#endif /* FIXED_POINT */
1.193 +
1.194 +/* De-normalise the energy to produce the synthesis from the unit-energy bands */
1.195 +void denormalise_bands(const CELTMode *m, const celt_norm * OPUS_RESTRICT X,
1.196 + celt_sig * OPUS_RESTRICT freq, const opus_val16 *bandLogE, int start, int end, int C, int M)
1.197 +{
1.198 + int i, c, N;
1.199 + const opus_int16 *eBands = m->eBands;
1.200 + N = M*m->shortMdctSize;
1.201 + celt_assert2(C<=2, "denormalise_bands() not implemented for >2 channels");
1.202 + c=0; do {
1.203 + celt_sig * OPUS_RESTRICT f;
1.204 + const celt_norm * OPUS_RESTRICT x;
1.205 + f = freq+c*N;
1.206 + x = X+c*N+M*eBands[start];
1.207 + for (i=0;i<M*eBands[start];i++)
1.208 + *f++ = 0;
1.209 + for (i=start;i<end;i++)
1.210 + {
1.211 + int j, band_end;
1.212 + opus_val16 g;
1.213 + opus_val16 lg;
1.214 +#ifdef FIXED_POINT
1.215 + int shift;
1.216 +#endif
1.217 + j=M*eBands[i];
1.218 + band_end = M*eBands[i+1];
1.219 + lg = ADD16(bandLogE[i+c*m->nbEBands], SHL16((opus_val16)eMeans[i],6));
1.220 +#ifndef FIXED_POINT
1.221 + g = celt_exp2(lg);
1.222 +#else
1.223 + /* Handle the integer part of the log energy */
1.224 + shift = 16-(lg>>DB_SHIFT);
1.225 + if (shift>31)
1.226 + {
1.227 + shift=0;
1.228 + g=0;
1.229 + } else {
1.230 + /* Handle the fractional part. */
1.231 + g = celt_exp2_frac(lg&((1<<DB_SHIFT)-1));
1.232 + }
1.233 + /* Handle extreme gains with negative shift. */
1.234 + if (shift<0)
1.235 + {
1.236 + /* For shift < -2 we'd be likely to overflow, so we're capping
1.237 + the gain here. This shouldn't happen unless the bitstream is
1.238 + already corrupted. */
1.239 + if (shift < -2)
1.240 + {
1.241 + g = 32767;
1.242 + shift = -2;
1.243 + }
1.244 + do {
1.245 + *f++ = SHL32(MULT16_16(*x++, g), -shift);
1.246 + } while (++j<band_end);
1.247 + } else
1.248 +#endif
1.249 + /* Be careful of the fixed-point "else" just above when changing this code */
1.250 + do {
1.251 + *f++ = SHR32(MULT16_16(*x++, g), shift);
1.252 + } while (++j<band_end);
1.253 + }
1.254 + celt_assert(start <= end);
1.255 + for (i=M*eBands[end];i<N;i++)
1.256 + *f++ = 0;
1.257 + } while (++c<C);
1.258 +}
1.259 +
1.260 +/* This prevents energy collapse for transients with multiple short MDCTs */
1.261 +void anti_collapse(const CELTMode *m, celt_norm *X_, unsigned char *collapse_masks, int LM, int C, int size,
1.262 + int start, int end, opus_val16 *logE, opus_val16 *prev1logE,
1.263 + opus_val16 *prev2logE, int *pulses, opus_uint32 seed)
1.264 +{
1.265 + int c, i, j, k;
1.266 + for (i=start;i<end;i++)
1.267 + {
1.268 + int N0;
1.269 + opus_val16 thresh, sqrt_1;
1.270 + int depth;
1.271 +#ifdef FIXED_POINT
1.272 + int shift;
1.273 + opus_val32 thresh32;
1.274 +#endif
1.275 +
1.276 + N0 = m->eBands[i+1]-m->eBands[i];
1.277 + /* depth in 1/8 bits */
1.278 + depth = (1+pulses[i])/((m->eBands[i+1]-m->eBands[i])<<LM);
1.279 +
1.280 +#ifdef FIXED_POINT
1.281 + thresh32 = SHR32(celt_exp2(-SHL16(depth, 10-BITRES)),1);
1.282 + thresh = MULT16_32_Q15(QCONST16(0.5f, 15), MIN32(32767,thresh32));
1.283 + {
1.284 + opus_val32 t;
1.285 + t = N0<<LM;
1.286 + shift = celt_ilog2(t)>>1;
1.287 + t = SHL32(t, (7-shift)<<1);
1.288 + sqrt_1 = celt_rsqrt_norm(t);
1.289 + }
1.290 +#else
1.291 + thresh = .5f*celt_exp2(-.125f*depth);
1.292 + sqrt_1 = celt_rsqrt(N0<<LM);
1.293 +#endif
1.294 +
1.295 + c=0; do
1.296 + {
1.297 + celt_norm *X;
1.298 + opus_val16 prev1;
1.299 + opus_val16 prev2;
1.300 + opus_val32 Ediff;
1.301 + opus_val16 r;
1.302 + int renormalize=0;
1.303 + prev1 = prev1logE[c*m->nbEBands+i];
1.304 + prev2 = prev2logE[c*m->nbEBands+i];
1.305 + if (C==1)
1.306 + {
1.307 + prev1 = MAX16(prev1,prev1logE[m->nbEBands+i]);
1.308 + prev2 = MAX16(prev2,prev2logE[m->nbEBands+i]);
1.309 + }
1.310 + Ediff = EXTEND32(logE[c*m->nbEBands+i])-EXTEND32(MIN16(prev1,prev2));
1.311 + Ediff = MAX32(0, Ediff);
1.312 +
1.313 +#ifdef FIXED_POINT
1.314 + if (Ediff < 16384)
1.315 + {
1.316 + opus_val32 r32 = SHR32(celt_exp2(-EXTRACT16(Ediff)),1);
1.317 + r = 2*MIN16(16383,r32);
1.318 + } else {
1.319 + r = 0;
1.320 + }
1.321 + if (LM==3)
1.322 + r = MULT16_16_Q14(23170, MIN32(23169, r));
1.323 + r = SHR16(MIN16(thresh, r),1);
1.324 + r = SHR32(MULT16_16_Q15(sqrt_1, r),shift);
1.325 +#else
1.326 + /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because
1.327 + short blocks don't have the same energy as long */
1.328 + r = 2.f*celt_exp2(-Ediff);
1.329 + if (LM==3)
1.330 + r *= 1.41421356f;
1.331 + r = MIN16(thresh, r);
1.332 + r = r*sqrt_1;
1.333 +#endif
1.334 + X = X_+c*size+(m->eBands[i]<<LM);
1.335 + for (k=0;k<1<<LM;k++)
1.336 + {
1.337 + /* Detect collapse */
1.338 + if (!(collapse_masks[i*C+c]&1<<k))
1.339 + {
1.340 + /* Fill with noise */
1.341 + for (j=0;j<N0;j++)
1.342 + {
1.343 + seed = celt_lcg_rand(seed);
1.344 + X[(j<<LM)+k] = (seed&0x8000 ? r : -r);
1.345 + }
1.346 + renormalize = 1;
1.347 + }
1.348 + }
1.349 + /* We just added some energy, so we need to renormalise */
1.350 + if (renormalize)
1.351 + renormalise_vector(X, N0<<LM, Q15ONE);
1.352 + } while (++c<C);
1.353 + }
1.354 +}
1.355 +
1.356 +static void intensity_stereo(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bandE, int bandID, int N)
1.357 +{
1.358 + int i = bandID;
1.359 + int j;
1.360 + opus_val16 a1, a2;
1.361 + opus_val16 left, right;
1.362 + opus_val16 norm;
1.363 +#ifdef FIXED_POINT
1.364 + int shift = celt_zlog2(MAX32(bandE[i], bandE[i+m->nbEBands]))-13;
1.365 +#endif
1.366 + left = VSHR32(bandE[i],shift);
1.367 + right = VSHR32(bandE[i+m->nbEBands],shift);
1.368 + norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
1.369 + a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
1.370 + a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);
1.371 + for (j=0;j<N;j++)
1.372 + {
1.373 + celt_norm r, l;
1.374 + l = X[j];
1.375 + r = Y[j];
1.376 + X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
1.377 + /* Side is not encoded, no need to calculate */
1.378 + }
1.379 +}
1.380 +
1.381 +static void stereo_split(celt_norm *X, celt_norm *Y, int N)
1.382 +{
1.383 + int j;
1.384 + for (j=0;j<N;j++)
1.385 + {
1.386 + celt_norm r, l;
1.387 + l = MULT16_16_Q15(QCONST16(.70710678f,15), X[j]);
1.388 + r = MULT16_16_Q15(QCONST16(.70710678f,15), Y[j]);
1.389 + X[j] = l+r;
1.390 + Y[j] = r-l;
1.391 + }
1.392 +}
1.393 +
1.394 +static void stereo_merge(celt_norm *X, celt_norm *Y, opus_val16 mid, int N)
1.395 +{
1.396 + int j;
1.397 + opus_val32 xp=0, side=0;
1.398 + opus_val32 El, Er;
1.399 + opus_val16 mid2;
1.400 +#ifdef FIXED_POINT
1.401 + int kl, kr;
1.402 +#endif
1.403 + opus_val32 t, lgain, rgain;
1.404 +
1.405 + /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
1.406 + dual_inner_prod(Y, X, Y, N, &xp, &side);
1.407 + /* Compensating for the mid normalization */
1.408 + xp = MULT16_32_Q15(mid, xp);
1.409 + /* mid and side are in Q15, not Q14 like X and Y */
1.410 + mid2 = SHR32(mid, 1);
1.411 + El = MULT16_16(mid2, mid2) + side - 2*xp;
1.412 + Er = MULT16_16(mid2, mid2) + side + 2*xp;
1.413 + if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28))
1.414 + {
1.415 + for (j=0;j<N;j++)
1.416 + Y[j] = X[j];
1.417 + return;
1.418 + }
1.419 +
1.420 +#ifdef FIXED_POINT
1.421 + kl = celt_ilog2(El)>>1;
1.422 + kr = celt_ilog2(Er)>>1;
1.423 +#endif
1.424 + t = VSHR32(El, (kl-7)<<1);
1.425 + lgain = celt_rsqrt_norm(t);
1.426 + t = VSHR32(Er, (kr-7)<<1);
1.427 + rgain = celt_rsqrt_norm(t);
1.428 +
1.429 +#ifdef FIXED_POINT
1.430 + if (kl < 7)
1.431 + kl = 7;
1.432 + if (kr < 7)
1.433 + kr = 7;
1.434 +#endif
1.435 +
1.436 + for (j=0;j<N;j++)
1.437 + {
1.438 + celt_norm r, l;
1.439 + /* Apply mid scaling (side is already scaled) */
1.440 + l = MULT16_16_Q15(mid, X[j]);
1.441 + r = Y[j];
1.442 + X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));
1.443 + Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));
1.444 + }
1.445 +}
1.446 +
1.447 +/* Decide whether we should spread the pulses in the current frame */
1.448 +int spreading_decision(const CELTMode *m, celt_norm *X, int *average,
1.449 + int last_decision, int *hf_average, int *tapset_decision, int update_hf,
1.450 + int end, int C, int M)
1.451 +{
1.452 + int i, c, N0;
1.453 + int sum = 0, nbBands=0;
1.454 + const opus_int16 * OPUS_RESTRICT eBands = m->eBands;
1.455 + int decision;
1.456 + int hf_sum=0;
1.457 +
1.458 + celt_assert(end>0);
1.459 +
1.460 + N0 = M*m->shortMdctSize;
1.461 +
1.462 + if (M*(eBands[end]-eBands[end-1]) <= 8)
1.463 + return SPREAD_NONE;
1.464 + c=0; do {
1.465 + for (i=0;i<end;i++)
1.466 + {
1.467 + int j, N, tmp=0;
1.468 + int tcount[3] = {0,0,0};
1.469 + celt_norm * OPUS_RESTRICT x = X+M*eBands[i]+c*N0;
1.470 + N = M*(eBands[i+1]-eBands[i]);
1.471 + if (N<=8)
1.472 + continue;
1.473 + /* Compute rough CDF of |x[j]| */
1.474 + for (j=0;j<N;j++)
1.475 + {
1.476 + opus_val32 x2N; /* Q13 */
1.477 +
1.478 + x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
1.479 + if (x2N < QCONST16(0.25f,13))
1.480 + tcount[0]++;
1.481 + if (x2N < QCONST16(0.0625f,13))
1.482 + tcount[1]++;
1.483 + if (x2N < QCONST16(0.015625f,13))
1.484 + tcount[2]++;
1.485 + }
1.486 +
1.487 + /* Only include four last bands (8 kHz and up) */
1.488 + if (i>m->nbEBands-4)
1.489 + hf_sum += 32*(tcount[1]+tcount[0])/N;
1.490 + tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
1.491 + sum += tmp*256;
1.492 + nbBands++;
1.493 + }
1.494 + } while (++c<C);
1.495 +
1.496 + if (update_hf)
1.497 + {
1.498 + if (hf_sum)
1.499 + hf_sum /= C*(4-m->nbEBands+end);
1.500 + *hf_average = (*hf_average+hf_sum)>>1;
1.501 + hf_sum = *hf_average;
1.502 + if (*tapset_decision==2)
1.503 + hf_sum += 4;
1.504 + else if (*tapset_decision==0)
1.505 + hf_sum -= 4;
1.506 + if (hf_sum > 22)
1.507 + *tapset_decision=2;
1.508 + else if (hf_sum > 18)
1.509 + *tapset_decision=1;
1.510 + else
1.511 + *tapset_decision=0;
1.512 + }
1.513 + /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/
1.514 + celt_assert(nbBands>0); /* end has to be non-zero */
1.515 + sum /= nbBands;
1.516 + /* Recursive averaging */
1.517 + sum = (sum+*average)>>1;
1.518 + *average = sum;
1.519 + /* Hysteresis */
1.520 + sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2;
1.521 + if (sum < 80)
1.522 + {
1.523 + decision = SPREAD_AGGRESSIVE;
1.524 + } else if (sum < 256)
1.525 + {
1.526 + decision = SPREAD_NORMAL;
1.527 + } else if (sum < 384)
1.528 + {
1.529 + decision = SPREAD_LIGHT;
1.530 + } else {
1.531 + decision = SPREAD_NONE;
1.532 + }
1.533 +#ifdef FUZZING
1.534 + decision = rand()&0x3;
1.535 + *tapset_decision=rand()%3;
1.536 +#endif
1.537 + return decision;
1.538 +}
1.539 +
1.540 +/* Indexing table for converting from natural Hadamard to ordery Hadamard
1.541 + This is essentially a bit-reversed Gray, on top of which we've added
1.542 + an inversion of the order because we want the DC at the end rather than
1.543 + the beginning. The lines are for N=2, 4, 8, 16 */
1.544 +static const int ordery_table[] = {
1.545 + 1, 0,
1.546 + 3, 0, 2, 1,
1.547 + 7, 0, 4, 3, 6, 1, 5, 2,
1.548 + 15, 0, 8, 7, 12, 3, 11, 4, 14, 1, 9, 6, 13, 2, 10, 5,
1.549 +};
1.550 +
1.551 +static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
1.552 +{
1.553 + int i,j;
1.554 + VARDECL(celt_norm, tmp);
1.555 + int N;
1.556 + SAVE_STACK;
1.557 + N = N0*stride;
1.558 + ALLOC(tmp, N, celt_norm);
1.559 + celt_assert(stride>0);
1.560 + if (hadamard)
1.561 + {
1.562 + const int *ordery = ordery_table+stride-2;
1.563 + for (i=0;i<stride;i++)
1.564 + {
1.565 + for (j=0;j<N0;j++)
1.566 + tmp[ordery[i]*N0+j] = X[j*stride+i];
1.567 + }
1.568 + } else {
1.569 + for (i=0;i<stride;i++)
1.570 + for (j=0;j<N0;j++)
1.571 + tmp[i*N0+j] = X[j*stride+i];
1.572 + }
1.573 + for (j=0;j<N;j++)
1.574 + X[j] = tmp[j];
1.575 + RESTORE_STACK;
1.576 +}
1.577 +
1.578 +static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
1.579 +{
1.580 + int i,j;
1.581 + VARDECL(celt_norm, tmp);
1.582 + int N;
1.583 + SAVE_STACK;
1.584 + N = N0*stride;
1.585 + ALLOC(tmp, N, celt_norm);
1.586 + if (hadamard)
1.587 + {
1.588 + const int *ordery = ordery_table+stride-2;
1.589 + for (i=0;i<stride;i++)
1.590 + for (j=0;j<N0;j++)
1.591 + tmp[j*stride+i] = X[ordery[i]*N0+j];
1.592 + } else {
1.593 + for (i=0;i<stride;i++)
1.594 + for (j=0;j<N0;j++)
1.595 + tmp[j*stride+i] = X[i*N0+j];
1.596 + }
1.597 + for (j=0;j<N;j++)
1.598 + X[j] = tmp[j];
1.599 + RESTORE_STACK;
1.600 +}
1.601 +
1.602 +void haar1(celt_norm *X, int N0, int stride)
1.603 +{
1.604 + int i, j;
1.605 + N0 >>= 1;
1.606 + for (i=0;i<stride;i++)
1.607 + for (j=0;j<N0;j++)
1.608 + {
1.609 + celt_norm tmp1, tmp2;
1.610 + tmp1 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*2*j+i]);
1.611 + tmp2 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*(2*j+1)+i]);
1.612 + X[stride*2*j+i] = tmp1 + tmp2;
1.613 + X[stride*(2*j+1)+i] = tmp1 - tmp2;
1.614 + }
1.615 +}
1.616 +
1.617 +static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo)
1.618 +{
1.619 + static const opus_int16 exp2_table8[8] =
1.620 + {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
1.621 + int qn, qb;
1.622 + int N2 = 2*N-1;
1.623 + if (stereo && N==2)
1.624 + N2--;
1.625 + /* The upper limit ensures that in a stereo split with itheta==16384, we'll
1.626 + always have enough bits left over to code at least one pulse in the
1.627 + side; otherwise it would collapse, since it doesn't get folded. */
1.628 + qb = IMIN(b-pulse_cap-(4<<BITRES), (b+N2*offset)/N2);
1.629 +
1.630 + qb = IMIN(8<<BITRES, qb);
1.631 +
1.632 + if (qb<(1<<BITRES>>1)) {
1.633 + qn = 1;
1.634 + } else {
1.635 + qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
1.636 + qn = (qn+1)>>1<<1;
1.637 + }
1.638 + celt_assert(qn <= 256);
1.639 + return qn;
1.640 +}
1.641 +
1.642 +struct band_ctx {
1.643 + int encode;
1.644 + const CELTMode *m;
1.645 + int i;
1.646 + int intensity;
1.647 + int spread;
1.648 + int tf_change;
1.649 + ec_ctx *ec;
1.650 + opus_int32 remaining_bits;
1.651 + const celt_ener *bandE;
1.652 + opus_uint32 seed;
1.653 +};
1.654 +
1.655 +struct split_ctx {
1.656 + int inv;
1.657 + int imid;
1.658 + int iside;
1.659 + int delta;
1.660 + int itheta;
1.661 + int qalloc;
1.662 +};
1.663 +
1.664 +static void compute_theta(struct band_ctx *ctx, struct split_ctx *sctx,
1.665 + celt_norm *X, celt_norm *Y, int N, int *b, int B, int B0,
1.666 + int LM,
1.667 + int stereo, int *fill)
1.668 +{
1.669 + int qn;
1.670 + int itheta=0;
1.671 + int delta;
1.672 + int imid, iside;
1.673 + int qalloc;
1.674 + int pulse_cap;
1.675 + int offset;
1.676 + opus_int32 tell;
1.677 + int inv=0;
1.678 + int encode;
1.679 + const CELTMode *m;
1.680 + int i;
1.681 + int intensity;
1.682 + ec_ctx *ec;
1.683 + const celt_ener *bandE;
1.684 +
1.685 + encode = ctx->encode;
1.686 + m = ctx->m;
1.687 + i = ctx->i;
1.688 + intensity = ctx->intensity;
1.689 + ec = ctx->ec;
1.690 + bandE = ctx->bandE;
1.691 +
1.692 + /* Decide on the resolution to give to the split parameter theta */
1.693 + pulse_cap = m->logN[i]+LM*(1<<BITRES);
1.694 + offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
1.695 + qn = compute_qn(N, *b, offset, pulse_cap, stereo);
1.696 + if (stereo && i>=intensity)
1.697 + qn = 1;
1.698 + if (encode)
1.699 + {
1.700 + /* theta is the atan() of the ratio between the (normalized)
1.701 + side and mid. With just that parameter, we can re-scale both
1.702 + mid and side because we know that 1) they have unit norm and
1.703 + 2) they are orthogonal. */
1.704 + itheta = stereo_itheta(X, Y, stereo, N);
1.705 + }
1.706 + tell = ec_tell_frac(ec);
1.707 + if (qn!=1)
1.708 + {
1.709 + if (encode)
1.710 + itheta = (itheta*qn+8192)>>14;
1.711 +
1.712 + /* Entropy coding of the angle. We use a uniform pdf for the
1.713 + time split, a step for stereo, and a triangular one for the rest. */
1.714 + if (stereo && N>2)
1.715 + {
1.716 + int p0 = 3;
1.717 + int x = itheta;
1.718 + int x0 = qn/2;
1.719 + int ft = p0*(x0+1) + x0;
1.720 + /* Use a probability of p0 up to itheta=8192 and then use 1 after */
1.721 + if (encode)
1.722 + {
1.723 + ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
1.724 + } else {
1.725 + int fs;
1.726 + fs=ec_decode(ec,ft);
1.727 + if (fs<(x0+1)*p0)
1.728 + x=fs/p0;
1.729 + else
1.730 + x=x0+1+(fs-(x0+1)*p0);
1.731 + 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);
1.732 + itheta = x;
1.733 + }
1.734 + } else if (B0>1 || stereo) {
1.735 + /* Uniform pdf */
1.736 + if (encode)
1.737 + ec_enc_uint(ec, itheta, qn+1);
1.738 + else
1.739 + itheta = ec_dec_uint(ec, qn+1);
1.740 + } else {
1.741 + int fs=1, ft;
1.742 + ft = ((qn>>1)+1)*((qn>>1)+1);
1.743 + if (encode)
1.744 + {
1.745 + int fl;
1.746 +
1.747 + fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
1.748 + fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
1.749 + ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
1.750 +
1.751 + ec_encode(ec, fl, fl+fs, ft);
1.752 + } else {
1.753 + /* Triangular pdf */
1.754 + int fl=0;
1.755 + int fm;
1.756 + fm = ec_decode(ec, ft);
1.757 +
1.758 + if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
1.759 + {
1.760 + itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;
1.761 + fs = itheta + 1;
1.762 + fl = itheta*(itheta + 1)>>1;
1.763 + }
1.764 + else
1.765 + {
1.766 + itheta = (2*(qn + 1)
1.767 + - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;
1.768 + fs = qn + 1 - itheta;
1.769 + fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
1.770 + }
1.771 +
1.772 + ec_dec_update(ec, fl, fl+fs, ft);
1.773 + }
1.774 + }
1.775 + itheta = (opus_int32)itheta*16384/qn;
1.776 + if (encode && stereo)
1.777 + {
1.778 + if (itheta==0)
1.779 + intensity_stereo(m, X, Y, bandE, i, N);
1.780 + else
1.781 + stereo_split(X, Y, N);
1.782 + }
1.783 + /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.
1.784 + Let's do that at higher complexity */
1.785 + } else if (stereo) {
1.786 + if (encode)
1.787 + {
1.788 + inv = itheta > 8192;
1.789 + if (inv)
1.790 + {
1.791 + int j;
1.792 + for (j=0;j<N;j++)
1.793 + Y[j] = -Y[j];
1.794 + }
1.795 + intensity_stereo(m, X, Y, bandE, i, N);
1.796 + }
1.797 + if (*b>2<<BITRES && ctx->remaining_bits > 2<<BITRES)
1.798 + {
1.799 + if (encode)
1.800 + ec_enc_bit_logp(ec, inv, 2);
1.801 + else
1.802 + inv = ec_dec_bit_logp(ec, 2);
1.803 + } else
1.804 + inv = 0;
1.805 + itheta = 0;
1.806 + }
1.807 + qalloc = ec_tell_frac(ec) - tell;
1.808 + *b -= qalloc;
1.809 +
1.810 + if (itheta == 0)
1.811 + {
1.812 + imid = 32767;
1.813 + iside = 0;
1.814 + *fill &= (1<<B)-1;
1.815 + delta = -16384;
1.816 + } else if (itheta == 16384)
1.817 + {
1.818 + imid = 0;
1.819 + iside = 32767;
1.820 + *fill &= ((1<<B)-1)<<B;
1.821 + delta = 16384;
1.822 + } else {
1.823 + imid = bitexact_cos((opus_int16)itheta);
1.824 + iside = bitexact_cos((opus_int16)(16384-itheta));
1.825 + /* This is the mid vs side allocation that minimizes squared error
1.826 + in that band. */
1.827 + delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid));
1.828 + }
1.829 +
1.830 + sctx->inv = inv;
1.831 + sctx->imid = imid;
1.832 + sctx->iside = iside;
1.833 + sctx->delta = delta;
1.834 + sctx->itheta = itheta;
1.835 + sctx->qalloc = qalloc;
1.836 +}
1.837 +static unsigned quant_band_n1(struct band_ctx *ctx, celt_norm *X, celt_norm *Y, int b,
1.838 + celt_norm *lowband_out)
1.839 +{
1.840 +#ifdef RESYNTH
1.841 + int resynth = 1;
1.842 +#else
1.843 + int resynth = !ctx->encode;
1.844 +#endif
1.845 + int c;
1.846 + int stereo;
1.847 + celt_norm *x = X;
1.848 + int encode;
1.849 + ec_ctx *ec;
1.850 +
1.851 + encode = ctx->encode;
1.852 + ec = ctx->ec;
1.853 +
1.854 + stereo = Y != NULL;
1.855 + c=0; do {
1.856 + int sign=0;
1.857 + if (ctx->remaining_bits>=1<<BITRES)
1.858 + {
1.859 + if (encode)
1.860 + {
1.861 + sign = x[0]<0;
1.862 + ec_enc_bits(ec, sign, 1);
1.863 + } else {
1.864 + sign = ec_dec_bits(ec, 1);
1.865 + }
1.866 + ctx->remaining_bits -= 1<<BITRES;
1.867 + b-=1<<BITRES;
1.868 + }
1.869 + if (resynth)
1.870 + x[0] = sign ? -NORM_SCALING : NORM_SCALING;
1.871 + x = Y;
1.872 + } while (++c<1+stereo);
1.873 + if (lowband_out)
1.874 + lowband_out[0] = SHR16(X[0],4);
1.875 + return 1;
1.876 +}
1.877 +
1.878 +/* This function is responsible for encoding and decoding a mono partition.
1.879 + It can split the band in two and transmit the energy difference with
1.880 + the two half-bands. It can be called recursively so bands can end up being
1.881 + split in 8 parts. */
1.882 +static unsigned quant_partition(struct band_ctx *ctx, celt_norm *X,
1.883 + int N, int b, int B, celt_norm *lowband,
1.884 + int LM,
1.885 + opus_val16 gain, int fill)
1.886 +{
1.887 + const unsigned char *cache;
1.888 + int q;
1.889 + int curr_bits;
1.890 + int imid=0, iside=0;
1.891 + int B0=B;
1.892 + opus_val16 mid=0, side=0;
1.893 + unsigned cm=0;
1.894 +#ifdef RESYNTH
1.895 + int resynth = 1;
1.896 +#else
1.897 + int resynth = !ctx->encode;
1.898 +#endif
1.899 + celt_norm *Y=NULL;
1.900 + int encode;
1.901 + const CELTMode *m;
1.902 + int i;
1.903 + int spread;
1.904 + ec_ctx *ec;
1.905 +
1.906 + encode = ctx->encode;
1.907 + m = ctx->m;
1.908 + i = ctx->i;
1.909 + spread = ctx->spread;
1.910 + ec = ctx->ec;
1.911 +
1.912 + /* If we need 1.5 more bit than we can produce, split the band in two. */
1.913 + cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
1.914 + if (LM != -1 && b > cache[cache[0]]+12 && N>2)
1.915 + {
1.916 + int mbits, sbits, delta;
1.917 + int itheta;
1.918 + int qalloc;
1.919 + struct split_ctx sctx;
1.920 + celt_norm *next_lowband2=NULL;
1.921 + opus_int32 rebalance;
1.922 +
1.923 + N >>= 1;
1.924 + Y = X+N;
1.925 + LM -= 1;
1.926 + if (B==1)
1.927 + fill = (fill&1)|(fill<<1);
1.928 + B = (B+1)>>1;
1.929 +
1.930 + compute_theta(ctx, &sctx, X, Y, N, &b, B, B0,
1.931 + LM, 0, &fill);
1.932 + imid = sctx.imid;
1.933 + iside = sctx.iside;
1.934 + delta = sctx.delta;
1.935 + itheta = sctx.itheta;
1.936 + qalloc = sctx.qalloc;
1.937 +#ifdef FIXED_POINT
1.938 + mid = imid;
1.939 + side = iside;
1.940 +#else
1.941 + mid = (1.f/32768)*imid;
1.942 + side = (1.f/32768)*iside;
1.943 +#endif
1.944 +
1.945 + /* Give more bits to low-energy MDCTs than they would otherwise deserve */
1.946 + if (B0>1 && (itheta&0x3fff))
1.947 + {
1.948 + if (itheta > 8192)
1.949 + /* Rough approximation for pre-echo masking */
1.950 + delta -= delta>>(4-LM);
1.951 + else
1.952 + /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
1.953 + delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
1.954 + }
1.955 + mbits = IMAX(0, IMIN(b, (b-delta)/2));
1.956 + sbits = b-mbits;
1.957 + ctx->remaining_bits -= qalloc;
1.958 +
1.959 + if (lowband)
1.960 + next_lowband2 = lowband+N; /* >32-bit split case */
1.961 +
1.962 + rebalance = ctx->remaining_bits;
1.963 + if (mbits >= sbits)
1.964 + {
1.965 + cm = quant_partition(ctx, X, N, mbits, B,
1.966 + lowband, LM,
1.967 + MULT16_16_P15(gain,mid), fill);
1.968 + rebalance = mbits - (rebalance-ctx->remaining_bits);
1.969 + if (rebalance > 3<<BITRES && itheta!=0)
1.970 + sbits += rebalance - (3<<BITRES);
1.971 + cm |= quant_partition(ctx, Y, N, sbits, B,
1.972 + next_lowband2, LM,
1.973 + MULT16_16_P15(gain,side), fill>>B)<<(B0>>1);
1.974 + } else {
1.975 + cm = quant_partition(ctx, Y, N, sbits, B,
1.976 + next_lowband2, LM,
1.977 + MULT16_16_P15(gain,side), fill>>B)<<(B0>>1);
1.978 + rebalance = sbits - (rebalance-ctx->remaining_bits);
1.979 + if (rebalance > 3<<BITRES && itheta!=16384)
1.980 + mbits += rebalance - (3<<BITRES);
1.981 + cm |= quant_partition(ctx, X, N, mbits, B,
1.982 + lowband, LM,
1.983 + MULT16_16_P15(gain,mid), fill);
1.984 + }
1.985 + } else {
1.986 + /* This is the basic no-split case */
1.987 + q = bits2pulses(m, i, LM, b);
1.988 + curr_bits = pulses2bits(m, i, LM, q);
1.989 + ctx->remaining_bits -= curr_bits;
1.990 +
1.991 + /* Ensures we can never bust the budget */
1.992 + while (ctx->remaining_bits < 0 && q > 0)
1.993 + {
1.994 + ctx->remaining_bits += curr_bits;
1.995 + q--;
1.996 + curr_bits = pulses2bits(m, i, LM, q);
1.997 + ctx->remaining_bits -= curr_bits;
1.998 + }
1.999 +
1.1000 + if (q!=0)
1.1001 + {
1.1002 + int K = get_pulses(q);
1.1003 +
1.1004 + /* Finally do the actual quantization */
1.1005 + if (encode)
1.1006 + {
1.1007 + cm = alg_quant(X, N, K, spread, B, ec
1.1008 +#ifdef RESYNTH
1.1009 + , gain
1.1010 +#endif
1.1011 + );
1.1012 + } else {
1.1013 + cm = alg_unquant(X, N, K, spread, B, ec, gain);
1.1014 + }
1.1015 + } else {
1.1016 + /* If there's no pulse, fill the band anyway */
1.1017 + int j;
1.1018 + if (resynth)
1.1019 + {
1.1020 + unsigned cm_mask;
1.1021 + /* B can be as large as 16, so this shift might overflow an int on a
1.1022 + 16-bit platform; use a long to get defined behavior.*/
1.1023 + cm_mask = (unsigned)(1UL<<B)-1;
1.1024 + fill &= cm_mask;
1.1025 + if (!fill)
1.1026 + {
1.1027 + for (j=0;j<N;j++)
1.1028 + X[j] = 0;
1.1029 + } else {
1.1030 + if (lowband == NULL)
1.1031 + {
1.1032 + /* Noise */
1.1033 + for (j=0;j<N;j++)
1.1034 + {
1.1035 + ctx->seed = celt_lcg_rand(ctx->seed);
1.1036 + X[j] = (celt_norm)((opus_int32)ctx->seed>>20);
1.1037 + }
1.1038 + cm = cm_mask;
1.1039 + } else {
1.1040 + /* Folded spectrum */
1.1041 + for (j=0;j<N;j++)
1.1042 + {
1.1043 + opus_val16 tmp;
1.1044 + ctx->seed = celt_lcg_rand(ctx->seed);
1.1045 + /* About 48 dB below the "normal" folding level */
1.1046 + tmp = QCONST16(1.0f/256, 10);
1.1047 + tmp = (ctx->seed)&0x8000 ? tmp : -tmp;
1.1048 + X[j] = lowband[j]+tmp;
1.1049 + }
1.1050 + cm = fill;
1.1051 + }
1.1052 + renormalise_vector(X, N, gain);
1.1053 + }
1.1054 + }
1.1055 + }
1.1056 + }
1.1057 +
1.1058 + return cm;
1.1059 +}
1.1060 +
1.1061 +
1.1062 +/* This function is responsible for encoding and decoding a band for the mono case. */
1.1063 +static unsigned quant_band(struct band_ctx *ctx, celt_norm *X,
1.1064 + int N, int b, int B, celt_norm *lowband,
1.1065 + int LM, celt_norm *lowband_out,
1.1066 + opus_val16 gain, celt_norm *lowband_scratch, int fill)
1.1067 +{
1.1068 + int N0=N;
1.1069 + int N_B=N;
1.1070 + int N_B0;
1.1071 + int B0=B;
1.1072 + int time_divide=0;
1.1073 + int recombine=0;
1.1074 + int longBlocks;
1.1075 + unsigned cm=0;
1.1076 +#ifdef RESYNTH
1.1077 + int resynth = 1;
1.1078 +#else
1.1079 + int resynth = !ctx->encode;
1.1080 +#endif
1.1081 + int k;
1.1082 + int encode;
1.1083 + int tf_change;
1.1084 +
1.1085 + encode = ctx->encode;
1.1086 + tf_change = ctx->tf_change;
1.1087 +
1.1088 + longBlocks = B0==1;
1.1089 +
1.1090 + N_B /= B;
1.1091 +
1.1092 + /* Special case for one sample */
1.1093 + if (N==1)
1.1094 + {
1.1095 + return quant_band_n1(ctx, X, NULL, b, lowband_out);
1.1096 + }
1.1097 +
1.1098 + if (tf_change>0)
1.1099 + recombine = tf_change;
1.1100 + /* Band recombining to increase frequency resolution */
1.1101 +
1.1102 + if (lowband_scratch && lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
1.1103 + {
1.1104 + int j;
1.1105 + for (j=0;j<N;j++)
1.1106 + lowband_scratch[j] = lowband[j];
1.1107 + lowband = lowband_scratch;
1.1108 + }
1.1109 +
1.1110 + for (k=0;k<recombine;k++)
1.1111 + {
1.1112 + static const unsigned char bit_interleave_table[16]={
1.1113 + 0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
1.1114 + };
1.1115 + if (encode)
1.1116 + haar1(X, N>>k, 1<<k);
1.1117 + if (lowband)
1.1118 + haar1(lowband, N>>k, 1<<k);
1.1119 + fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
1.1120 + }
1.1121 + B>>=recombine;
1.1122 + N_B<<=recombine;
1.1123 +
1.1124 + /* Increasing the time resolution */
1.1125 + while ((N_B&1) == 0 && tf_change<0)
1.1126 + {
1.1127 + if (encode)
1.1128 + haar1(X, N_B, B);
1.1129 + if (lowband)
1.1130 + haar1(lowband, N_B, B);
1.1131 + fill |= fill<<B;
1.1132 + B <<= 1;
1.1133 + N_B >>= 1;
1.1134 + time_divide++;
1.1135 + tf_change++;
1.1136 + }
1.1137 + B0=B;
1.1138 + N_B0 = N_B;
1.1139 +
1.1140 + /* Reorganize the samples in time order instead of frequency order */
1.1141 + if (B0>1)
1.1142 + {
1.1143 + if (encode)
1.1144 + deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1.1145 + if (lowband)
1.1146 + deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
1.1147 + }
1.1148 +
1.1149 + cm = quant_partition(ctx, X, N, b, B, lowband,
1.1150 + LM, gain, fill);
1.1151 +
1.1152 + /* This code is used by the decoder and by the resynthesis-enabled encoder */
1.1153 + if (resynth)
1.1154 + {
1.1155 + /* Undo the sample reorganization going from time order to frequency order */
1.1156 + if (B0>1)
1.1157 + interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1.1158 +
1.1159 + /* Undo time-freq changes that we did earlier */
1.1160 + N_B = N_B0;
1.1161 + B = B0;
1.1162 + for (k=0;k<time_divide;k++)
1.1163 + {
1.1164 + B >>= 1;
1.1165 + N_B <<= 1;
1.1166 + cm |= cm>>B;
1.1167 + haar1(X, N_B, B);
1.1168 + }
1.1169 +
1.1170 + for (k=0;k<recombine;k++)
1.1171 + {
1.1172 + static const unsigned char bit_deinterleave_table[16]={
1.1173 + 0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1.1174 + 0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1.1175 + };
1.1176 + cm = bit_deinterleave_table[cm];
1.1177 + haar1(X, N0>>k, 1<<k);
1.1178 + }
1.1179 + B<<=recombine;
1.1180 +
1.1181 + /* Scale output for later folding */
1.1182 + if (lowband_out)
1.1183 + {
1.1184 + int j;
1.1185 + opus_val16 n;
1.1186 + n = celt_sqrt(SHL32(EXTEND32(N0),22));
1.1187 + for (j=0;j<N0;j++)
1.1188 + lowband_out[j] = MULT16_16_Q15(n,X[j]);
1.1189 + }
1.1190 + cm &= (1<<B)-1;
1.1191 + }
1.1192 + return cm;
1.1193 +}
1.1194 +
1.1195 +
1.1196 +/* This function is responsible for encoding and decoding a band for the stereo case. */
1.1197 +static unsigned quant_band_stereo(struct band_ctx *ctx, celt_norm *X, celt_norm *Y,
1.1198 + int N, int b, int B, celt_norm *lowband,
1.1199 + int LM, celt_norm *lowband_out,
1.1200 + celt_norm *lowband_scratch, int fill)
1.1201 +{
1.1202 + int imid=0, iside=0;
1.1203 + int inv = 0;
1.1204 + opus_val16 mid=0, side=0;
1.1205 + unsigned cm=0;
1.1206 +#ifdef RESYNTH
1.1207 + int resynth = 1;
1.1208 +#else
1.1209 + int resynth = !ctx->encode;
1.1210 +#endif
1.1211 + int mbits, sbits, delta;
1.1212 + int itheta;
1.1213 + int qalloc;
1.1214 + struct split_ctx sctx;
1.1215 + int orig_fill;
1.1216 + int encode;
1.1217 + ec_ctx *ec;
1.1218 +
1.1219 + encode = ctx->encode;
1.1220 + ec = ctx->ec;
1.1221 +
1.1222 + /* Special case for one sample */
1.1223 + if (N==1)
1.1224 + {
1.1225 + return quant_band_n1(ctx, X, Y, b, lowband_out);
1.1226 + }
1.1227 +
1.1228 + orig_fill = fill;
1.1229 +
1.1230 + compute_theta(ctx, &sctx, X, Y, N, &b, B, B,
1.1231 + LM, 1, &fill);
1.1232 + inv = sctx.inv;
1.1233 + imid = sctx.imid;
1.1234 + iside = sctx.iside;
1.1235 + delta = sctx.delta;
1.1236 + itheta = sctx.itheta;
1.1237 + qalloc = sctx.qalloc;
1.1238 +#ifdef FIXED_POINT
1.1239 + mid = imid;
1.1240 + side = iside;
1.1241 +#else
1.1242 + mid = (1.f/32768)*imid;
1.1243 + side = (1.f/32768)*iside;
1.1244 +#endif
1.1245 +
1.1246 + /* This is a special case for N=2 that only works for stereo and takes
1.1247 + advantage of the fact that mid and side are orthogonal to encode
1.1248 + the side with just one bit. */
1.1249 + if (N==2)
1.1250 + {
1.1251 + int c;
1.1252 + int sign=0;
1.1253 + celt_norm *x2, *y2;
1.1254 + mbits = b;
1.1255 + sbits = 0;
1.1256 + /* Only need one bit for the side. */
1.1257 + if (itheta != 0 && itheta != 16384)
1.1258 + sbits = 1<<BITRES;
1.1259 + mbits -= sbits;
1.1260 + c = itheta > 8192;
1.1261 + ctx->remaining_bits -= qalloc+sbits;
1.1262 +
1.1263 + x2 = c ? Y : X;
1.1264 + y2 = c ? X : Y;
1.1265 + if (sbits)
1.1266 + {
1.1267 + if (encode)
1.1268 + {
1.1269 + /* Here we only need to encode a sign for the side. */
1.1270 + sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
1.1271 + ec_enc_bits(ec, sign, 1);
1.1272 + } else {
1.1273 + sign = ec_dec_bits(ec, 1);
1.1274 + }
1.1275 + }
1.1276 + sign = 1-2*sign;
1.1277 + /* We use orig_fill here because we want to fold the side, but if
1.1278 + itheta==16384, we'll have cleared the low bits of fill. */
1.1279 + cm = quant_band(ctx, x2, N, mbits, B, lowband,
1.1280 + LM, lowband_out, Q15ONE, lowband_scratch, orig_fill);
1.1281 + /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
1.1282 + and there's no need to worry about mixing with the other channel. */
1.1283 + y2[0] = -sign*x2[1];
1.1284 + y2[1] = sign*x2[0];
1.1285 + if (resynth)
1.1286 + {
1.1287 + celt_norm tmp;
1.1288 + X[0] = MULT16_16_Q15(mid, X[0]);
1.1289 + X[1] = MULT16_16_Q15(mid, X[1]);
1.1290 + Y[0] = MULT16_16_Q15(side, Y[0]);
1.1291 + Y[1] = MULT16_16_Q15(side, Y[1]);
1.1292 + tmp = X[0];
1.1293 + X[0] = SUB16(tmp,Y[0]);
1.1294 + Y[0] = ADD16(tmp,Y[0]);
1.1295 + tmp = X[1];
1.1296 + X[1] = SUB16(tmp,Y[1]);
1.1297 + Y[1] = ADD16(tmp,Y[1]);
1.1298 + }
1.1299 + } else {
1.1300 + /* "Normal" split code */
1.1301 + opus_int32 rebalance;
1.1302 +
1.1303 + mbits = IMAX(0, IMIN(b, (b-delta)/2));
1.1304 + sbits = b-mbits;
1.1305 + ctx->remaining_bits -= qalloc;
1.1306 +
1.1307 + rebalance = ctx->remaining_bits;
1.1308 + if (mbits >= sbits)
1.1309 + {
1.1310 + /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1.1311 + mid for folding later. */
1.1312 + cm = quant_band(ctx, X, N, mbits, B,
1.1313 + lowband, LM, lowband_out,
1.1314 + Q15ONE, lowband_scratch, fill);
1.1315 + rebalance = mbits - (rebalance-ctx->remaining_bits);
1.1316 + if (rebalance > 3<<BITRES && itheta!=0)
1.1317 + sbits += rebalance - (3<<BITRES);
1.1318 +
1.1319 + /* For a stereo split, the high bits of fill are always zero, so no
1.1320 + folding will be done to the side. */
1.1321 + cm |= quant_band(ctx, Y, N, sbits, B,
1.1322 + NULL, LM, NULL,
1.1323 + side, NULL, fill>>B);
1.1324 + } else {
1.1325 + /* For a stereo split, the high bits of fill are always zero, so no
1.1326 + folding will be done to the side. */
1.1327 + cm = quant_band(ctx, Y, N, sbits, B,
1.1328 + NULL, LM, NULL,
1.1329 + side, NULL, fill>>B);
1.1330 + rebalance = sbits - (rebalance-ctx->remaining_bits);
1.1331 + if (rebalance > 3<<BITRES && itheta!=16384)
1.1332 + mbits += rebalance - (3<<BITRES);
1.1333 + /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1.1334 + mid for folding later. */
1.1335 + cm |= quant_band(ctx, X, N, mbits, B,
1.1336 + lowband, LM, lowband_out,
1.1337 + Q15ONE, lowband_scratch, fill);
1.1338 + }
1.1339 + }
1.1340 +
1.1341 +
1.1342 + /* This code is used by the decoder and by the resynthesis-enabled encoder */
1.1343 + if (resynth)
1.1344 + {
1.1345 + if (N!=2)
1.1346 + stereo_merge(X, Y, mid, N);
1.1347 + if (inv)
1.1348 + {
1.1349 + int j;
1.1350 + for (j=0;j<N;j++)
1.1351 + Y[j] = -Y[j];
1.1352 + }
1.1353 + }
1.1354 + return cm;
1.1355 +}
1.1356 +
1.1357 +
1.1358 +void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1.1359 + celt_norm *X_, celt_norm *Y_, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1.1360 + int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res,
1.1361 + opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed)
1.1362 +{
1.1363 + int i;
1.1364 + opus_int32 remaining_bits;
1.1365 + const opus_int16 * OPUS_RESTRICT eBands = m->eBands;
1.1366 + celt_norm * OPUS_RESTRICT norm, * OPUS_RESTRICT norm2;
1.1367 + VARDECL(celt_norm, _norm);
1.1368 + celt_norm *lowband_scratch;
1.1369 + int B;
1.1370 + int M;
1.1371 + int lowband_offset;
1.1372 + int update_lowband = 1;
1.1373 + int C = Y_ != NULL ? 2 : 1;
1.1374 + int norm_offset;
1.1375 +#ifdef RESYNTH
1.1376 + int resynth = 1;
1.1377 +#else
1.1378 + int resynth = !encode;
1.1379 +#endif
1.1380 + struct band_ctx ctx;
1.1381 + SAVE_STACK;
1.1382 +
1.1383 + M = 1<<LM;
1.1384 + B = shortBlocks ? M : 1;
1.1385 + norm_offset = M*eBands[start];
1.1386 + /* No need to allocate norm for the last band because we don't need an
1.1387 + output in that band. */
1.1388 + ALLOC(_norm, C*(M*eBands[m->nbEBands-1]-norm_offset), celt_norm);
1.1389 + norm = _norm;
1.1390 + norm2 = norm + M*eBands[m->nbEBands-1]-norm_offset;
1.1391 + /* We can use the last band as scratch space because we don't need that
1.1392 + scratch space for the last band. */
1.1393 + lowband_scratch = X_+M*eBands[m->nbEBands-1];
1.1394 +
1.1395 + lowband_offset = 0;
1.1396 + ctx.bandE = bandE;
1.1397 + ctx.ec = ec;
1.1398 + ctx.encode = encode;
1.1399 + ctx.intensity = intensity;
1.1400 + ctx.m = m;
1.1401 + ctx.seed = *seed;
1.1402 + ctx.spread = spread;
1.1403 + for (i=start;i<end;i++)
1.1404 + {
1.1405 + opus_int32 tell;
1.1406 + int b;
1.1407 + int N;
1.1408 + opus_int32 curr_balance;
1.1409 + int effective_lowband=-1;
1.1410 + celt_norm * OPUS_RESTRICT X, * OPUS_RESTRICT Y;
1.1411 + int tf_change=0;
1.1412 + unsigned x_cm;
1.1413 + unsigned y_cm;
1.1414 + int last;
1.1415 +
1.1416 + ctx.i = i;
1.1417 + last = (i==end-1);
1.1418 +
1.1419 + X = X_+M*eBands[i];
1.1420 + if (Y_!=NULL)
1.1421 + Y = Y_+M*eBands[i];
1.1422 + else
1.1423 + Y = NULL;
1.1424 + N = M*eBands[i+1]-M*eBands[i];
1.1425 + tell = ec_tell_frac(ec);
1.1426 +
1.1427 + /* Compute how many bits we want to allocate to this band */
1.1428 + if (i != start)
1.1429 + balance -= tell;
1.1430 + remaining_bits = total_bits-tell-1;
1.1431 + ctx.remaining_bits = remaining_bits;
1.1432 + if (i <= codedBands-1)
1.1433 + {
1.1434 + curr_balance = balance / IMIN(3, codedBands-i);
1.1435 + b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1.1436 + } else {
1.1437 + b = 0;
1.1438 + }
1.1439 +
1.1440 + if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1.1441 + lowband_offset = i;
1.1442 +
1.1443 + tf_change = tf_res[i];
1.1444 + ctx.tf_change = tf_change;
1.1445 + if (i>=m->effEBands)
1.1446 + {
1.1447 + X=norm;
1.1448 + if (Y_!=NULL)
1.1449 + Y = norm;
1.1450 + lowband_scratch = NULL;
1.1451 + }
1.1452 + if (i==end-1)
1.1453 + lowband_scratch = NULL;
1.1454 +
1.1455 + /* Get a conservative estimate of the collapse_mask's for the bands we're
1.1456 + going to be folding from. */
1.1457 + if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1.1458 + {
1.1459 + int fold_start;
1.1460 + int fold_end;
1.1461 + int fold_i;
1.1462 + /* This ensures we never repeat spectral content within one band */
1.1463 + effective_lowband = IMAX(0, M*eBands[lowband_offset]-norm_offset-N);
1.1464 + fold_start = lowband_offset;
1.1465 + while(M*eBands[--fold_start] > effective_lowband+norm_offset);
1.1466 + fold_end = lowband_offset-1;
1.1467 + while(M*eBands[++fold_end] < effective_lowband+norm_offset+N);
1.1468 + x_cm = y_cm = 0;
1.1469 + fold_i = fold_start; do {
1.1470 + x_cm |= collapse_masks[fold_i*C+0];
1.1471 + y_cm |= collapse_masks[fold_i*C+C-1];
1.1472 + } while (++fold_i<fold_end);
1.1473 + }
1.1474 + /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1.1475 + always) be non-zero. */
1.1476 + else
1.1477 + x_cm = y_cm = (1<<B)-1;
1.1478 +
1.1479 + if (dual_stereo && i==intensity)
1.1480 + {
1.1481 + int j;
1.1482 +
1.1483 + /* Switch off dual stereo to do intensity. */
1.1484 + dual_stereo = 0;
1.1485 + if (resynth)
1.1486 + for (j=0;j<M*eBands[i]-norm_offset;j++)
1.1487 + norm[j] = HALF32(norm[j]+norm2[j]);
1.1488 + }
1.1489 + if (dual_stereo)
1.1490 + {
1.1491 + x_cm = quant_band(&ctx, X, N, b/2, B,
1.1492 + effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1.1493 + last?NULL:norm+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, x_cm);
1.1494 + y_cm = quant_band(&ctx, Y, N, b/2, B,
1.1495 + effective_lowband != -1 ? norm2+effective_lowband : NULL, LM,
1.1496 + last?NULL:norm2+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, y_cm);
1.1497 + } else {
1.1498 + if (Y!=NULL)
1.1499 + {
1.1500 + x_cm = quant_band_stereo(&ctx, X, Y, N, b, B,
1.1501 + effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1.1502 + last?NULL:norm+M*eBands[i]-norm_offset, lowband_scratch, x_cm|y_cm);
1.1503 + } else {
1.1504 + x_cm = quant_band(&ctx, X, N, b, B,
1.1505 + effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1.1506 + last?NULL:norm+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, x_cm|y_cm);
1.1507 + }
1.1508 + y_cm = x_cm;
1.1509 + }
1.1510 + collapse_masks[i*C+0] = (unsigned char)x_cm;
1.1511 + collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1.1512 + balance += pulses[i] + tell;
1.1513 +
1.1514 + /* Update the folding position only as long as we have 1 bit/sample depth. */
1.1515 + update_lowband = b>(N<<BITRES);
1.1516 + }
1.1517 + *seed = ctx.seed;
1.1518 +
1.1519 + RESTORE_STACK;
1.1520 +}
1.1521 +
