opus_pvq: minor cleanups
[ffmpeg.git] / libavcodec / opus_pvq.c
1 /*
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>
8  *
9  * This file is part of FFmpeg.
10  *
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.
15  *
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.
20  *
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
24  */
25
26 #include "opustab.h"
27 #include "opus_pvq.h"
28
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))
31
32 static inline int16_t celt_cos(int16_t x)
33 {
34     x = (MUL16(x, x) + 4096) >> 13;
35     x = (32767-x) + ROUND_MUL16(x, (-7651 + ROUND_MUL16(x, (8277 + ROUND_MUL16(-626, x)))));
36     return x + 1;
37 }
38
39 static inline int celt_log2tan(int isin, int icos)
40 {
41     int lc, ls;
42     lc = opus_ilog(icos);
43     ls = opus_ilog(isin);
44     icos <<= 15 - lc;
45     isin <<= 15 - ls;
46     return (ls << 11) - (lc << 11) +
47            ROUND_MUL16(isin, ROUND_MUL16(isin, -2597) + 7932) -
48            ROUND_MUL16(icos, ROUND_MUL16(icos, -2597) + 7932);
49 }
50
51 static inline int celt_bits2pulses(const uint8_t *cache, int bits)
52 {
53     // TODO: Find the size of cache and make it into an array in the parameters list
54     int i, low = 0, high;
55
56     high = cache[0];
57     bits--;
58
59     for (i = 0; i < 6; i++) {
60         int center = (low + high + 1) >> 1;
61         if (cache[center] >= bits)
62             high = center;
63         else
64             low = center;
65     }
66
67     return (bits - (low == 0 ? -1 : cache[low]) <= cache[high] - bits) ? low : high;
68 }
69
70 static inline int celt_pulses2bits(const uint8_t *cache, int pulses)
71 {
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;
74 }
75
76 static inline void celt_normalize_residual(const int * av_restrict iy, float * av_restrict X,
77                                            int N, float g)
78 {
79     int i;
80     for (i = 0; i < N; i++)
81         X[i] = g * iy[i];
82 }
83
84 static void celt_exp_rotation_impl(float *X, uint32_t len, uint32_t stride,
85                                    float c, float s)
86 {
87     float *Xptr;
88     int i;
89
90     Xptr = X;
91     for (i = 0; i < len - stride; i++) {
92         float x1, x2;
93         x1           = Xptr[0];
94         x2           = Xptr[stride];
95         Xptr[stride] = c * x2 + s * x1;
96         *Xptr++      = c * x1 - s * x2;
97     }
98
99     Xptr = &X[len - 2 * stride - 1];
100     for (i = len - 2 * stride - 1; i >= 0; i--) {
101         float x1, x2;
102         x1           = Xptr[0];
103         x2           = Xptr[stride];
104         Xptr[stride] = c * x2 + s * x1;
105         *Xptr--      = c * x1 - s * x2;
106     }
107 }
108
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)
112 {
113     uint32_t stride2 = 0;
114     float c, s;
115     float gain, theta;
116     int i;
117
118     if (2*K >= len || spread == CELT_SPREAD_NONE)
119         return;
120
121     gain = (float)len / (len + (20 - 5*spread) * K);
122     theta = M_PI * gain * gain / 4;
123
124     c = cosf(theta);
125     s = sinf(theta);
126
127     if (len >= stride << 3) {
128         stride2 = 1;
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)
132             stride2++;
133     }
134
135     /*NOTE: As a minor optimization, we could be passing around log2(B), not B, for both this and for
136     extract_collapse_mask().*/
137     len /= stride;
138     for (i = 0; i < stride; i++) {
139         if (encode) {
140             celt_exp_rotation_impl(X + i * len, len, 1, c, -s);
141             if (stride2)
142                 celt_exp_rotation_impl(X + i * len, len, stride2, s, -c);
143         } else {
144             if (stride2)
145                 celt_exp_rotation_impl(X + i * len, len, stride2, s, c);
146             celt_exp_rotation_impl(X + i * len, len, 1, c, s);
147         }
148     }
149 }
150
151 static inline uint32_t celt_extract_collapse_mask(const int *iy, uint32_t N, uint32_t B)
152 {
153     uint32_t collapse_mask;
154     int N0;
155     int i, j;
156
157     if (B <= 1)
158         return 1;
159
160     /*NOTE: As a minor optimization, we could be passing around log2(B), not B, for both this and for
161     exp_rotation().*/
162     N0 = N/B;
163     collapse_mask = 0;
164     for (i = 0; i < B; i++)
165         for (j = 0; j < N0; j++)
166             collapse_mask |= (!!iy[i*N0+j]) << i;
167     return collapse_mask;
168 }
169
170 static inline void celt_stereo_merge(float *X, float *Y, float mid, int N)
171 {
172     int i;
173     float xp = 0, side = 0;
174     float E[2];
175     float mid2;
176     float gain[2];
177
178     /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
179     for (i = 0; i < N; i++) {
180         xp   += X[i] * Y[i];
181         side += Y[i] * Y[i];
182     }
183
184     /* Compensating for the mid normalization */
185     xp *= mid;
186     mid2 = mid;
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++)
191             Y[i] = X[i];
192         return;
193     }
194
195     gain[0] = 1.0f / sqrtf(E[0]);
196     gain[1] = 1.0f / sqrtf(E[1]);
197
198     for (i = 0; i < N; i++) {
199         float value[2];
200         /* Apply mid scaling (side is already scaled) */
201         value[0] = mid * X[i];
202         value[1] = Y[i];
203         X[i] = gain[0] * (value[0] - value[1]);
204         Y[i] = gain[1] * (value[0] + value[1]);
205     }
206 }
207
208 static void celt_interleave_hadamard(float *tmp, float *X, int N0,
209                                      int stride, int hadamard)
210 {
211     int i, j, N = N0*stride;
212     const uint8_t *order = &ff_celt_hadamard_order[hadamard ? stride - 2 : 30];
213
214     for (i = 0; i < stride; i++)
215         for (j = 0; j < N0; j++)
216             tmp[j*stride+i] = X[order[i]*N0+j];
217
218     memcpy(X, tmp, N*sizeof(float));
219 }
220
221 static void celt_deinterleave_hadamard(float *tmp, float *X, int N0,
222                                        int stride, int hadamard)
223 {
224     int i, j, N = N0*stride;
225     const uint8_t *order = &ff_celt_hadamard_order[hadamard ? stride - 2 : 30];
226
227     for (i = 0; i < stride; i++)
228         for (j = 0; j < N0; j++)
229             tmp[order[i]*N0+j] = X[j*stride+i];
230
231     memcpy(X, tmp, N*sizeof(float));
232 }
233
234 static void celt_haar1(float *X, int N0, int stride)
235 {
236     int i, j;
237     N0 >>= 1;
238     for (i = 0; i < stride; i++) {
239         for (j = 0; j < N0; j++) {
240             float x0 = X[stride * (2 * j + 0) + i];
241             float x1 = X[stride * (2 * j + 1) + i];
242             X[stride * (2 * j + 0) + i] = (x0 + x1) * M_SQRT1_2;
243             X[stride * (2 * j + 1) + i] = (x0 - x1) * M_SQRT1_2;
244         }
245     }
246 }
247
248 static inline int celt_compute_qn(int N, int b, int offset, int pulse_cap,
249                                   int stereo)
250 {
251     int qn, qb;
252     int N2 = 2 * N - 1;
253     if (stereo && N == 2)
254         N2--;
255
256     /* The upper limit ensures that in a stereo split with itheta==16384, we'll
257      * always have enough bits left over to code at least one pulse in the
258      * side; otherwise it would collapse, since it doesn't get folded. */
259     qb = FFMIN3(b - pulse_cap - (4 << 3), (b + N2 * offset) / N2, 8 << 3);
260     qn = (qb < (1 << 3 >> 1)) ? 1 : ((ff_celt_qn_exp2[qb & 0x7] >> (14 - (qb >> 3))) + 1) >> 1 << 1;
261     return qn;
262 }
263
264 /* Convert the quantized vector to an index */
265 static inline uint32_t celt_icwrsi(uint32_t N, uint32_t K, const int *y)
266 {
267     int i, idx = 0, sum = 0;
268     for (i = N - 1; i >= 0; i--) {
269         const uint32_t i_s = CELT_PVQ_U(N - i, sum + FFABS(y[i]) + 1);
270         idx += CELT_PVQ_U(N - i, sum) + (y[i] < 0)*i_s;
271         sum += FFABS(y[i]);
272     }
273     return idx;
274 }
275
276 // this code was adapted from libopus
277 static inline uint64_t celt_cwrsi(uint32_t N, uint32_t K, uint32_t i, int *y)
278 {
279     uint64_t norm = 0;
280     uint32_t p;
281     int s, val;
282     int k0;
283
284     while (N > 2) {
285         uint32_t q;
286
287         /*Lots of pulses case:*/
288         if (K >= N) {
289             const uint32_t *row = ff_celt_pvq_u_row[N];
290
291             /* Are the pulses in this dimension negative? */
292             p  = row[K + 1];
293             s  = -(i >= p);
294             i -= p & s;
295
296             /*Count how many pulses were placed in this dimension.*/
297             k0 = K;
298             q = row[N];
299             if (q > i) {
300                 K = N;
301                 do {
302                     p = ff_celt_pvq_u_row[--K][N];
303                 } while (p > i);
304             } else
305                 for (p = row[K]; p > i; p = row[K])
306                     K--;
307
308             i    -= p;
309             val   = (k0 - K + s) ^ s;
310             norm += val * val;
311             *y++  = val;
312         } else { /*Lots of dimensions case:*/
313             /*Are there any pulses in this dimension at all?*/
314             p = ff_celt_pvq_u_row[K    ][N];
315             q = ff_celt_pvq_u_row[K + 1][N];
316
317             if (p <= i && i < q) {
318                 i -= p;
319                 *y++ = 0;
320             } else {
321                 /*Are the pulses in this dimension negative?*/
322                 s  = -(i >= q);
323                 i -= q & s;
324
325                 /*Count how many pulses were placed in this dimension.*/
326                 k0 = K;
327                 do p = ff_celt_pvq_u_row[--K][N];
328                 while (p > i);
329
330                 i    -= p;
331                 val   = (k0 - K + s) ^ s;
332                 norm += val * val;
333                 *y++  = val;
334             }
335         }
336         N--;
337     }
338
339     /* N == 2 */
340     p  = 2 * K + 1;
341     s  = -(i >= p);
342     i -= p & s;
343     k0 = K;
344     K  = (i + 1) / 2;
345
346     if (K)
347         i -= 2 * K - 1;
348
349     val   = (k0 - K + s) ^ s;
350     norm += val * val;
351     *y++  = val;
352
353     /* N==1 */
354     s     = -i;
355     val   = (K + s) ^ s;
356     norm += val * val;
357     *y    = val;
358
359     return norm;
360 }
361
362 static inline void celt_encode_pulses(OpusRangeCoder *rc, int *y, uint32_t N, uint32_t K)
363 {
364     ff_opus_rc_enc_uint(rc, celt_icwrsi(N, K, y), CELT_PVQ_V(N, K));
365 }
366
367 static inline float celt_decode_pulses(OpusRangeCoder *rc, int *y, uint32_t N, uint32_t K)
368 {
369     const uint32_t idx = ff_opus_rc_dec_uint(rc, CELT_PVQ_V(N, K));
370     return celt_cwrsi(N, K, idx, y);
371 }
372
373 /*
374  * Faster than libopus's search, operates entirely in the signed domain.
375  * Slightly worse/better depending on N, K and the input vector.
376  */
377 static int celt_pvq_search(float *X, int *y, int K, int N)
378 {
379     int i, y_norm = 0;
380     float res = 0.0f, xy_norm = 0.0f;
381
382     for (i = 0; i < N; i++)
383         res += FFABS(X[i]);
384
385     res = K/(res + FLT_EPSILON);
386
387     for (i = 0; i < N; i++) {
388         y[i] = lrintf(res*X[i]);
389         y_norm  += y[i]*y[i];
390         xy_norm += y[i]*X[i];
391         K -= FFABS(y[i]);
392     }
393
394     while (K) {
395         int max_idx = 0, max_den = 1, phase = FFSIGN(K);
396         float max_num = 0.0f;
397         y_norm += 1.0f;
398
399         for (i = 0; i < N; i++) {
400             /* If the sum has been overshot and the best place has 0 pulses allocated
401              * to it, attempting to decrease it further will actually increase the
402              * sum. Prevent this by disregarding any 0 positions when decrementing. */
403             const int ca = 1 ^ ((y[i] == 0) & (phase < 0));
404             const int y_new = y_norm  + 2*phase*FFABS(y[i]);
405             float xy_new = xy_norm + 1*phase*FFABS(X[i]);
406             xy_new = xy_new * xy_new;
407             if (ca && (max_den*xy_new) > (y_new*max_num)) {
408                 max_den = y_new;
409                 max_num = xy_new;
410                 max_idx = i;
411             }
412         }
413
414         K -= phase;
415
416         phase *= FFSIGN(X[max_idx]);
417         xy_norm += 1*phase*X[max_idx];
418         y_norm  += 2*phase*y[max_idx];
419         y[max_idx] += phase;
420     }
421
422     return y_norm;
423 }
424
425 static uint32_t celt_alg_quant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K,
426                                enum CeltSpread spread, uint32_t blocks, float gain)
427 {
428     int y[176];
429
430     celt_exp_rotation(X, N, blocks, K, spread, 1);
431     gain /= sqrtf(celt_pvq_search(X, y, K, N));
432     celt_encode_pulses(rc, y,  N, K);
433     celt_normalize_residual(y, X, N, gain);
434     celt_exp_rotation(X, N, blocks, K, spread, 0);
435     return celt_extract_collapse_mask(y, N, blocks);
436 }
437
438 /** Decode pulse vector and combine the result with the pitch vector to produce
439     the final normalised signal in the current band. */
440 static uint32_t celt_alg_unquant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K,
441                                  enum CeltSpread spread, uint32_t blocks, float gain)
442 {
443     int y[176];
444
445     gain /= sqrtf(celt_decode_pulses(rc, y, N, K));
446     celt_normalize_residual(y, X, N, gain);
447     celt_exp_rotation(X, N, blocks, K, spread, 0);
448     return celt_extract_collapse_mask(y, N, blocks);
449 }
450
451 uint32_t ff_celt_decode_band(CeltFrame *f, OpusRangeCoder *rc, const int band,
452                              float *X, float *Y, int N, int b, uint32_t blocks,
453                              float *lowband, int duration, float *lowband_out, int level,
454                              float gain, float *lowband_scratch, int fill)
455 {
456     int i;
457     const uint8_t *cache;
458     int stereo = !!Y, split = !!Y;
459     int imid = 0, iside = 0;
460     uint32_t N0 = N;
461     int N_B = N / blocks;
462     int N_B0 = N_B;
463     int B0 = blocks;
464     int time_divide = 0;
465     int recombine = 0;
466     int inv = 0;
467     float mid = 0, side = 0;
468     int longblocks = (B0 == 1);
469     uint32_t cm = 0;
470
471     if (N == 1) {
472         /* special case for one sample */
473         float *x = X;
474         for (i = 0; i <= stereo; i++) {
475             int sign = 0;
476             if (f->remaining2 >= 1<<3) {
477                 sign           = ff_opus_rc_get_raw(rc, 1);
478                 f->remaining2 -= 1 << 3;
479                 b             -= 1 << 3;
480             }
481             x[0] = sign ? -1.0f : 1.0f;
482             x = Y;
483         }
484         if (lowband_out)
485             lowband_out[0] = X[0];
486         return 1;
487     }
488
489     if (!stereo && level == 0) {
490         int tf_change = f->tf_change[band];
491         int k;
492         if (tf_change > 0)
493             recombine = tf_change;
494         /* Band recombining to increase frequency resolution */
495
496         if (lowband &&
497             (recombine || ((N_B & 1) == 0 && tf_change < 0) || B0 > 1)) {
498             for (i = 0; i < N; i++)
499                 lowband_scratch[i] = lowband[i];
500             lowband = lowband_scratch;
501         }
502
503         for (k = 0; k < recombine; k++) {
504             if (lowband)
505                 celt_haar1(lowband, N >> k, 1 << k);
506             fill = ff_celt_bit_interleave[fill & 0xF] | ff_celt_bit_interleave[fill >> 4] << 2;
507         }
508         blocks >>= recombine;
509         N_B <<= recombine;
510
511         /* Increasing the time resolution */
512         while ((N_B & 1) == 0 && tf_change < 0) {
513             if (lowband)
514                 celt_haar1(lowband, N_B, blocks);
515             fill |= fill << blocks;
516             blocks <<= 1;
517             N_B >>= 1;
518             time_divide++;
519             tf_change++;
520         }
521         B0 = blocks;
522         N_B0 = N_B;
523
524         /* Reorganize the samples in time order instead of frequency order */
525         if (B0 > 1 && lowband)
526             celt_deinterleave_hadamard(f->scratch, lowband, N_B >> recombine,
527                                        B0 << recombine, longblocks);
528     }
529
530     /* If we need 1.5 more bit than we can produce, split the band in two. */
531     cache = ff_celt_cache_bits +
532             ff_celt_cache_index[(duration + 1) * CELT_MAX_BANDS + band];
533     if (!stereo && duration >= 0 && b > cache[cache[0]] + 12 && N > 2) {
534         N >>= 1;
535         Y = X + N;
536         split = 1;
537         duration -= 1;
538         if (blocks == 1)
539             fill = (fill & 1) | (fill << 1);
540         blocks = (blocks + 1) >> 1;
541     }
542
543     if (split) {
544         int qn;
545         int itheta = 0;
546         int mbits, sbits, delta;
547         int qalloc;
548         int pulse_cap;
549         int offset;
550         int orig_fill;
551         int tell;
552
553         /* Decide on the resolution to give to the split parameter theta */
554         pulse_cap = ff_celt_log_freq_range[band] + duration * 8;
555         offset = (pulse_cap >> 1) - (stereo && N == 2 ? CELT_QTHETA_OFFSET_TWOPHASE :
556                                                           CELT_QTHETA_OFFSET);
557         qn = (stereo && band >= f->intensity_stereo) ? 1 :
558              celt_compute_qn(N, b, offset, pulse_cap, stereo);
559         tell = opus_rc_tell_frac(rc);
560         if (qn != 1) {
561             /* Entropy coding of the angle. We use a uniform pdf for the
562             time split, a step for stereo, and a triangular one for the rest. */
563             if (stereo && N > 2)
564                 itheta = ff_opus_rc_dec_uint_step(rc, qn/2);
565             else if (stereo || B0 > 1)
566                 itheta = ff_opus_rc_dec_uint(rc, qn+1);
567             else
568                 itheta = ff_opus_rc_dec_uint_tri(rc, qn);
569             itheta = itheta * 16384 / qn;
570             /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.
571             Let's do that at higher complexity */
572         } else if (stereo) {
573             inv = (b > 2 << 3 && f->remaining2 > 2 << 3) ? ff_opus_rc_dec_log(rc, 2) : 0;
574             itheta = 0;
575         }
576         qalloc = opus_rc_tell_frac(rc) - tell;
577         b -= qalloc;
578
579         orig_fill = fill;
580         if (itheta == 0) {
581             imid = 32767;
582             iside = 0;
583             fill = av_mod_uintp2(fill, blocks);
584             delta = -16384;
585         } else if (itheta == 16384) {
586             imid = 0;
587             iside = 32767;
588             fill &= ((1 << blocks) - 1) << blocks;
589             delta = 16384;
590         } else {
591             imid = celt_cos(itheta);
592             iside = celt_cos(16384-itheta);
593             /* This is the mid vs side allocation that minimizes squared error
594             in that band. */
595             delta = ROUND_MUL16((N - 1) << 7, celt_log2tan(iside, imid));
596         }
597
598         mid  = imid  / 32768.0f;
599         side = iside / 32768.0f;
600
601         /* This is a special case for N=2 that only works for stereo and takes
602         advantage of the fact that mid and side are orthogonal to encode
603         the side with just one bit. */
604         if (N == 2 && stereo) {
605             int c;
606             int sign = 0;
607             float tmp;
608             float *x2, *y2;
609             mbits = b;
610             /* Only need one bit for the side */
611             sbits = (itheta != 0 && itheta != 16384) ? 1 << 3 : 0;
612             mbits -= sbits;
613             c = (itheta > 8192);
614             f->remaining2 -= qalloc+sbits;
615
616             x2 = c ? Y : X;
617             y2 = c ? X : Y;
618             if (sbits)
619                 sign = ff_opus_rc_get_raw(rc, 1);
620             sign = 1 - 2 * sign;
621             /* We use orig_fill here because we want to fold the side, but if
622             itheta==16384, we'll have cleared the low bits of fill. */
623             cm = ff_celt_decode_band(f, rc, band, x2, NULL, N, mbits, blocks,
624                                      lowband, duration, lowband_out, level, gain,
625                                      lowband_scratch, orig_fill);
626             /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
627             and there's no need to worry about mixing with the other channel. */
628             y2[0] = -sign * x2[1];
629             y2[1] =  sign * x2[0];
630             X[0] *= mid;
631             X[1] *= mid;
632             Y[0] *= side;
633             Y[1] *= side;
634             tmp = X[0];
635             X[0] = tmp - Y[0];
636             Y[0] = tmp + Y[0];
637             tmp = X[1];
638             X[1] = tmp - Y[1];
639             Y[1] = tmp + Y[1];
640         } else {
641             /* "Normal" split code */
642             float *next_lowband2     = NULL;
643             float *next_lowband_out1 = NULL;
644             int next_level = 0;
645             int rebalance;
646
647             /* Give more bits to low-energy MDCTs than they would
648              * otherwise deserve */
649             if (B0 > 1 && !stereo && (itheta & 0x3fff)) {
650                 if (itheta > 8192)
651                     /* Rough approximation for pre-echo masking */
652                     delta -= delta >> (4 - duration);
653                 else
654                     /* Corresponds to a forward-masking slope of
655                      * 1.5 dB per 10 ms */
656                     delta = FFMIN(0, delta + (N << 3 >> (5 - duration)));
657             }
658             mbits = av_clip((b - delta) / 2, 0, b);
659             sbits = b - mbits;
660             f->remaining2 -= qalloc;
661
662             if (lowband && !stereo)
663                 next_lowband2 = lowband + N; /* >32-bit split case */
664
665             /* Only stereo needs to pass on lowband_out.
666              * Otherwise, it's handled at the end */
667             if (stereo)
668                 next_lowband_out1 = lowband_out;
669             else
670                 next_level = level + 1;
671
672             rebalance = f->remaining2;
673             if (mbits >= sbits) {
674                 /* In stereo mode, we do not apply a scaling to the mid
675                  * because we need the normalized mid for folding later */
676                 cm = ff_celt_decode_band(f, rc, band, X, NULL, N, mbits, blocks,
677                                          lowband, duration, next_lowband_out1,
678                                          next_level, stereo ? 1.0f : (gain * mid),
679                                          lowband_scratch, fill);
680
681                 rebalance = mbits - (rebalance - f->remaining2);
682                 if (rebalance > 3 << 3 && itheta != 0)
683                     sbits += rebalance - (3 << 3);
684
685                 /* For a stereo split, the high bits of fill are always zero,
686                  * so no folding will be done to the side. */
687                 cm |= ff_celt_decode_band(f, rc, band, Y, NULL, N, sbits, blocks,
688                                           next_lowband2, duration, NULL,
689                                           next_level, gain * side, NULL,
690                                           fill >> blocks) << ((B0 >> 1) & (stereo - 1));
691             } else {
692                 /* For a stereo split, the high bits of fill are always zero,
693                  * so no folding will be done to the side. */
694                 cm = ff_celt_decode_band(f, rc, band, Y, NULL, N, sbits, blocks,
695                                          next_lowband2, duration, NULL,
696                                          next_level, gain * side, NULL,
697                                          fill >> blocks) << ((B0 >> 1) & (stereo - 1));
698
699                 rebalance = sbits - (rebalance - f->remaining2);
700                 if (rebalance > 3 << 3 && itheta != 16384)
701                     mbits += rebalance - (3 << 3);
702
703                 /* In stereo mode, we do not apply a scaling to the mid because
704                  * we need the normalized mid for folding later */
705                 cm |= ff_celt_decode_band(f, rc, band, X, NULL, N, mbits, blocks,
706                                           lowband, duration, next_lowband_out1,
707                                           next_level, stereo ? 1.0f : (gain * mid),
708                                           lowband_scratch, fill);
709             }
710         }
711     } else {
712         /* This is the basic no-split case */
713         uint32_t q         = celt_bits2pulses(cache, b);
714         uint32_t curr_bits = celt_pulses2bits(cache, q);
715         f->remaining2 -= curr_bits;
716
717         /* Ensures we can never bust the budget */
718         while (f->remaining2 < 0 && q > 0) {
719             f->remaining2 += curr_bits;
720             curr_bits      = celt_pulses2bits(cache, --q);
721             f->remaining2 -= curr_bits;
722         }
723
724         if (q != 0) {
725             /* Finally do the actual quantization */
726             cm = celt_alg_unquant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
727                                   f->spread, blocks, gain);
728         } else {
729             /* If there's no pulse, fill the band anyway */
730             uint32_t cm_mask = (1 << blocks) - 1;
731             fill &= cm_mask;
732             if (fill) {
733                 if (!lowband) {
734                     /* Noise */
735                     for (i = 0; i < N; i++)
736                         X[i] = (((int32_t)celt_rng(f)) >> 20);
737                     cm = cm_mask;
738                 } else {
739                     /* Folded spectrum */
740                     for (i = 0; i < N; i++) {
741                         /* About 48 dB below the "normal" folding level */
742                         X[i] = lowband[i] + (((celt_rng(f)) & 0x8000) ? 1.0f / 256 : -1.0f / 256);
743                     }
744                     cm = fill;
745                 }
746                 celt_renormalize_vector(X, N, gain);
747             } else {
748                 memset(X, 0, N*sizeof(float));
749             }
750         }
751     }
752
753     /* This code is used by the decoder and by the resynthesis-enabled encoder */
754     if (stereo) {
755         if (N > 2)
756             celt_stereo_merge(X, Y, mid, N);
757         if (inv) {
758             for (i = 0; i < N; i++)
759                 Y[i] *= -1;
760         }
761     } else if (level == 0) {
762         int k;
763
764         /* Undo the sample reorganization going from time order to frequency order */
765         if (B0 > 1)
766             celt_interleave_hadamard(f->scratch, X, N_B >> recombine,
767                                      B0 << recombine, longblocks);
768
769         /* Undo time-freq changes that we did earlier */
770         N_B = N_B0;
771         blocks = B0;
772         for (k = 0; k < time_divide; k++) {
773             blocks >>= 1;
774             N_B <<= 1;
775             cm |= cm >> blocks;
776             celt_haar1(X, N_B, blocks);
777         }
778
779         for (k = 0; k < recombine; k++) {
780             cm = ff_celt_bit_deinterleave[cm];
781             celt_haar1(X, N0>>k, 1<<k);
782         }
783         blocks <<= recombine;
784
785         /* Scale output for later folding */
786         if (lowband_out) {
787             float n = sqrtf(N0);
788             for (i = 0; i < N0; i++)
789                 lowband_out[i] = n * X[i];
790         }
791         cm = av_mod_uintp2(cm, blocks);
792     }
793
794     return cm;
795 }
796
797 /* This has to be, AND MUST BE done by the psychoacoustic system, this has a very
798  * big impact on the entire quantization and especially huge on transients */
799 static int celt_calc_theta(const float *X, const float *Y, int coupling, int N)
800 {
801     int i;
802     float e[2] = { 0.0f, 0.0f };
803     if (coupling) { /* Coupling case */
804         for (i = 0; i < N; i++) {
805             e[0] += (X[i] + Y[i])*(X[i] + Y[i]);
806             e[1] += (X[i] - Y[i])*(X[i] - Y[i]);
807         }
808     } else {
809         for (i = 0; i < N; i++) {
810             e[0] += X[i]*X[i];
811             e[1] += Y[i]*Y[i];
812         }
813     }
814     return lrintf(32768.0f*atan2f(sqrtf(e[1]), sqrtf(e[0]))/M_PI);
815 }
816
817 static void celt_stereo_is_decouple(float *X, float *Y, float e_l, float e_r, int N)
818 {
819     int i;
820     const float energy_n = 1.0f/(sqrtf(e_l*e_l + e_r*e_r) + FLT_EPSILON);
821     e_l *= energy_n;
822     e_r *= energy_n;
823     for (i = 0; i < N; i++)
824         X[i] = e_l*X[i] + e_r*Y[i];
825 }
826
827 static void celt_stereo_ms_decouple(float *X, float *Y, int N)
828 {
829     int i;
830     for (i = 0; i < N; i++) {
831         const float Xret = X[i];
832         X[i] = (X[i] + Y[i])*M_SQRT1_2;
833         Y[i] = (Y[i] - Xret)*M_SQRT1_2;
834     }
835 }
836
837 uint32_t ff_celt_encode_band(CeltFrame *f, OpusRangeCoder *rc, const int band,
838                              float *X, float *Y, int N, int b, uint32_t blocks,
839                              float *lowband, int duration, float *lowband_out, int level,
840                              float gain, float *lowband_scratch, int fill)
841 {
842     int i;
843     const uint8_t *cache;
844     int stereo = !!Y, split = !!Y;
845     int imid = 0, iside = 0;
846     uint32_t N0 = N;
847     int N_B = N / blocks;
848     int N_B0 = N_B;
849     int B0 = blocks;
850     int time_divide = 0;
851     int recombine = 0;
852     int inv = 0;
853     float mid = 0, side = 0;
854     int longblocks = (B0 == 1);
855     uint32_t cm = 0;
856
857     if (N == 1) {
858         /* special case for one sample - the decoder's output will be +- 1.0f!!! */
859         float *x = X;
860         for (i = 0; i <= stereo; i++) {
861             if (f->remaining2 >= 1<<3) {
862                 ff_opus_rc_put_raw(rc, x[0] < 0, 1);
863                 f->remaining2 -= 1 << 3;
864                 b             -= 1 << 3;
865             }
866             x[0] = 1.0f - 2.0f*(x[0] < 0);
867             x = Y;
868         }
869         if (lowband_out)
870             lowband_out[0] = X[0];
871         return 1;
872     }
873
874     if (!stereo && level == 0) {
875         int tf_change = f->tf_change[band];
876         int k;
877         if (tf_change > 0)
878             recombine = tf_change;
879         /* Band recombining to increase frequency resolution */
880
881         if (lowband &&
882             (recombine || ((N_B & 1) == 0 && tf_change < 0) || B0 > 1)) {
883             for (i = 0; i < N; i++)
884                 lowband_scratch[i] = lowband[i];
885             lowband = lowband_scratch;
886         }
887
888         for (k = 0; k < recombine; k++) {
889             celt_haar1(X, N >> k, 1 << k);
890             fill = ff_celt_bit_interleave[fill & 0xF] | ff_celt_bit_interleave[fill >> 4] << 2;
891         }
892         blocks >>= recombine;
893         N_B <<= recombine;
894
895         /* Increasing the time resolution */
896         while ((N_B & 1) == 0 && tf_change < 0) {
897             celt_haar1(X, N_B, blocks);
898             fill |= fill << blocks;
899             blocks <<= 1;
900             N_B >>= 1;
901             time_divide++;
902             tf_change++;
903         }
904         B0 = blocks;
905         N_B0 = N_B;
906
907         /* Reorganize the samples in time order instead of frequency order */
908         if (B0 > 1)
909             celt_deinterleave_hadamard(f->scratch, X, N_B >> recombine,
910                                        B0 << recombine, longblocks);
911     }
912
913     /* If we need 1.5 more bit than we can produce, split the band in two. */
914     cache = ff_celt_cache_bits +
915             ff_celt_cache_index[(duration + 1) * CELT_MAX_BANDS + band];
916     if (!stereo && duration >= 0 && b > cache[cache[0]] + 12 && N > 2) {
917         N >>= 1;
918         Y = X + N;
919         split = 1;
920         duration -= 1;
921         if (blocks == 1)
922             fill = (fill & 1) | (fill << 1);
923         blocks = (blocks + 1) >> 1;
924     }
925
926     if (split) {
927         int qn;
928         int itheta = celt_calc_theta(X, Y, stereo, N);
929         int mbits, sbits, delta;
930         int qalloc;
931         int pulse_cap;
932         int offset;
933         int orig_fill;
934         int tell;
935
936         /* Decide on the resolution to give to the split parameter theta */
937         pulse_cap = ff_celt_log_freq_range[band] + duration * 8;
938         offset = (pulse_cap >> 1) - (stereo && N == 2 ? CELT_QTHETA_OFFSET_TWOPHASE :
939                                                           CELT_QTHETA_OFFSET);
940         qn = (stereo && band >= f->intensity_stereo) ? 1 :
941              celt_compute_qn(N, b, offset, pulse_cap, stereo);
942         tell = opus_rc_tell_frac(rc);
943
944         if (qn != 1) {
945
946             itheta = (itheta*qn + 8192) >> 14;
947
948             /* Entropy coding of the angle. We use a uniform pdf for the
949              * time split, a step for stereo, and a triangular one for the rest. */
950             if (stereo && N > 2)
951                 ff_opus_rc_enc_uint_step(rc, itheta, qn / 2);
952             else if (stereo || B0 > 1)
953                 ff_opus_rc_enc_uint(rc, itheta, qn + 1);
954             else
955                 ff_opus_rc_enc_uint_tri(rc, itheta, qn);
956             itheta = itheta * 16384 / qn;
957
958             if (stereo) {
959                 if (itheta == 0)
960                     celt_stereo_is_decouple(X, Y, f->block[0].lin_energy[band],
961                                             f->block[1].lin_energy[band], N);
962                 else
963                     celt_stereo_ms_decouple(X, Y, N);
964             }
965         } else if (stereo) {
966              inv = itheta > 8192;
967              if (inv) {
968                 for (i = 0; i < N; i++)
969                    Y[i] *= -1;
970              }
971              celt_stereo_is_decouple(X, Y, f->block[0].lin_energy[band],
972                                      f->block[1].lin_energy[band], N);
973
974             if (b > 2 << 3 && f->remaining2 > 2 << 3) {
975                 ff_opus_rc_enc_log(rc, inv, 2);
976             } else {
977                 inv = 0;
978             }
979
980             itheta = 0;
981         }
982         qalloc = opus_rc_tell_frac(rc) - tell;
983         b -= qalloc;
984
985         orig_fill = fill;
986         if (itheta == 0) {
987             imid = 32767;
988             iside = 0;
989             fill = av_mod_uintp2(fill, blocks);
990             delta = -16384;
991         } else if (itheta == 16384) {
992             imid = 0;
993             iside = 32767;
994             fill &= ((1 << blocks) - 1) << blocks;
995             delta = 16384;
996         } else {
997             imid = celt_cos(itheta);
998             iside = celt_cos(16384-itheta);
999             /* This is the mid vs side allocation that minimizes squared error
1000             in that band. */
1001             delta = ROUND_MUL16((N - 1) << 7, celt_log2tan(iside, imid));
1002         }
1003
1004         mid  = imid  / 32768.0f;
1005         side = iside / 32768.0f;
1006
1007         /* This is a special case for N=2 that only works for stereo and takes
1008         advantage of the fact that mid and side are orthogonal to encode
1009         the side with just one bit. */
1010         if (N == 2 && stereo) {
1011             int c;
1012             int sign = 0;
1013             float tmp;
1014             float *x2, *y2;
1015             mbits = b;
1016             /* Only need one bit for the side */
1017             sbits = (itheta != 0 && itheta != 16384) ? 1 << 3 : 0;
1018             mbits -= sbits;
1019             c = (itheta > 8192);
1020             f->remaining2 -= qalloc+sbits;
1021
1022             x2 = c ? Y : X;
1023             y2 = c ? X : Y;
1024             if (sbits) {
1025                 sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
1026                 ff_opus_rc_put_raw(rc, sign, 1);
1027             }
1028             sign = 1 - 2 * sign;
1029             /* We use orig_fill here because we want to fold the side, but if
1030             itheta==16384, we'll have cleared the low bits of fill. */
1031             cm = ff_celt_encode_band(f, rc, band, x2, NULL, N, mbits, blocks,
1032                                      lowband, duration, lowband_out, level, gain,
1033                                      lowband_scratch, orig_fill);
1034             /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
1035             and there's no need to worry about mixing with the other channel. */
1036             y2[0] = -sign * x2[1];
1037             y2[1] =  sign * x2[0];
1038             X[0] *= mid;
1039             X[1] *= mid;
1040             Y[0] *= side;
1041             Y[1] *= side;
1042             tmp = X[0];
1043             X[0] = tmp - Y[0];
1044             Y[0] = tmp + Y[0];
1045             tmp = X[1];
1046             X[1] = tmp - Y[1];
1047             Y[1] = tmp + Y[1];
1048         } else {
1049             /* "Normal" split code */
1050             float *next_lowband2     = NULL;
1051             float *next_lowband_out1 = NULL;
1052             int next_level = 0;
1053             int rebalance;
1054
1055             /* Give more bits to low-energy MDCTs than they would
1056              * otherwise deserve */
1057             if (B0 > 1 && !stereo && (itheta & 0x3fff)) {
1058                 if (itheta > 8192)
1059                     /* Rough approximation for pre-echo masking */
1060                     delta -= delta >> (4 - duration);
1061                 else
1062                     /* Corresponds to a forward-masking slope of
1063                      * 1.5 dB per 10 ms */
1064                     delta = FFMIN(0, delta + (N << 3 >> (5 - duration)));
1065             }
1066             mbits = av_clip((b - delta) / 2, 0, b);
1067             sbits = b - mbits;
1068             f->remaining2 -= qalloc;
1069
1070             if (lowband && !stereo)
1071                 next_lowband2 = lowband + N; /* >32-bit split case */
1072
1073             /* Only stereo needs to pass on lowband_out.
1074              * Otherwise, it's handled at the end */
1075             if (stereo)
1076                 next_lowband_out1 = lowband_out;
1077             else
1078                 next_level = level + 1;
1079
1080             rebalance = f->remaining2;
1081             if (mbits >= sbits) {
1082                 /* In stereo mode, we do not apply a scaling to the mid
1083                  * because we need the normalized mid for folding later */
1084                 cm = ff_celt_encode_band(f, rc, band, X, NULL, N, mbits, blocks,
1085                                          lowband, duration, next_lowband_out1,
1086                                          next_level, stereo ? 1.0f : (gain * mid),
1087                                          lowband_scratch, fill);
1088
1089                 rebalance = mbits - (rebalance - f->remaining2);
1090                 if (rebalance > 3 << 3 && itheta != 0)
1091                     sbits += rebalance - (3 << 3);
1092
1093                 /* For a stereo split, the high bits of fill are always zero,
1094                  * so no folding will be done to the side. */
1095                 cm |= ff_celt_encode_band(f, rc, band, Y, NULL, N, sbits, blocks,
1096                                           next_lowband2, duration, NULL,
1097                                           next_level, gain * side, NULL,
1098                                           fill >> blocks) << ((B0 >> 1) & (stereo - 1));
1099             } else {
1100                 /* For a stereo split, the high bits of fill are always zero,
1101                  * so no folding will be done to the side. */
1102                 cm = ff_celt_encode_band(f, rc, band, Y, NULL, N, sbits, blocks,
1103                                          next_lowband2, duration, NULL,
1104                                          next_level, gain * side, NULL,
1105                                          fill >> blocks) << ((B0 >> 1) & (stereo - 1));
1106
1107                 rebalance = sbits - (rebalance - f->remaining2);
1108                 if (rebalance > 3 << 3 && itheta != 16384)
1109                     mbits += rebalance - (3 << 3);
1110
1111                 /* In stereo mode, we do not apply a scaling to the mid because
1112                  * 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, stereo ? 1.0f : (gain * mid),
1116                                           lowband_scratch, fill);
1117             }
1118         }
1119     } else {
1120         /* This is the basic no-split case */
1121         uint32_t q         = celt_bits2pulses(cache, b);
1122         uint32_t curr_bits = celt_pulses2bits(cache, q);
1123         f->remaining2 -= curr_bits;
1124
1125         /* Ensures we can never bust the budget */
1126         while (f->remaining2 < 0 && q > 0) {
1127             f->remaining2 += curr_bits;
1128             curr_bits      = celt_pulses2bits(cache, --q);
1129             f->remaining2 -= curr_bits;
1130         }
1131
1132         if (q != 0) {
1133             /* Finally do the actual quantization */
1134             cm = celt_alg_quant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
1135                                 f->spread, blocks, gain);
1136         } else {
1137             /* If there's no pulse, fill the band anyway */
1138             uint32_t cm_mask = (1 << blocks) - 1;
1139             fill &= cm_mask;
1140             if (fill) {
1141                 if (!lowband) {
1142                     /* Noise */
1143                     for (i = 0; i < N; i++)
1144                         X[i] = (((int32_t)celt_rng(f)) >> 20);
1145                     cm = cm_mask;
1146                 } else {
1147                     /* Folded spectrum */
1148                     for (i = 0; i < N; i++) {
1149                         /* About 48 dB below the "normal" folding level */
1150                         X[i] = lowband[i] + (((celt_rng(f)) & 0x8000) ? 1.0f / 256 : -1.0f / 256);
1151                     }
1152                     cm = fill;
1153                 }
1154                 celt_renormalize_vector(X, N, gain);
1155             } else {
1156                 memset(X, 0, N*sizeof(float));
1157             }
1158         }
1159     }
1160
1161     /* This code is used by the decoder and by the resynthesis-enabled encoder */
1162     if (stereo) {
1163         if (N > 2)
1164             celt_stereo_merge(X, Y, mid, N);
1165         if (inv) {
1166             for (i = 0; i < N; i++)
1167                 Y[i] *= -1;
1168         }
1169     } else if (level == 0) {
1170         int k;
1171
1172         /* Undo the sample reorganization going from time order to frequency order */
1173         if (B0 > 1)
1174             celt_interleave_hadamard(f->scratch, X, N_B >> recombine,
1175                                      B0 << recombine, longblocks);
1176
1177         /* Undo time-freq changes that we did earlier */
1178         N_B = N_B0;
1179         blocks = B0;
1180         for (k = 0; k < time_divide; k++) {
1181             blocks >>= 1;
1182             N_B <<= 1;
1183             cm |= cm >> blocks;
1184             celt_haar1(X, N_B, blocks);
1185         }
1186
1187         for (k = 0; k < recombine; k++) {
1188             cm = ff_celt_bit_deinterleave[cm];
1189             celt_haar1(X, N0>>k, 1<<k);
1190         }
1191         blocks <<= recombine;
1192
1193         /* Scale output for later folding */
1194         if (lowband_out) {
1195             float n = sqrtf(N0);
1196             for (i = 0; i < N0; i++)
1197                 lowband_out[i] = n * X[i];
1198         }
1199         cm = av_mod_uintp2(cm, blocks);
1200     }
1201
1202     return cm;
1203 }
1204
1205 float ff_celt_quant_band_cost(CeltFrame *f, OpusRangeCoder *rc, int band, float *bits,
1206                               float lambda)
1207 {
1208     int i, b = 0;
1209     uint32_t cm[2] = { (1 << f->blocks) - 1, (1 << f->blocks) - 1 };
1210     const int band_size = ff_celt_freq_range[band] << f->size;
1211     float buf[352], lowband_scratch[176], norm1[176], norm2[176];
1212     float dist, cost, err_x = 0.0f, err_y = 0.0f;
1213     float *X = buf;
1214     float *X_orig = f->block[0].coeffs + (ff_celt_freq_bands[band] << f->size);
1215     float *Y = (f->channels == 2) ? &buf[176] : NULL;
1216     float *Y_orig = f->block[1].coeffs + (ff_celt_freq_bands[band] << f->size);
1217     OPUS_RC_CHECKPOINT_SPAWN(rc);
1218
1219     memcpy(X, X_orig, band_size*sizeof(float));
1220     if (Y)
1221         memcpy(Y, Y_orig, band_size*sizeof(float));
1222
1223     f->remaining2 = ((f->framebits << 3) - f->anticollapse_needed) - opus_rc_tell_frac(rc) - 1;
1224     if (band <= f->coded_bands - 1) {
1225         int curr_balance = f->remaining / FFMIN(3, f->coded_bands - band);
1226         b = av_clip_uintp2(FFMIN(f->remaining2 + 1, f->pulses[band] + curr_balance), 14);
1227     }
1228
1229     if (f->dual_stereo) {
1230         ff_celt_encode_band(f, rc, band, X, NULL, band_size, b / 2, f->blocks, NULL,
1231                             f->size, norm1, 0, 1.0f, lowband_scratch, cm[0]);
1232
1233         ff_celt_encode_band(f, rc, band, Y, NULL, band_size, b / 2, f->blocks, NULL,
1234                             f->size, norm2, 0, 1.0f, lowband_scratch, cm[1]);
1235     } else {
1236         ff_celt_encode_band(f, rc, band, X, Y, band_size, b, f->blocks, NULL, f->size,
1237                             norm1, 0, 1.0f, lowband_scratch, cm[0] | cm[1]);
1238     }
1239
1240     for (i = 0; i < band_size; i++) {
1241         err_x += (X[i] - X_orig[i])*(X[i] - X_orig[i]);
1242         err_y += (Y[i] - Y_orig[i])*(Y[i] - Y_orig[i]);
1243     }
1244
1245     dist = sqrtf(err_x) + sqrtf(err_y);
1246     cost = OPUS_RC_CHECKPOINT_BITS(rc)/8.0f;
1247     *bits += cost;
1248
1249     OPUS_RC_CHECKPOINT_ROLLBACK(rc);
1250
1251     return lambda*dist*cost;
1252 }