opus_pvq.c 29.7 KB
Newer Older
1
/*
2 3 4
 * Copyright (c) 2007-2008 CSIRO
 * Copyright (c) 2007-2009 Xiph.Org Foundation
 * Copyright (c) 2008-2009 Gregory Maxwell
5 6
 * Copyright (c) 2012 Andrew D'Addesio
 * Copyright (c) 2013-2014 Mozilla Corporation
7
 * Copyright (c) 2017 Rostislav Pehlivanov <atomnuker@gmail.com>
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
 *
 * This file is part of FFmpeg.
 *
 * FFmpeg is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * FFmpeg is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with FFmpeg; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */

#include "opustab.h"
#include "opus_pvq.h"

#define CELT_PVQ_U(n, k) (ff_celt_pvq_u_row[FFMIN(n, k)][FFMAX(n, k)])
#define CELT_PVQ_V(n, k) (CELT_PVQ_U(n, k) + CELT_PVQ_U(n, (k) + 1))

static inline int16_t celt_cos(int16_t x)
{
    x = (MUL16(x, x) + 4096) >> 13;
    x = (32767-x) + ROUND_MUL16(x, (-7651 + ROUND_MUL16(x, (8277 + ROUND_MUL16(-626, x)))));
36
    return x + 1;
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
}

static inline int celt_log2tan(int isin, int icos)
{
    int lc, ls;
    lc = opus_ilog(icos);
    ls = opus_ilog(isin);
    icos <<= 15 - lc;
    isin <<= 15 - ls;
    return (ls << 11) - (lc << 11) +
           ROUND_MUL16(isin, ROUND_MUL16(isin, -2597) + 7932) -
           ROUND_MUL16(icos, ROUND_MUL16(icos, -2597) + 7932);
}

static inline int celt_bits2pulses(const uint8_t *cache, int bits)
{
    // TODO: Find the size of cache and make it into an array in the parameters list
    int i, low = 0, high;

    high = cache[0];
    bits--;

    for (i = 0; i < 6; i++) {
        int center = (low + high + 1) >> 1;
        if (cache[center] >= bits)
            high = center;
        else
            low = center;
    }

    return (bits - (low == 0 ? -1 : cache[low]) <= cache[high] - bits) ? low : high;
}

static inline int celt_pulses2bits(const uint8_t *cache, int pulses)
{
    // TODO: Find the size of cache and make it into an array in the parameters list
   return (pulses == 0) ? 0 : cache[pulses] + 1;
}

static inline void celt_normalize_residual(const int * av_restrict iy, float * av_restrict X,
                                           int N, float g)
{
    int i;
    for (i = 0; i < N; i++)
        X[i] = g * iy[i];
}

84 85
static void celt_exp_rotation_impl(float *X, uint32_t len, uint32_t stride,
                                   float c, float s)
86 87 88 89 90 91
{
    float *Xptr;
    int i;

    Xptr = X;
    for (i = 0; i < len - stride; i++) {
92 93
        float x1     = Xptr[0];
        float x2     = Xptr[stride];
94 95 96 97 98 99
        Xptr[stride] = c * x2 + s * x1;
        *Xptr++      = c * x1 - s * x2;
    }

    Xptr = &X[len - 2 * stride - 1];
    for (i = len - 2 * stride - 1; i >= 0; i--) {
100 101
        float x1     = Xptr[0];
        float x2     = Xptr[stride];
102 103 104 105 106 107 108
        Xptr[stride] = c * x2 + s * x1;
        *Xptr--      = c * x1 - s * x2;
    }
}

static inline void celt_exp_rotation(float *X, uint32_t len,
                                     uint32_t stride, uint32_t K,
109
                                     enum CeltSpread spread, const int encode)
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
{
    uint32_t stride2 = 0;
    float c, s;
    float gain, theta;
    int i;

    if (2*K >= len || spread == CELT_SPREAD_NONE)
        return;

    gain = (float)len / (len + (20 - 5*spread) * K);
    theta = M_PI * gain * gain / 4;

    c = cosf(theta);
    s = sinf(theta);

    if (len >= stride << 3) {
        stride2 = 1;
        /* This is just a simple (equivalent) way of computing sqrt(len/stride) with rounding.
        It's basically incrementing long as (stride2+0.5)^2 < len/stride. */
        while ((stride2 * stride2 + stride2) * stride + (stride >> 2) < len)
            stride2++;
    }

    len /= stride;
    for (i = 0; i < stride; i++) {
135 136 137 138 139 140 141 142 143
        if (encode) {
            celt_exp_rotation_impl(X + i * len, len, 1, c, -s);
            if (stride2)
                celt_exp_rotation_impl(X + i * len, len, stride2, s, -c);
        } else {
            if (stride2)
                celt_exp_rotation_impl(X + i * len, len, stride2, s, c);
            celt_exp_rotation_impl(X + i * len, len, 1, c, s);
        }
144 145 146 147 148
    }
}

static inline uint32_t celt_extract_collapse_mask(const int *iy, uint32_t N, uint32_t B)
{
149 150
    int i, j, N0 = N / B;
    uint32_t collapse_mask = 0;
151 152 153 154 155 156

    if (B <= 1)
        return 1;

    for (i = 0; i < B; i++)
        for (j = 0; j < N0; j++)
157
            collapse_mask |= (!!iy[i*N0+j]) << i;
158 159 160 161 162 163 164 165 166
    return collapse_mask;
}

static inline void celt_stereo_merge(float *X, float *Y, float mid, int N)
{
    int i;
    float xp = 0, side = 0;
    float E[2];
    float mid2;
167
    float gain[2];
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185

    /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
    for (i = 0; i < N; i++) {
        xp   += X[i] * Y[i];
        side += Y[i] * Y[i];
    }

    /* Compensating for the mid normalization */
    xp *= mid;
    mid2 = mid;
    E[0] = mid2 * mid2 + side - 2 * xp;
    E[1] = mid2 * mid2 + side + 2 * xp;
    if (E[0] < 6e-4f || E[1] < 6e-4f) {
        for (i = 0; i < N; i++)
            Y[i] = X[i];
        return;
    }

186 187
    gain[0] = 1.0f / sqrtf(E[0]);
    gain[1] = 1.0f / sqrtf(E[1]);
188 189 190 191 192 193 194 195 196 197 198 199 200 201

    for (i = 0; i < N; i++) {
        float value[2];
        /* Apply mid scaling (side is already scaled) */
        value[0] = mid * X[i];
        value[1] = Y[i];
        X[i] = gain[0] * (value[0] - value[1]);
        Y[i] = gain[1] * (value[0] + value[1]);
    }
}

static void celt_interleave_hadamard(float *tmp, float *X, int N0,
                                     int stride, int hadamard)
{
202 203
    int i, j, N = N0*stride;
    const uint8_t *order = &ff_celt_hadamard_order[hadamard ? stride - 2 : 30];
204

205 206 207
    for (i = 0; i < stride; i++)
        for (j = 0; j < N0; j++)
            tmp[j*stride+i] = X[order[i]*N0+j];
208

209
    memcpy(X, tmp, N*sizeof(float));
210 211 212 213 214
}

static void celt_deinterleave_hadamard(float *tmp, float *X, int N0,
                                       int stride, int hadamard)
{
215 216
    int i, j, N = N0*stride;
    const uint8_t *order = &ff_celt_hadamard_order[hadamard ? stride - 2 : 30];
217

218 219 220
    for (i = 0; i < stride; i++)
        for (j = 0; j < N0; j++)
            tmp[order[i]*N0+j] = X[j*stride+i];
221

222
    memcpy(X, tmp, N*sizeof(float));
223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
}

static void celt_haar1(float *X, int N0, int stride)
{
    int i, j;
    N0 >>= 1;
    for (i = 0; i < stride; i++) {
        for (j = 0; j < N0; j++) {
            float x0 = X[stride * (2 * j + 0) + i];
            float x1 = X[stride * (2 * j + 1) + i];
            X[stride * (2 * j + 0) + i] = (x0 + x1) * M_SQRT1_2;
            X[stride * (2 * j + 1) + i] = (x0 - x1) * M_SQRT1_2;
        }
    }
}

static inline int celt_compute_qn(int N, int b, int offset, int pulse_cap,
240
                                  int stereo)
241 242 243
{
    int qn, qb;
    int N2 = 2 * N - 1;
244
    if (stereo && N == 2)
245 246 247 248 249 250 251 252 253 254
        N2--;

    /* The upper limit ensures that in a stereo split with itheta==16384, we'll
     * always have enough bits left over to code at least one pulse in the
     * side; otherwise it would collapse, since it doesn't get folded. */
    qb = FFMIN3(b - pulse_cap - (4 << 3), (b + N2 * offset) / N2, 8 << 3);
    qn = (qb < (1 << 3 >> 1)) ? 1 : ((ff_celt_qn_exp2[qb & 0x7] >> (14 - (qb >> 3))) + 1) >> 1 << 1;
    return qn;
}

255
/* Convert the quantized vector to an index */
256
static inline uint32_t celt_icwrsi(uint32_t N, uint32_t K, const int *y)
257 258 259 260 261 262 263 264 265 266
{
    int i, idx = 0, sum = 0;
    for (i = N - 1; i >= 0; i--) {
        const uint32_t i_s = CELT_PVQ_U(N - i, sum + FFABS(y[i]) + 1);
        idx += CELT_PVQ_U(N - i, sum) + (y[i] < 0)*i_s;
        sum += FFABS(y[i]);
    }
    return idx;
}

267 268 269 270
// this code was adapted from libopus
static inline uint64_t celt_cwrsi(uint32_t N, uint32_t K, uint32_t i, int *y)
{
    uint64_t norm = 0;
271
    uint32_t q, p;
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350
    int s, val;
    int k0;

    while (N > 2) {
        /*Lots of pulses case:*/
        if (K >= N) {
            const uint32_t *row = ff_celt_pvq_u_row[N];

            /* Are the pulses in this dimension negative? */
            p  = row[K + 1];
            s  = -(i >= p);
            i -= p & s;

            /*Count how many pulses were placed in this dimension.*/
            k0 = K;
            q = row[N];
            if (q > i) {
                K = N;
                do {
                    p = ff_celt_pvq_u_row[--K][N];
                } while (p > i);
            } else
                for (p = row[K]; p > i; p = row[K])
                    K--;

            i    -= p;
            val   = (k0 - K + s) ^ s;
            norm += val * val;
            *y++  = val;
        } else { /*Lots of dimensions case:*/
            /*Are there any pulses in this dimension at all?*/
            p = ff_celt_pvq_u_row[K    ][N];
            q = ff_celt_pvq_u_row[K + 1][N];

            if (p <= i && i < q) {
                i -= p;
                *y++ = 0;
            } else {
                /*Are the pulses in this dimension negative?*/
                s  = -(i >= q);
                i -= q & s;

                /*Count how many pulses were placed in this dimension.*/
                k0 = K;
                do p = ff_celt_pvq_u_row[--K][N];
                while (p > i);

                i    -= p;
                val   = (k0 - K + s) ^ s;
                norm += val * val;
                *y++  = val;
            }
        }
        N--;
    }

    /* N == 2 */
    p  = 2 * K + 1;
    s  = -(i >= p);
    i -= p & s;
    k0 = K;
    K  = (i + 1) / 2;

    if (K)
        i -= 2 * K - 1;

    val   = (k0 - K + s) ^ s;
    norm += val * val;
    *y++  = val;

    /* N==1 */
    s     = -i;
    val   = (K + s) ^ s;
    norm += val * val;
    *y    = val;

    return norm;
}

351 352
static inline void celt_encode_pulses(OpusRangeCoder *rc, int *y, uint32_t N, uint32_t K)
{
353
    ff_opus_rc_enc_uint(rc, celt_icwrsi(N, K, y), CELT_PVQ_V(N, K));
354 355
}

356 357 358 359 360 361
static inline float celt_decode_pulses(OpusRangeCoder *rc, int *y, uint32_t N, uint32_t K)
{
    const uint32_t idx = ff_opus_rc_dec_uint(rc, CELT_PVQ_V(N, K));
    return celt_cwrsi(N, K, idx, y);
}

362 363 364 365
/*
 * Faster than libopus's search, operates entirely in the signed domain.
 * Slightly worse/better depending on N, K and the input vector.
 */
366
static float ppp_pvq_search_c(float *X, int *y, int K, int N)
367
{
368 369
    int i, y_norm = 0;
    float res = 0.0f, xy_norm = 0.0f;
370 371 372 373

    for (i = 0; i < N; i++)
        res += FFABS(X[i]);

374
    res = K/(res + FLT_EPSILON);
375 376 377 378 379 380 381 382 383

    for (i = 0; i < N; i++) {
        y[i] = lrintf(res*X[i]);
        y_norm  += y[i]*y[i];
        xy_norm += y[i]*X[i];
        K -= FFABS(y[i]);
    }

    while (K) {
384
        int max_idx = 0, phase = FFSIGN(K);
385
        float max_num = 0.0f;
386
        float max_den = 1.0f;
387 388 389
        y_norm += 1.0f;

        for (i = 0; i < N; i++) {
390 391 392 393
            /* If the sum has been overshot and the best place has 0 pulses allocated
             * to it, attempting to decrease it further will actually increase the
             * sum. Prevent this by disregarding any 0 positions when decrementing. */
            const int ca = 1 ^ ((y[i] == 0) & (phase < 0));
394
            const int y_new = y_norm  + 2*phase*FFABS(y[i]);
395 396
            float xy_new = xy_norm + 1*phase*FFABS(X[i]);
            xy_new = xy_new * xy_new;
397
            if (ca && (max_den*xy_new) > (y_new*max_num)) {
398 399 400 401 402 403 404 405 406 407 408 409 410
                max_den = y_new;
                max_num = xy_new;
                max_idx = i;
            }
        }

        K -= phase;

        phase *= FFSIGN(X[max_idx]);
        xy_norm += 1*phase*X[max_idx];
        y_norm  += 2*phase*y[max_idx];
        y[max_idx] += phase;
    }
411

412
    return (float)y_norm;
413 414 415
}

static uint32_t celt_alg_quant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K,
416
                               enum CeltSpread spread, uint32_t blocks, float gain,
417
                               CeltPVQ *pvq)
418
{
419
    int *y = pvq->qcoeff;
420 421

    celt_exp_rotation(X, N, blocks, K, spread, 1);
422
    gain /= sqrtf(pvq->pvq_search(X, y, K, N));
423
    celt_encode_pulses(rc, y,  N, K);
424 425
    celt_normalize_residual(y, X, N, gain);
    celt_exp_rotation(X, N, blocks, K, spread, 0);
426 427 428
    return celt_extract_collapse_mask(y, N, blocks);
}

429 430 431
/** Decode pulse vector and combine the result with the pitch vector to produce
    the final normalised signal in the current band. */
static uint32_t celt_alg_unquant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K,
432
                                 enum CeltSpread spread, uint32_t blocks, float gain,
433
                                 CeltPVQ *pvq)
434
{
435
    int *y = pvq->qcoeff;
436 437 438

    gain /= sqrtf(celt_decode_pulses(rc, y, N, K));
    celt_normalize_residual(y, X, N, gain);
439
    celt_exp_rotation(X, N, blocks, K, spread, 0);
440 441 442
    return celt_extract_collapse_mask(y, N, blocks);
}

443 444
static int celt_calc_theta(const float *X, const float *Y, int coupling, int N)
{
445
    int i;
446
    float e[2] = { 0.0f, 0.0f };
447 448 449 450 451 452 453 454 455
    if (coupling) { /* Coupling case */
        for (i = 0; i < N; i++) {
            e[0] += (X[i] + Y[i])*(X[i] + Y[i]);
            e[1] += (X[i] - Y[i])*(X[i] - Y[i]);
        }
    } else {
        for (i = 0; i < N; i++) {
            e[0] += X[i]*X[i];
            e[1] += Y[i]*Y[i];
456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475
        }
    }
    return lrintf(32768.0f*atan2f(sqrtf(e[1]), sqrtf(e[0]))/M_PI);
}

static void celt_stereo_is_decouple(float *X, float *Y, float e_l, float e_r, int N)
{
    int i;
    const float energy_n = 1.0f/(sqrtf(e_l*e_l + e_r*e_r) + FLT_EPSILON);
    e_l *= energy_n;
    e_r *= energy_n;
    for (i = 0; i < N; i++)
        X[i] = e_l*X[i] + e_r*Y[i];
}

static void celt_stereo_ms_decouple(float *X, float *Y, int N)
{
    int i;
    for (i = 0; i < N; i++) {
        const float Xret = X[i];
476 477
        X[i] = (X[i] + Y[i])*M_SQRT1_2;
        Y[i] = (Y[i] - Xret)*M_SQRT1_2;
478 479 480
    }
}

481 482 483 484 485 486 487 488
static av_always_inline uint32_t quant_band_template(CeltPVQ *pvq, CeltFrame *f,
                                                     OpusRangeCoder *rc,
                                                     const int band, float *X,
                                                     float *Y, int N, int b,
                                                     uint32_t blocks, float *lowband,
                                                     int duration, float *lowband_out,
                                                     int level, float gain,
                                                     float *lowband_scratch,
489
                                                     int fill, int quant)
490
{
491
    int i;
492
    const uint8_t *cache;
493
    int stereo = !!Y, split = stereo;
494
    int imid = 0, iside = 0;
495
    uint32_t N0 = N;
496
    int N_B = N / blocks;
497
    int N_B0 = N_B;
498 499 500 501 502 503 504 505 506 507
    int B0 = blocks;
    int time_divide = 0;
    int recombine = 0;
    int inv = 0;
    float mid = 0, side = 0;
    int longblocks = (B0 == 1);
    uint32_t cm = 0;

    if (N == 1) {
        float *x = X;
508
        for (i = 0; i <= stereo; i++) {
509 510 511 512 513 514 515 516
            int sign = 0;
            if (f->remaining2 >= 1 << 3) {
                if (quant) {
                    sign = x[0] < 0;
                    ff_opus_rc_put_raw(rc, sign, 1);
                } else {
                    sign = ff_opus_rc_get_raw(rc, 1);
                }
517 518
                f->remaining2 -= 1 << 3;
            }
519
            x[0] = 1.0f - 2.0f*sign;
520 521 522 523 524 525 526
            x = Y;
        }
        if (lowband_out)
            lowband_out[0] = X[0];
        return 1;
    }

527
    if (!stereo && level == 0) {
528 529 530 531 532 533 534 535
        int tf_change = f->tf_change[band];
        int k;
        if (tf_change > 0)
            recombine = tf_change;
        /* Band recombining to increase frequency resolution */

        if (lowband &&
            (recombine || ((N_B & 1) == 0 && tf_change < 0) || B0 > 1)) {
536 537
            for (i = 0; i < N; i++)
                lowband_scratch[i] = lowband[i];
538 539 540 541
            lowband = lowband_scratch;
        }

        for (k = 0; k < recombine; k++) {
542 543
            if (quant || lowband)
                celt_haar1(quant ? X : lowband, N >> k, 1 << k);
544 545 546 547 548 549 550
            fill = ff_celt_bit_interleave[fill & 0xF] | ff_celt_bit_interleave[fill >> 4] << 2;
        }
        blocks >>= recombine;
        N_B <<= recombine;

        /* Increasing the time resolution */
        while ((N_B & 1) == 0 && tf_change < 0) {
551 552
            if (quant || lowband)
                celt_haar1(quant ? X : lowband, N_B, blocks);
553 554 555 556 557 558 559
            fill |= fill << blocks;
            blocks <<= 1;
            N_B >>= 1;
            time_divide++;
            tf_change++;
        }
        B0 = blocks;
560
        N_B0 = N_B;
561 562

        /* Reorganize the samples in time order instead of frequency order */
563
        if (B0 > 1 && (quant || lowband))
564
            celt_deinterleave_hadamard(pvq->hadamard_tmp, quant ? X : lowband,
565 566
                                       N_B >> recombine, B0 << recombine,
                                       longblocks);
567 568 569 570 571
    }

    /* If we need 1.5 more bit than we can produce, split the band in two. */
    cache = ff_celt_cache_bits +
            ff_celt_cache_index[(duration + 1) * CELT_MAX_BANDS + band];
572
    if (!stereo && duration >= 0 && b > cache[cache[0]] + 12 && N > 2) {
573 574 575 576 577 578 579 580 581 582 583
        N >>= 1;
        Y = X + N;
        split = 1;
        duration -= 1;
        if (blocks == 1)
            fill = (fill & 1) | (fill << 1);
        blocks = (blocks + 1) >> 1;
    }

    if (split) {
        int qn;
584
        int itheta = quant ? celt_calc_theta(X, Y, stereo, N) : 0;
585 586 587 588 589 590 591 592 593
        int mbits, sbits, delta;
        int qalloc;
        int pulse_cap;
        int offset;
        int orig_fill;
        int tell;

        /* Decide on the resolution to give to the split parameter theta */
        pulse_cap = ff_celt_log_freq_range[band] + duration * 8;
594
        offset = (pulse_cap >> 1) - (stereo && N == 2 ? CELT_QTHETA_OFFSET_TWOPHASE :
595
                                                          CELT_QTHETA_OFFSET);
596 597
        qn = (stereo && band >= f->intensity_stereo) ? 1 :
             celt_compute_qn(N, b, offset, pulse_cap, stereo);
598 599
        tell = opus_rc_tell_frac(rc);
        if (qn != 1) {
600 601
            if (quant)
                itheta = (itheta*qn + 8192) >> 14;
602 603
            /* Entropy coding of the angle. We use a uniform pdf for the
             * time split, a step for stereo, and a triangular one for the rest. */
604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623
            if (quant) {
                if (stereo && N > 2)
                    ff_opus_rc_enc_uint_step(rc, itheta, qn / 2);
                else if (stereo || B0 > 1)
                    ff_opus_rc_enc_uint(rc, itheta, qn + 1);
                else
                    ff_opus_rc_enc_uint_tri(rc, itheta, qn);
                itheta = itheta * 16384 / qn;
                if (stereo) {
                    if (itheta == 0)
                        celt_stereo_is_decouple(X, Y, f->block[0].lin_energy[band],
                                                f->block[1].lin_energy[band], N);
                    else
                        celt_stereo_ms_decouple(X, Y, N);
                }
            } else {
                if (stereo && N > 2)
                    itheta = ff_opus_rc_dec_uint_step(rc, qn / 2);
                else if (stereo || B0 > 1)
                    itheta = ff_opus_rc_dec_uint(rc, qn+1);
624
                else
625 626
                    itheta = ff_opus_rc_dec_uint_tri(rc, qn);
                itheta = itheta * 16384 / qn;
627
            }
628
        } else if (stereo) {
629 630 631 632 633 634 635 636 637 638 639 640 641 642
            if (quant) {
                inv = itheta > 8192;
                 if (inv) {
                    for (i = 0; i < N; i++)
                       Y[i] *= -1;
                 }
                 celt_stereo_is_decouple(X, Y, f->block[0].lin_energy[band],
                                         f->block[1].lin_energy[band], N);

                if (b > 2 << 3 && f->remaining2 > 2 << 3) {
                    ff_opus_rc_enc_log(rc, inv, 2);
                } else {
                    inv = 0;
                }
643
            } else {
644
                inv = (b > 2 << 3 && f->remaining2 > 2 << 3) ? ff_opus_rc_dec_log(rc, 2) : 0;
645
                inv = f->apply_phase_inv ? inv : 0;
646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676
            }
            itheta = 0;
        }
        qalloc = opus_rc_tell_frac(rc) - tell;
        b -= qalloc;

        orig_fill = fill;
        if (itheta == 0) {
            imid = 32767;
            iside = 0;
            fill = av_mod_uintp2(fill, blocks);
            delta = -16384;
        } else if (itheta == 16384) {
            imid = 0;
            iside = 32767;
            fill &= ((1 << blocks) - 1) << blocks;
            delta = 16384;
        } else {
            imid = celt_cos(itheta);
            iside = celt_cos(16384-itheta);
            /* This is the mid vs side allocation that minimizes squared error
            in that band. */
            delta = ROUND_MUL16((N - 1) << 7, celt_log2tan(iside, imid));
        }

        mid  = imid  / 32768.0f;
        side = iside / 32768.0f;

        /* This is a special case for N=2 that only works for stereo and takes
        advantage of the fact that mid and side are orthogonal to encode
        the side with just one bit. */
677
        if (N == 2 && stereo) {
678 679 680 681 682 683 684 685 686 687 688 689 690 691
            int c;
            int sign = 0;
            float tmp;
            float *x2, *y2;
            mbits = b;
            /* Only need one bit for the side */
            sbits = (itheta != 0 && itheta != 16384) ? 1 << 3 : 0;
            mbits -= sbits;
            c = (itheta > 8192);
            f->remaining2 -= qalloc+sbits;

            x2 = c ? Y : X;
            y2 = c ? X : Y;
            if (sbits) {
692 693 694 695 696 697
                if (quant) {
                    sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
                    ff_opus_rc_put_raw(rc, sign, 1);
                } else {
                    sign = ff_opus_rc_get_raw(rc, 1);
                }
698 699 700 701
            }
            sign = 1 - 2 * sign;
            /* We use orig_fill here because we want to fold the side, but if
            itheta==16384, we'll have cleared the low bits of fill. */
702 703
            cm = pvq->quant_band(pvq, f, rc, band, x2, NULL, N, mbits, blocks, lowband, duration,
                                 lowband_out, level, gain, lowband_scratch, orig_fill);
704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723
            /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
            and there's no need to worry about mixing with the other channel. */
            y2[0] = -sign * x2[1];
            y2[1] =  sign * x2[0];
            X[0] *= mid;
            X[1] *= mid;
            Y[0] *= side;
            Y[1] *= side;
            tmp = X[0];
            X[0] = tmp - Y[0];
            Y[0] = tmp + Y[0];
            tmp = X[1];
            X[1] = tmp - Y[1];
            Y[1] = tmp + Y[1];
        } else {
            /* "Normal" split code */
            float *next_lowband2     = NULL;
            float *next_lowband_out1 = NULL;
            int next_level = 0;
            int rebalance;
724
            uint32_t cmt;
725 726 727

            /* Give more bits to low-energy MDCTs than they would
             * otherwise deserve */
728
            if (B0 > 1 && !stereo && (itheta & 0x3fff)) {
729 730 731 732 733 734 735 736 737 738 739 740
                if (itheta > 8192)
                    /* Rough approximation for pre-echo masking */
                    delta -= delta >> (4 - duration);
                else
                    /* Corresponds to a forward-masking slope of
                     * 1.5 dB per 10 ms */
                    delta = FFMIN(0, delta + (N << 3 >> (5 - duration)));
            }
            mbits = av_clip((b - delta) / 2, 0, b);
            sbits = b - mbits;
            f->remaining2 -= qalloc;

741
            if (lowband && !stereo)
742 743 744 745
                next_lowband2 = lowband + N; /* >32-bit split case */

            /* Only stereo needs to pass on lowband_out.
             * Otherwise, it's handled at the end */
746
            if (stereo)
747 748 749 750 751 752 753 754
                next_lowband_out1 = lowband_out;
            else
                next_level = level + 1;

            rebalance = f->remaining2;
            if (mbits >= sbits) {
                /* In stereo mode, we do not apply a scaling to the mid
                 * because we need the normalized mid for folding later */
755 756 757
                cm = pvq->quant_band(pvq, f, rc, band, X, NULL, N, mbits, blocks,
                                     lowband, duration, next_lowband_out1, next_level,
                                     stereo ? 1.0f : (gain * mid), lowband_scratch, fill);
758 759 760 761 762 763
                rebalance = mbits - (rebalance - f->remaining2);
                if (rebalance > 3 << 3 && itheta != 0)
                    sbits += rebalance - (3 << 3);

                /* For a stereo split, the high bits of fill are always zero,
                 * so no folding will be done to the side. */
764 765 766
                cmt = pvq->quant_band(pvq, f, rc, band, Y, NULL, N, sbits, blocks,
                                      next_lowband2, duration, NULL, next_level,
                                      gain * side, NULL, fill >> blocks);
767
                cm |= cmt << ((B0 >> 1) & (stereo - 1));
768 769 770
            } else {
                /* For a stereo split, the high bits of fill are always zero,
                 * so no folding will be done to the side. */
771 772 773
                cm = pvq->quant_band(pvq, f, rc, band, Y, NULL, N, sbits, blocks,
                                     next_lowband2, duration, NULL, next_level,
                                     gain * side, NULL, fill >> blocks);
774
                cm <<= ((B0 >> 1) & (stereo - 1));
775 776 777 778 779 780
                rebalance = sbits - (rebalance - f->remaining2);
                if (rebalance > 3 << 3 && itheta != 16384)
                    mbits += rebalance - (3 << 3);

                /* In stereo mode, we do not apply a scaling to the mid because
                 * we need the normalized mid for folding later */
781 782 783
                cm |= pvq->quant_band(pvq, f, rc, band, X, NULL, N, mbits, blocks,
                                      lowband, duration, next_lowband_out1, next_level,
                                      stereo ? 1.0f : (gain * mid), lowband_scratch, fill);
784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799
            }
        }
    } else {
        /* This is the basic no-split case */
        uint32_t q         = celt_bits2pulses(cache, b);
        uint32_t curr_bits = celt_pulses2bits(cache, q);
        f->remaining2 -= curr_bits;

        /* Ensures we can never bust the budget */
        while (f->remaining2 < 0 && q > 0) {
            f->remaining2 += curr_bits;
            curr_bits      = celt_pulses2bits(cache, --q);
            f->remaining2 -= curr_bits;
        }

        if (q != 0) {
800 801 802
            /* Finally do the actual (de)quantization */
            if (quant) {
                cm = celt_alg_quant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
803
                                    f->spread, blocks, gain, pvq);
804 805
            } else {
                cm = celt_alg_unquant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
806
                                      f->spread, blocks, gain, pvq);
807
            }
808 809 810 811
        } else {
            /* If there's no pulse, fill the band anyway */
            uint32_t cm_mask = (1 << blocks) - 1;
            fill &= cm_mask;
812
            if (fill) {
813 814
                if (!lowband) {
                    /* Noise */
815 816
                    for (i = 0; i < N; i++)
                        X[i] = (((int32_t)celt_rng(f)) >> 20);
817 818 819
                    cm = cm_mask;
                } else {
                    /* Folded spectrum */
820
                    for (i = 0; i < N; i++) {
821
                        /* About 48 dB below the "normal" folding level */
822
                        X[i] = lowband[i] + (((celt_rng(f)) & 0x8000) ? 1.0f / 256 : -1.0f / 256);
823 824 825 826
                    }
                    cm = fill;
                }
                celt_renormalize_vector(X, N, gain);
827 828
            } else {
                memset(X, 0, N*sizeof(float));
829 830 831 832 833
            }
        }
    }

    /* This code is used by the decoder and by the resynthesis-enabled encoder */
834 835
    if (stereo) {
        if (N > 2)
836 837
            celt_stereo_merge(X, Y, mid, N);
        if (inv) {
838 839
            for (i = 0; i < N; i++)
                Y[i] *= -1;
840 841 842 843 844 845
        }
    } else if (level == 0) {
        int k;

        /* Undo the sample reorganization going from time order to frequency order */
        if (B0 > 1)
846
            celt_interleave_hadamard(pvq->hadamard_tmp, X, N_B >> recombine,
847
                                     B0 << recombine, longblocks);
848 849 850 851 852 853 854 855 856

        /* Undo time-freq changes that we did earlier */
        N_B = N_B0;
        blocks = B0;
        for (k = 0; k < time_divide; k++) {
            blocks >>= 1;
            N_B <<= 1;
            cm |= cm >> blocks;
            celt_haar1(X, N_B, blocks);
857
        }
858 859 860 861 862 863 864 865 866 867

        for (k = 0; k < recombine; k++) {
            cm = ff_celt_bit_deinterleave[cm];
            celt_haar1(X, N0>>k, 1<<k);
        }
        blocks <<= recombine;

        /* Scale output for later folding */
        if (lowband_out) {
            float n = sqrtf(N0);
868 869
            for (i = 0; i < N0; i++)
                lowband_out[i] = n * X[i];
870 871
        }
        cm = av_mod_uintp2(cm, blocks);
872 873
    }

874 875
    return cm;
}
876

877
static QUANT_FN(pvq_decode_band)
878
{
879
#if CONFIG_OPUS_DECODER
880
    return quant_band_template(pvq, f, rc, band, X, Y, N, b, blocks, lowband, duration,
881
                               lowband_out, level, gain, lowband_scratch, fill, 0);
882 883 884
#else
    return 0;
#endif
885 886
}

887
static QUANT_FN(pvq_encode_band)
888
{
889
#if CONFIG_OPUS_ENCODER
890
    return quant_band_template(pvq, f, rc, band, X, Y, N, b, blocks, lowband, duration,
891
                               lowband_out, level, gain, lowband_scratch, fill, 1);
892 893 894
#else
    return 0;
#endif
895
}
896

897
int av_cold ff_celt_pvq_init(CeltPVQ **pvq, int encode)
898 899 900 901 902
{
    CeltPVQ *s = av_malloc(sizeof(CeltPVQ));
    if (!s)
        return AVERROR(ENOMEM);

903 904
    s->pvq_search = ppp_pvq_search_c;
    s->quant_band = encode ? pvq_encode_band : pvq_decode_band;
905

906
    if (CONFIG_OPUS_ENCODER && ARCH_X86)
907
        ff_celt_pvq_init_x86(s);
908

909 910 911 912 913 914 915 916 917
    *pvq = s;

    return 0;
}

void av_cold ff_celt_pvq_uninit(CeltPVQ **pvq)
{
    av_freep(pvq);
}