2 * Copyright (c) 2007-2008 CSIRO
3 * Copyright (c) 2007-2009 Xiph.Org Foundation
4 * Copyright (c) 2008-2009 Gregory Maxwell
5 * Copyright (c) 2012 Andrew D'Addesio
6 * Copyright (c) 2013-2014 Mozilla Corporation
7 * Copyright (c) 2017 Rostislav Pehlivanov <atomnuker@gmail.com>
9 * This file is part of FFmpeg.
11 * FFmpeg is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public
13 * License as published by the Free Software Foundation; either
14 * version 2.1 of the License, or (at your option) any later version.
16 * FFmpeg is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with FFmpeg; if not, write to the Free Software
23 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
29 #define CELT_PVQ_U(n, k) (ff_celt_pvq_u_row[FFMIN(n, k)][FFMAX(n, k)])
30 #define CELT_PVQ_V(n, k) (CELT_PVQ_U(n, k) + CELT_PVQ_U(n, (k) + 1))
32 static inline int16_t celt_cos(int16_t x)
34 x = (MUL16(x, x) + 4096) >> 13;
35 x = (32767-x) + ROUND_MUL16(x, (-7651 + ROUND_MUL16(x, (8277 + ROUND_MUL16(-626, x)))));
39 static inline int celt_log2tan(int isin, int icos)
46 return (ls << 11) - (lc << 11) +
47 ROUND_MUL16(isin, ROUND_MUL16(isin, -2597) + 7932) -
48 ROUND_MUL16(icos, ROUND_MUL16(icos, -2597) + 7932);
51 static inline int celt_bits2pulses(const uint8_t *cache, int bits)
53 // TODO: Find the size of cache and make it into an array in the parameters list
59 for (i = 0; i < 6; i++) {
60 int center = (low + high + 1) >> 1;
61 if (cache[center] >= bits)
67 return (bits - (low == 0 ? -1 : cache[low]) <= cache[high] - bits) ? low : high;
70 static inline int celt_pulses2bits(const uint8_t *cache, int pulses)
72 // TODO: Find the size of cache and make it into an array in the parameters list
73 return (pulses == 0) ? 0 : cache[pulses] + 1;
76 static inline void celt_normalize_residual(const int * av_restrict iy, float * av_restrict X,
80 for (i = 0; i < N; i++)
84 static void celt_exp_rotation_impl(float *X, uint32_t len, uint32_t stride,
91 for (i = 0; i < len - stride; i++) {
95 Xptr[stride] = c * x2 + s * x1;
96 *Xptr++ = c * x1 - s * x2;
99 Xptr = &X[len - 2 * stride - 1];
100 for (i = len - 2 * stride - 1; i >= 0; i--) {
104 Xptr[stride] = c * x2 + s * x1;
105 *Xptr-- = c * x1 - s * x2;
109 static inline void celt_exp_rotation(float *X, uint32_t len,
110 uint32_t stride, uint32_t K,
111 enum CeltSpread spread, const int encode)
113 uint32_t stride2 = 0;
118 if (2*K >= len || spread == CELT_SPREAD_NONE)
121 gain = (float)len / (len + (20 - 5*spread) * K);
122 theta = M_PI * gain * gain / 4;
127 if (len >= stride << 3) {
129 /* This is just a simple (equivalent) way of computing sqrt(len/stride) with rounding.
130 It's basically incrementing long as (stride2+0.5)^2 < len/stride. */
131 while ((stride2 * stride2 + stride2) * stride + (stride >> 2) < len)
135 /*NOTE: As a minor optimization, we could be passing around log2(B), not B, for both this and for
136 extract_collapse_mask().*/
138 for (i = 0; i < stride; i++) {
140 celt_exp_rotation_impl(X + i * len, len, 1, c, -s);
142 celt_exp_rotation_impl(X + i * len, len, stride2, s, -c);
145 celt_exp_rotation_impl(X + i * len, len, stride2, s, c);
146 celt_exp_rotation_impl(X + i * len, len, 1, c, s);
151 static inline uint32_t celt_extract_collapse_mask(const int *iy, uint32_t N, uint32_t B)
153 uint32_t collapse_mask;
160 /*NOTE: As a minor optimization, we could be passing around log2(B), not B, for both this and for
164 for (i = 0; i < B; i++)
165 for (j = 0; j < N0; j++)
166 collapse_mask |= (iy[i*N0+j]!=0)<<i;
167 return collapse_mask;
170 static inline void celt_stereo_merge(float *X, float *Y, float mid, int N)
173 float xp = 0, side = 0;
178 /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
179 for (i = 0; i < N; i++) {
184 /* Compensating for the mid normalization */
187 E[0] = mid2 * mid2 + side - 2 * xp;
188 E[1] = mid2 * mid2 + side + 2 * xp;
189 if (E[0] < 6e-4f || E[1] < 6e-4f) {
190 for (i = 0; i < N; i++)
196 gain[0] = 1.0f / sqrtf(t);
198 gain[1] = 1.0f / sqrtf(t);
200 for (i = 0; i < N; i++) {
202 /* Apply mid scaling (side is already scaled) */
203 value[0] = mid * X[i];
205 X[i] = gain[0] * (value[0] - value[1]);
206 Y[i] = gain[1] * (value[0] + value[1]);
210 static void celt_interleave_hadamard(float *tmp, float *X, int N0,
211 int stride, int hadamard)
217 const uint8_t *ordery = ff_celt_hadamard_ordery + stride - 2;
218 for (i = 0; i < stride; i++)
219 for (j = 0; j < N0; j++)
220 tmp[j*stride+i] = X[ordery[i]*N0+j];
222 for (i = 0; i < stride; i++)
223 for (j = 0; j < N0; j++)
224 tmp[j*stride+i] = X[i*N0+j];
227 for (i = 0; i < N; i++)
231 static void celt_deinterleave_hadamard(float *tmp, float *X, int N0,
232 int stride, int hadamard)
238 const uint8_t *ordery = ff_celt_hadamard_ordery + stride - 2;
239 for (i = 0; i < stride; i++)
240 for (j = 0; j < N0; j++)
241 tmp[ordery[i]*N0+j] = X[j*stride+i];
243 for (i = 0; i < stride; i++)
244 for (j = 0; j < N0; j++)
245 tmp[i*N0+j] = X[j*stride+i];
248 for (i = 0; i < N; i++)
252 static void celt_haar1(float *X, int N0, int stride)
256 for (i = 0; i < stride; i++) {
257 for (j = 0; j < N0; j++) {
258 float x0 = X[stride * (2 * j + 0) + i];
259 float x1 = X[stride * (2 * j + 1) + i];
260 X[stride * (2 * j + 0) + i] = (x0 + x1) * M_SQRT1_2;
261 X[stride * (2 * j + 1) + i] = (x0 - x1) * M_SQRT1_2;
266 static inline int celt_compute_qn(int N, int b, int offset, int pulse_cap,
271 if (dualstereo && N == 2)
274 /* The upper limit ensures that in a stereo split with itheta==16384, we'll
275 * always have enough bits left over to code at least one pulse in the
276 * side; otherwise it would collapse, since it doesn't get folded. */
277 qb = FFMIN3(b - pulse_cap - (4 << 3), (b + N2 * offset) / N2, 8 << 3);
278 qn = (qb < (1 << 3 >> 1)) ? 1 : ((ff_celt_qn_exp2[qb & 0x7] >> (14 - (qb >> 3))) + 1) >> 1 << 1;
282 /* Convert the quantized vector to an index */
283 static inline uint32_t celt_icwrsi(uint32_t N, uint32_t K, const int *y)
285 int i, idx = 0, sum = 0;
286 for (i = N - 1; i >= 0; i--) {
287 const uint32_t i_s = CELT_PVQ_U(N - i, sum + FFABS(y[i]) + 1);
288 idx += CELT_PVQ_U(N - i, sum) + (y[i] < 0)*i_s;
294 // this code was adapted from libopus
295 static inline uint64_t celt_cwrsi(uint32_t N, uint32_t K, uint32_t i, int *y)
305 /*Lots of pulses case:*/
307 const uint32_t *row = ff_celt_pvq_u_row[N];
309 /* Are the pulses in this dimension negative? */
314 /*Count how many pulses were placed in this dimension.*/
320 p = ff_celt_pvq_u_row[--K][N];
323 for (p = row[K]; p > i; p = row[K])
327 val = (k0 - K + s) ^ s;
330 } else { /*Lots of dimensions case:*/
331 /*Are there any pulses in this dimension at all?*/
332 p = ff_celt_pvq_u_row[K ][N];
333 q = ff_celt_pvq_u_row[K + 1][N];
335 if (p <= i && i < q) {
339 /*Are the pulses in this dimension negative?*/
343 /*Count how many pulses were placed in this dimension.*/
345 do p = ff_celt_pvq_u_row[--K][N];
349 val = (k0 - K + s) ^ s;
367 val = (k0 - K + s) ^ s;
380 static inline void celt_encode_pulses(OpusRangeCoder *rc, int *y, uint32_t N, uint32_t K)
382 ff_opus_rc_enc_uint(rc, celt_icwrsi(N, K, y), CELT_PVQ_V(N, K));
385 static inline float celt_decode_pulses(OpusRangeCoder *rc, int *y, uint32_t N, uint32_t K)
387 const uint32_t idx = ff_opus_rc_dec_uint(rc, CELT_PVQ_V(N, K));
388 return celt_cwrsi(N, K, idx, y);
392 * Faster than libopus's search, operates entirely in the signed domain.
393 * Slightly worse/better depending on N, K and the input vector.
395 static int celt_pvq_search(float *X, int *y, int K, int N)
398 float res = 0.0f, xy_norm = 0.0f;
400 for (i = 0; i < N; i++)
403 res = K/(res + FLT_EPSILON);
405 for (i = 0; i < N; i++) {
406 y[i] = lrintf(res*X[i]);
408 xy_norm += y[i]*X[i];
413 int max_idx = 0, max_den = 1, phase = FFSIGN(K);
414 float max_num = 0.0f;
417 for (i = 0; i < N; i++) {
418 /* If the sum has been overshot and the best place has 0 pulses allocated
419 * to it, attempting to decrease it further will actually increase the
420 * sum. Prevent this by disregarding any 0 positions when decrementing. */
421 const int ca = 1 ^ ((y[i] == 0) & (phase < 0));
422 const int y_new = y_norm + 2*phase*FFABS(y[i]);
423 float xy_new = xy_norm + 1*phase*FFABS(X[i]);
424 xy_new = xy_new * xy_new;
425 if (ca && (max_den*xy_new) > (y_new*max_num)) {
434 phase *= FFSIGN(X[max_idx]);
435 xy_norm += 1*phase*X[max_idx];
436 y_norm += 2*phase*y[max_idx];
443 static uint32_t celt_alg_quant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K,
444 enum CeltSpread spread, uint32_t blocks, float gain)
448 celt_exp_rotation(X, N, blocks, K, spread, 1);
449 gain /= sqrtf(celt_pvq_search(X, y, K, N));
450 celt_encode_pulses(rc, y, N, K);
451 celt_normalize_residual(y, X, N, gain);
452 celt_exp_rotation(X, N, blocks, K, spread, 0);
453 return celt_extract_collapse_mask(y, N, blocks);
456 /** Decode pulse vector and combine the result with the pitch vector to produce
457 the final normalised signal in the current band. */
458 static uint32_t celt_alg_unquant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K,
459 enum CeltSpread spread, uint32_t blocks, float gain)
463 gain /= sqrtf(celt_decode_pulses(rc, y, N, K));
464 celt_normalize_residual(y, X, N, gain);
465 celt_exp_rotation(X, N, blocks, K, spread, 0);
466 return celt_extract_collapse_mask(y, N, blocks);
469 uint32_t ff_celt_decode_band(CeltFrame *f, OpusRangeCoder *rc, const int band,
470 float *X, float *Y, int N, int b, uint32_t blocks,
471 float *lowband, int duration, float *lowband_out, int level,
472 float gain, float *lowband_scratch, int fill)
474 const uint8_t *cache;
475 int dualstereo, split;
476 int imid = 0, iside = 0;
484 float mid = 0, side = 0;
485 int longblocks = (B0 == 1);
488 N_B0 = N_B = N / blocks;
489 split = dualstereo = (Y != NULL);
492 /* special case for one sample */
495 for (i = 0; i <= dualstereo; i++) {
497 if (f->remaining2 >= 1<<3) {
498 sign = ff_opus_rc_get_raw(rc, 1);
499 f->remaining2 -= 1 << 3;
502 x[0] = sign ? -1.0f : 1.0f;
506 lowband_out[0] = X[0];
510 if (!dualstereo && level == 0) {
511 int tf_change = f->tf_change[band];
514 recombine = tf_change;
515 /* Band recombining to increase frequency resolution */
518 (recombine || ((N_B & 1) == 0 && tf_change < 0) || B0 > 1)) {
520 for (j = 0; j < N; j++)
521 lowband_scratch[j] = lowband[j];
522 lowband = lowband_scratch;
525 for (k = 0; k < recombine; k++) {
527 celt_haar1(lowband, N >> k, 1 << k);
528 fill = ff_celt_bit_interleave[fill & 0xF] | ff_celt_bit_interleave[fill >> 4] << 2;
530 blocks >>= recombine;
533 /* Increasing the time resolution */
534 while ((N_B & 1) == 0 && tf_change < 0) {
536 celt_haar1(lowband, N_B, blocks);
537 fill |= fill << blocks;
546 /* Reorganize the samples in time order instead of frequency order */
547 if (B0 > 1 && lowband)
548 celt_deinterleave_hadamard(f->scratch, lowband, N_B >> recombine,
549 B0 << recombine, longblocks);
552 /* If we need 1.5 more bit than we can produce, split the band in two. */
553 cache = ff_celt_cache_bits +
554 ff_celt_cache_index[(duration + 1) * CELT_MAX_BANDS + band];
555 if (!dualstereo && duration >= 0 && b > cache[cache[0]] + 12 && N > 2) {
561 fill = (fill & 1) | (fill << 1);
562 blocks = (blocks + 1) >> 1;
568 int mbits, sbits, delta;
575 /* Decide on the resolution to give to the split parameter theta */
576 pulse_cap = ff_celt_log_freq_range[band] + duration * 8;
577 offset = (pulse_cap >> 1) - (dualstereo && N == 2 ? CELT_QTHETA_OFFSET_TWOPHASE :
579 qn = (dualstereo && band >= f->intensity_stereo) ? 1 :
580 celt_compute_qn(N, b, offset, pulse_cap, dualstereo);
581 tell = opus_rc_tell_frac(rc);
583 /* Entropy coding of the angle. We use a uniform pdf for the
584 time split, a step for stereo, and a triangular one for the rest. */
585 if (dualstereo && N > 2)
586 itheta = ff_opus_rc_dec_uint_step(rc, qn/2);
587 else if (dualstereo || B0 > 1)
588 itheta = ff_opus_rc_dec_uint(rc, qn+1);
590 itheta = ff_opus_rc_dec_uint_tri(rc, qn);
591 itheta = itheta * 16384 / qn;
592 /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.
593 Let's do that at higher complexity */
594 } else if (dualstereo) {
595 inv = (b > 2 << 3 && f->remaining2 > 2 << 3) ? ff_opus_rc_dec_log(rc, 2) : 0;
598 qalloc = opus_rc_tell_frac(rc) - tell;
605 fill = av_mod_uintp2(fill, blocks);
607 } else if (itheta == 16384) {
610 fill &= ((1 << blocks) - 1) << blocks;
613 imid = celt_cos(itheta);
614 iside = celt_cos(16384-itheta);
615 /* This is the mid vs side allocation that minimizes squared error
617 delta = ROUND_MUL16((N - 1) << 7, celt_log2tan(iside, imid));
620 mid = imid / 32768.0f;
621 side = iside / 32768.0f;
623 /* This is a special case for N=2 that only works for stereo and takes
624 advantage of the fact that mid and side are orthogonal to encode
625 the side with just one bit. */
626 if (N == 2 && dualstereo) {
632 /* Only need one bit for the side */
633 sbits = (itheta != 0 && itheta != 16384) ? 1 << 3 : 0;
636 f->remaining2 -= qalloc+sbits;
641 sign = ff_opus_rc_get_raw(rc, 1);
643 /* We use orig_fill here because we want to fold the side, but if
644 itheta==16384, we'll have cleared the low bits of fill. */
645 cm = ff_celt_decode_band(f, rc, band, x2, NULL, N, mbits, blocks,
646 lowband, duration, lowband_out, level, gain,
647 lowband_scratch, orig_fill);
648 /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
649 and there's no need to worry about mixing with the other channel. */
650 y2[0] = -sign * x2[1];
651 y2[1] = sign * x2[0];
663 /* "Normal" split code */
664 float *next_lowband2 = NULL;
665 float *next_lowband_out1 = NULL;
669 /* Give more bits to low-energy MDCTs than they would
670 * otherwise deserve */
671 if (B0 > 1 && !dualstereo && (itheta & 0x3fff)) {
673 /* Rough approximation for pre-echo masking */
674 delta -= delta >> (4 - duration);
676 /* Corresponds to a forward-masking slope of
677 * 1.5 dB per 10 ms */
678 delta = FFMIN(0, delta + (N << 3 >> (5 - duration)));
680 mbits = av_clip((b - delta) / 2, 0, b);
682 f->remaining2 -= qalloc;
684 if (lowband && !dualstereo)
685 next_lowband2 = lowband + N; /* >32-bit split case */
687 /* Only stereo needs to pass on lowband_out.
688 * Otherwise, it's handled at the end */
690 next_lowband_out1 = lowband_out;
692 next_level = level + 1;
694 rebalance = f->remaining2;
695 if (mbits >= sbits) {
696 /* In stereo mode, we do not apply a scaling to the mid
697 * because we need the normalized mid for folding later */
698 cm = ff_celt_decode_band(f, rc, band, X, NULL, N, mbits, blocks,
699 lowband, duration, next_lowband_out1,
700 next_level, dualstereo ? 1.0f : (gain * mid),
701 lowband_scratch, fill);
703 rebalance = mbits - (rebalance - f->remaining2);
704 if (rebalance > 3 << 3 && itheta != 0)
705 sbits += rebalance - (3 << 3);
707 /* For a stereo split, the high bits of fill are always zero,
708 * so no folding will be done to the side. */
709 cm |= ff_celt_decode_band(f, rc, band, Y, NULL, N, sbits, blocks,
710 next_lowband2, duration, NULL,
711 next_level, gain * side, NULL,
712 fill >> blocks) << ((B0 >> 1) & (dualstereo - 1));
714 /* For a stereo split, the high bits of fill are always zero,
715 * so no folding will be done to the side. */
716 cm = ff_celt_decode_band(f, rc, band, Y, NULL, N, sbits, blocks,
717 next_lowband2, duration, NULL,
718 next_level, gain * side, NULL,
719 fill >> blocks) << ((B0 >> 1) & (dualstereo - 1));
721 rebalance = sbits - (rebalance - f->remaining2);
722 if (rebalance > 3 << 3 && itheta != 16384)
723 mbits += rebalance - (3 << 3);
725 /* In stereo mode, we do not apply a scaling to the mid because
726 * we need the normalized mid for folding later */
727 cm |= ff_celt_decode_band(f, rc, band, X, NULL, N, mbits, blocks,
728 lowband, duration, next_lowband_out1,
729 next_level, dualstereo ? 1.0f : (gain * mid),
730 lowband_scratch, fill);
734 /* This is the basic no-split case */
735 uint32_t q = celt_bits2pulses(cache, b);
736 uint32_t curr_bits = celt_pulses2bits(cache, q);
737 f->remaining2 -= curr_bits;
739 /* Ensures we can never bust the budget */
740 while (f->remaining2 < 0 && q > 0) {
741 f->remaining2 += curr_bits;
742 curr_bits = celt_pulses2bits(cache, --q);
743 f->remaining2 -= curr_bits;
747 /* Finally do the actual quantization */
748 cm = celt_alg_unquant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
749 f->spread, blocks, gain);
751 /* If there's no pulse, fill the band anyway */
753 uint32_t cm_mask = (1 << blocks) - 1;
756 for (j = 0; j < N; j++)
761 for (j = 0; j < N; j++)
762 X[j] = (((int32_t)celt_rng(f)) >> 20);
765 /* Folded spectrum */
766 for (j = 0; j < N; j++) {
767 /* About 48 dB below the "normal" folding level */
768 X[j] = lowband[j] + (((celt_rng(f)) & 0x8000) ? 1.0f / 256 : -1.0f / 256);
772 celt_renormalize_vector(X, N, gain);
777 /* This code is used by the decoder and by the resynthesis-enabled encoder */
781 celt_stereo_merge(X, Y, mid, N);
783 for (j = 0; j < N; j++)
786 } else if (level == 0) {
789 /* Undo the sample reorganization going from time order to frequency order */
791 celt_interleave_hadamard(f->scratch, X, N_B>>recombine,
792 B0<<recombine, longblocks);
794 /* Undo time-freq changes that we did earlier */
797 for (k = 0; k < time_divide; k++) {
801 celt_haar1(X, N_B, blocks);
804 for (k = 0; k < recombine; k++) {
805 cm = ff_celt_bit_deinterleave[cm];
806 celt_haar1(X, N0>>k, 1<<k);
808 blocks <<= recombine;
810 /* Scale output for later folding */
814 for (j = 0; j < N0; j++)
815 lowband_out[j] = n * X[j];
817 cm = av_mod_uintp2(cm, blocks);
823 /* This has to be, AND MUST BE done by the psychoacoustic system, this has a very
824 * big impact on the entire quantization and especially huge on transients */
825 static int celt_calc_theta(const float *X, const float *Y, int coupling, int N)
828 float e[2] = { 0.0f, 0.0f };
829 for (j = 0; j < N; j++) {
830 if (coupling) { /* Coupling case */
831 e[0] += (X[j] + Y[j])*(X[j] + Y[j]);
832 e[1] += (X[j] - Y[j])*(X[j] - Y[j]);
838 return lrintf(32768.0f*atan2f(sqrtf(e[1]), sqrtf(e[0]))/M_PI);
841 static void celt_stereo_is_decouple(float *X, float *Y, float e_l, float e_r, int N)
844 const float energy_n = 1.0f/(sqrtf(e_l*e_l + e_r*e_r) + FLT_EPSILON);
847 for (i = 0; i < N; i++)
848 X[i] = e_l*X[i] + e_r*Y[i];
851 static void celt_stereo_ms_decouple(float *X, float *Y, int N)
854 const float decouple_norm = 1.0f/sqrtf(1.0f + 1.0f);
855 for (i = 0; i < N; i++) {
856 const float Xret = X[i];
857 X[i] = (X[i] + Y[i])*decouple_norm;
858 Y[i] = (Y[i] - Xret)*decouple_norm;
862 uint32_t ff_celt_encode_band(CeltFrame *f, OpusRangeCoder *rc, const int band,
863 float *X, float *Y, int N, int b, uint32_t blocks,
864 float *lowband, int duration, float *lowband_out, int level,
865 float gain, float *lowband_scratch, int fill)
867 const uint8_t *cache;
868 int dualstereo, split;
869 int imid = 0, iside = 0;
871 int N_B = N / blocks;
877 float mid = 0, side = 0;
878 int longblocks = (B0 == 1);
881 split = dualstereo = (Y != NULL);
884 /* special case for one sample - the decoder's output will be +- 1.0f!!! */
887 for (i = 0; i <= dualstereo; i++) {
888 if (f->remaining2 >= 1<<3) {
889 ff_opus_rc_put_raw(rc, x[0] < 0, 1);
890 f->remaining2 -= 1 << 3;
893 x[0] = 1.0f - 2.0f*(x[0] < 0);
897 lowband_out[0] = X[0];
901 if (!dualstereo && level == 0) {
902 int tf_change = f->tf_change[band];
905 recombine = tf_change;
906 /* Band recombining to increase frequency resolution */
909 (recombine || ((N_B & 1) == 0 && tf_change < 0) || B0 > 1)) {
911 for (j = 0; j < N; j++)
912 lowband_scratch[j] = lowband[j];
913 lowband = lowband_scratch;
916 for (k = 0; k < recombine; k++) {
917 celt_haar1(X, N >> k, 1 << k);
918 fill = ff_celt_bit_interleave[fill & 0xF] | ff_celt_bit_interleave[fill >> 4] << 2;
920 blocks >>= recombine;
923 /* Increasing the time resolution */
924 while ((N_B & 1) == 0 && tf_change < 0) {
925 celt_haar1(X, N_B, blocks);
926 fill |= fill << blocks;
935 /* Reorganize the samples in time order instead of frequency order */
937 celt_deinterleave_hadamard(f->scratch, X, N_B >> recombine,
938 B0 << recombine, longblocks);
941 /* If we need 1.5 more bit than we can produce, split the band in two. */
942 cache = ff_celt_cache_bits +
943 ff_celt_cache_index[(duration + 1) * CELT_MAX_BANDS + band];
944 if (!dualstereo && duration >= 0 && b > cache[cache[0]] + 12 && N > 2) {
950 fill = (fill & 1) | (fill << 1);
951 blocks = (blocks + 1) >> 1;
956 int itheta = celt_calc_theta(X, Y, dualstereo, N);
957 int mbits, sbits, delta;
964 /* Decide on the resolution to give to the split parameter theta */
965 pulse_cap = ff_celt_log_freq_range[band] + duration * 8;
966 offset = (pulse_cap >> 1) - (dualstereo && N == 2 ? CELT_QTHETA_OFFSET_TWOPHASE :
968 qn = (dualstereo && band >= f->intensity_stereo) ? 1 :
969 celt_compute_qn(N, b, offset, pulse_cap, dualstereo);
970 tell = opus_rc_tell_frac(rc);
974 itheta = (itheta*qn + 8192) >> 14;
976 /* Entropy coding of the angle. We use a uniform pdf for the
977 * time split, a step for stereo, and a triangular one for the rest. */
978 if (dualstereo && N > 2)
979 ff_opus_rc_enc_uint_step(rc, itheta, qn / 2);
980 else if (dualstereo || B0 > 1)
981 ff_opus_rc_enc_uint(rc, itheta, qn + 1);
983 ff_opus_rc_enc_uint_tri(rc, itheta, qn);
984 itheta = itheta * 16384 / qn;
988 celt_stereo_is_decouple(X, Y, f->block[0].lin_energy[band],
989 f->block[1].lin_energy[band], N);
991 celt_stereo_ms_decouple(X, Y, N);
993 } else if (dualstereo) {
997 for (j = 0; j < N; j++)
1000 celt_stereo_is_decouple(X, Y, f->block[0].lin_energy[band],
1001 f->block[1].lin_energy[band], N);
1003 if (b > 2 << 3 && f->remaining2 > 2 << 3) {
1004 ff_opus_rc_enc_log(rc, inv, 2);
1011 qalloc = opus_rc_tell_frac(rc) - tell;
1018 fill = av_mod_uintp2(fill, blocks);
1020 } else if (itheta == 16384) {
1023 fill &= ((1 << blocks) - 1) << blocks;
1026 imid = celt_cos(itheta);
1027 iside = celt_cos(16384-itheta);
1028 /* This is the mid vs side allocation that minimizes squared error
1030 delta = ROUND_MUL16((N - 1) << 7, celt_log2tan(iside, imid));
1033 mid = imid / 32768.0f;
1034 side = iside / 32768.0f;
1036 /* This is a special case for N=2 that only works for stereo and takes
1037 advantage of the fact that mid and side are orthogonal to encode
1038 the side with just one bit. */
1039 if (N == 2 && dualstereo) {
1045 /* Only need one bit for the side */
1046 sbits = (itheta != 0 && itheta != 16384) ? 1 << 3 : 0;
1048 c = (itheta > 8192);
1049 f->remaining2 -= qalloc+sbits;
1054 sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
1055 ff_opus_rc_put_raw(rc, sign, 1);
1057 sign = 1 - 2 * sign;
1058 /* We use orig_fill here because we want to fold the side, but if
1059 itheta==16384, we'll have cleared the low bits of fill. */
1060 cm = ff_celt_encode_band(f, rc, band, x2, NULL, N, mbits, blocks,
1061 lowband, duration, lowband_out, level, gain,
1062 lowband_scratch, orig_fill);
1063 /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
1064 and there's no need to worry about mixing with the other channel. */
1065 y2[0] = -sign * x2[1];
1066 y2[1] = sign * x2[0];
1078 /* "Normal" split code */
1079 float *next_lowband2 = NULL;
1080 float *next_lowband_out1 = NULL;
1084 /* Give more bits to low-energy MDCTs than they would
1085 * otherwise deserve */
1086 if (B0 > 1 && !dualstereo && (itheta & 0x3fff)) {
1088 /* Rough approximation for pre-echo masking */
1089 delta -= delta >> (4 - duration);
1091 /* Corresponds to a forward-masking slope of
1092 * 1.5 dB per 10 ms */
1093 delta = FFMIN(0, delta + (N << 3 >> (5 - duration)));
1095 mbits = av_clip((b - delta) / 2, 0, b);
1097 f->remaining2 -= qalloc;
1099 if (lowband && !dualstereo)
1100 next_lowband2 = lowband + N; /* >32-bit split case */
1102 /* Only stereo needs to pass on lowband_out.
1103 * Otherwise, it's handled at the end */
1105 next_lowband_out1 = lowband_out;
1107 next_level = level + 1;
1109 rebalance = f->remaining2;
1110 if (mbits >= sbits) {
1111 /* In stereo mode, we do not apply a scaling to the mid
1112 * because we need the normalized mid for folding later */
1113 cm = ff_celt_encode_band(f, rc, band, X, NULL, N, mbits, blocks,
1114 lowband, duration, next_lowband_out1,
1115 next_level, dualstereo ? 1.0f : (gain * mid),
1116 lowband_scratch, fill);
1118 rebalance = mbits - (rebalance - f->remaining2);
1119 if (rebalance > 3 << 3 && itheta != 0)
1120 sbits += rebalance - (3 << 3);
1122 /* For a stereo split, the high bits of fill are always zero,
1123 * so no folding will be done to the side. */
1124 cm |= ff_celt_encode_band(f, rc, band, Y, NULL, N, sbits, blocks,
1125 next_lowband2, duration, NULL,
1126 next_level, gain * side, NULL,
1127 fill >> blocks) << ((B0 >> 1) & (dualstereo - 1));
1129 /* For a stereo split, the high bits of fill are always zero,
1130 * so no folding will be done to the side. */
1131 cm = ff_celt_encode_band(f, rc, band, Y, NULL, N, sbits, blocks,
1132 next_lowband2, duration, NULL,
1133 next_level, gain * side, NULL,
1134 fill >> blocks) << ((B0 >> 1) & (dualstereo - 1));
1136 rebalance = sbits - (rebalance - f->remaining2);
1137 if (rebalance > 3 << 3 && itheta != 16384)
1138 mbits += rebalance - (3 << 3);
1140 /* In stereo mode, we do not apply a scaling to the mid because
1141 * we need the normalized mid for folding later */
1142 cm |= ff_celt_encode_band(f, rc, band, X, NULL, N, mbits, blocks,
1143 lowband, duration, next_lowband_out1,
1144 next_level, dualstereo ? 1.0f : (gain * mid),
1145 lowband_scratch, fill);
1149 /* This is the basic no-split case */
1150 uint32_t q = celt_bits2pulses(cache, b);
1151 uint32_t curr_bits = celt_pulses2bits(cache, q);
1152 f->remaining2 -= curr_bits;
1154 /* Ensures we can never bust the budget */
1155 while (f->remaining2 < 0 && q > 0) {
1156 f->remaining2 += curr_bits;
1157 curr_bits = celt_pulses2bits(cache, --q);
1158 f->remaining2 -= curr_bits;
1162 /* Finally do the actual quantization */
1163 cm = celt_alg_quant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
1164 f->spread, blocks, gain);
1166 /* If there's no pulse, fill the band anyway */
1168 uint32_t cm_mask = (1 << blocks) - 1;
1171 for (j = 0; j < N; j++)
1176 for (j = 0; j < N; j++)
1177 X[j] = (((int32_t)celt_rng(f)) >> 20);
1180 /* Folded spectrum */
1181 for (j = 0; j < N; j++) {
1182 /* About 48 dB below the "normal" folding level */
1183 X[j] = lowband[j] + (((celt_rng(f)) & 0x8000) ? 1.0f / 256 : -1.0f / 256);
1187 celt_renormalize_vector(X, N, gain);
1192 /* This code is used by the decoder and by the resynthesis-enabled encoder */
1196 celt_stereo_merge(X, Y, mid, N);
1198 for (j = 0; j < N; j++)
1201 } else if (level == 0) {
1204 /* Undo the sample reorganization going from time order to frequency order */
1206 celt_interleave_hadamard(f->scratch, X, N_B >> recombine,
1207 B0<<recombine, longblocks);
1209 /* Undo time-freq changes that we did earlier */
1212 for (k = 0; k < time_divide; k++) {
1216 celt_haar1(X, N_B, blocks);
1219 for (k = 0; k < recombine; k++) {
1220 cm = ff_celt_bit_deinterleave[cm];
1221 celt_haar1(X, N0>>k, 1<<k);
1223 blocks <<= recombine;
1225 /* Scale output for later folding */
1228 float n = sqrtf(N0);
1229 for (j = 0; j < N0; j++)
1230 lowband_out[j] = n * X[j];
1232 cm = av_mod_uintp2(cm, blocks);
1238 float ff_celt_quant_band_cost(CeltFrame *f, OpusRangeCoder *rc, int band, float *bits,
1242 uint32_t cm[2] = { (1 << f->blocks) - 1, (1 << f->blocks) - 1 };
1243 const int band_size = ff_celt_freq_range[band] << f->size;
1244 float buf[352], lowband_scratch[176], norm1[176], norm2[176];
1245 float dist, cost, err_x = 0.0f, err_y = 0.0f;
1247 float *X_orig = f->block[0].coeffs + (ff_celt_freq_bands[band] << f->size);
1248 float *Y = (f->channels == 2) ? &buf[176] : NULL;
1249 float *Y_orig = f->block[1].coeffs + (ff_celt_freq_bands[band] << f->size);
1250 OPUS_RC_CHECKPOINT_SPAWN(rc);
1252 memcpy(X, X_orig, band_size*sizeof(float));
1254 memcpy(Y, Y_orig, band_size*sizeof(float));
1256 f->remaining2 = ((f->framebits << 3) - f->anticollapse_needed) - opus_rc_tell_frac(rc) - 1;
1257 if (band <= f->coded_bands - 1) {
1258 int curr_balance = f->remaining / FFMIN(3, f->coded_bands - band);
1259 b = av_clip_uintp2(FFMIN(f->remaining2 + 1, f->pulses[band] + curr_balance), 14);
1262 if (f->dual_stereo) {
1263 ff_celt_encode_band(f, rc, band, X, NULL, band_size, b / 2, f->blocks, NULL,
1264 f->size, norm1, 0, 1.0f, lowband_scratch, cm[0]);
1266 ff_celt_encode_band(f, rc, band, Y, NULL, band_size, b / 2, f->blocks, NULL,
1267 f->size, norm2, 0, 1.0f, lowband_scratch, cm[1]);
1269 ff_celt_encode_band(f, rc, band, X, Y, band_size, b, f->blocks, NULL, f->size,
1270 norm1, 0, 1.0f, lowband_scratch, cm[0] | cm[1]);
1273 for (i = 0; i < band_size; i++) {
1274 err_x += (X[i] - X_orig[i])*(X[i] - X_orig[i]);
1275 err_y += (Y[i] - Y_orig[i])*(Y[i] - Y_orig[i]);
1278 dist = sqrtf(err_x) + sqrtf(err_y);
1279 cost = OPUS_RC_CHECKPOINT_BITS(rc)/8.0f;
1282 OPUS_RC_CHECKPOINT_ROLLBACK(rc);
1284 return lambda*dist*cost;