Save more integer divisions on ARM when we know the operands are positive
[opus.git] / celt / bands.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Copyright (c) 2008-2009 Gregory Maxwell
4    Written by Jean-Marc Valin and Gregory Maxwell */
5 /*
6    Redistribution and use in source and binary forms, with or without
7    modification, are permitted provided that the following conditions
8    are met:
9
10    - Redistributions of source code must retain the above copyright
11    notice, this list of conditions and the following disclaimer.
12
13    - Redistributions in binary form must reproduce the above copyright
14    notice, this list of conditions and the following disclaimer in the
15    documentation and/or other materials provided with the distribution.
16
17    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
18    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
19    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
20    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
21    OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
22    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
23    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
24    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
25    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
26    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
27    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30 #ifdef HAVE_CONFIG_H
31 #include "config.h"
32 #endif
33
34 #include <math.h>
35 #include "bands.h"
36 #include "modes.h"
37 #include "vq.h"
38 #include "cwrs.h"
39 #include "stack_alloc.h"
40 #include "os_support.h"
41 #include "mathops.h"
42 #include "rate.h"
43 #include "quant_bands.h"
44 #include "pitch.h"
45
46 int hysteresis_decision(opus_val16 val, const opus_val16 *thresholds, const opus_val16 *hysteresis, int N, int prev)
47 {
48    int i;
49    for (i=0;i<N;i++)
50    {
51       if (val < thresholds[i])
52          break;
53    }
54    if (i>prev && val < thresholds[prev]+hysteresis[prev])
55       i=prev;
56    if (i<prev && val > thresholds[prev-1]-hysteresis[prev-1])
57       i=prev;
58    return i;
59 }
60
61 opus_uint32 celt_lcg_rand(opus_uint32 seed)
62 {
63    return 1664525 * seed + 1013904223;
64 }
65
66 /* This is a cos() approximation designed to be bit-exact on any platform. Bit exactness
67    with this approximation is important because it has an impact on the bit allocation */
68 static opus_int16 bitexact_cos(opus_int16 x)
69 {
70    opus_int32 tmp;
71    opus_int16 x2;
72    tmp = (4096+((opus_int32)(x)*(x)))>>13;
73    celt_assert(tmp<=32767);
74    x2 = tmp;
75    x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));
76    celt_assert(x2<=32766);
77    return 1+x2;
78 }
79
80 static int bitexact_log2tan(int isin,int icos)
81 {
82    int lc;
83    int ls;
84    lc=EC_ILOG(icos);
85    ls=EC_ILOG(isin);
86    icos<<=15-lc;
87    isin<<=15-ls;
88    return (ls-lc)*(1<<11)
89          +FRAC_MUL16(isin, FRAC_MUL16(isin, -2597) + 7932)
90          -FRAC_MUL16(icos, FRAC_MUL16(icos, -2597) + 7932);
91 }
92
93 #ifdef FIXED_POINT
94 /* Compute the amplitude (sqrt energy) in each of the bands */
95 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int LM)
96 {
97    int i, c, N;
98    const opus_int16 *eBands = m->eBands;
99    N = m->shortMdctSize<<LM;
100    c=0; do {
101       for (i=0;i<end;i++)
102       {
103          int j;
104          opus_val32 maxval=0;
105          opus_val32 sum = 0;
106
107          maxval = celt_maxabs32(&X[c*N+(eBands[i]<<LM)], (eBands[i+1]-eBands[i])<<LM);
108          if (maxval > 0)
109          {
110             int shift = celt_ilog2(maxval) - 14 + (((m->logN[i]>>BITRES)+LM+1)>>1);
111             j=eBands[i]<<LM;
112             if (shift>0)
113             {
114                do {
115                   sum = MAC16_16(sum, EXTRACT16(SHR32(X[j+c*N],shift)),
116                         EXTRACT16(SHR32(X[j+c*N],shift)));
117                } while (++j<eBands[i+1]<<LM);
118             } else {
119                do {
120                   sum = MAC16_16(sum, EXTRACT16(SHL32(X[j+c*N],-shift)),
121                         EXTRACT16(SHL32(X[j+c*N],-shift)));
122                } while (++j<eBands[i+1]<<LM);
123             }
124             /* We're adding one here to ensure the normalized band isn't larger than unity norm */
125             bandE[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
126          } else {
127             bandE[i+c*m->nbEBands] = EPSILON;
128          }
129          /*printf ("%f ", bandE[i+c*m->nbEBands]);*/
130       }
131    } while (++c<C);
132    /*printf ("\n");*/
133 }
134
135 /* Normalise each band such that the energy is one. */
136 void normalise_bands(const CELTMode *m, const celt_sig * OPUS_RESTRICT freq, celt_norm * OPUS_RESTRICT X, const celt_ener *bandE, int end, int C, int M)
137 {
138    int i, c, N;
139    const opus_int16 *eBands = m->eBands;
140    N = M*m->shortMdctSize;
141    c=0; do {
142       i=0; do {
143          opus_val16 g;
144          int j,shift;
145          opus_val16 E;
146          shift = celt_zlog2(bandE[i+c*m->nbEBands])-13;
147          E = VSHR32(bandE[i+c*m->nbEBands], shift);
148          g = EXTRACT16(celt_rcp(SHL32(E,3)));
149          j=M*eBands[i]; do {
150             X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
151          } while (++j<M*eBands[i+1]);
152       } while (++i<end);
153    } while (++c<C);
154 }
155
156 #else /* FIXED_POINT */
157 /* Compute the amplitude (sqrt energy) in each of the bands */
158 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int LM)
159 {
160    int i, c, N;
161    const opus_int16 *eBands = m->eBands;
162    N = m->shortMdctSize<<LM;
163    c=0; do {
164       for (i=0;i<end;i++)
165       {
166          opus_val32 sum;
167          sum = 1e-27f + celt_inner_prod(&X[c*N+(eBands[i]<<LM)], &X[c*N+(eBands[i]<<LM)], (eBands[i+1]-eBands[i])<<LM);
168          bandE[i+c*m->nbEBands] = celt_sqrt(sum);
169          /*printf ("%f ", bandE[i+c*m->nbEBands]);*/
170       }
171    } while (++c<C);
172    /*printf ("\n");*/
173 }
174
175 /* Normalise each band such that the energy is one. */
176 void normalise_bands(const CELTMode *m, const celt_sig * OPUS_RESTRICT freq, celt_norm * OPUS_RESTRICT X, const celt_ener *bandE, int end, int C, int M)
177 {
178    int i, c, N;
179    const opus_int16 *eBands = m->eBands;
180    N = M*m->shortMdctSize;
181    c=0; do {
182       for (i=0;i<end;i++)
183       {
184          int j;
185          opus_val16 g = 1.f/(1e-27f+bandE[i+c*m->nbEBands]);
186          for (j=M*eBands[i];j<M*eBands[i+1];j++)
187             X[j+c*N] = freq[j+c*N]*g;
188       }
189    } while (++c<C);
190 }
191
192 #endif /* FIXED_POINT */
193
194 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
195 void denormalise_bands(const CELTMode *m, const celt_norm * OPUS_RESTRICT X,
196       celt_sig * OPUS_RESTRICT freq, const opus_val16 *bandLogE, int start,
197       int end, int M, int downsample, int silence)
198 {
199    int i, N;
200    int bound;
201    celt_sig * OPUS_RESTRICT f;
202    const celt_norm * OPUS_RESTRICT x;
203    const opus_int16 *eBands = m->eBands;
204    N = M*m->shortMdctSize;
205    bound = M*eBands[end];
206    if (downsample!=1)
207       bound = IMIN(bound, N/downsample);
208    if (silence)
209    {
210       bound = 0;
211       start = end = 0;
212    }
213    f = freq;
214    x = X+M*eBands[start];
215    for (i=0;i<M*eBands[start];i++)
216       *f++ = 0;
217    for (i=start;i<end;i++)
218    {
219       int j, band_end;
220       opus_val16 g;
221       opus_val16 lg;
222 #ifdef FIXED_POINT
223       int shift;
224 #endif
225       j=M*eBands[i];
226       band_end = M*eBands[i+1];
227       lg = ADD16(bandLogE[i], SHL16((opus_val16)eMeans[i],6));
228 #ifndef FIXED_POINT
229       g = celt_exp2(lg);
230 #else
231       /* Handle the integer part of the log energy */
232       shift = 16-(lg>>DB_SHIFT);
233       if (shift>31)
234       {
235          shift=0;
236          g=0;
237       } else {
238          /* Handle the fractional part. */
239          g = celt_exp2_frac(lg&((1<<DB_SHIFT)-1));
240       }
241       /* Handle extreme gains with negative shift. */
242       if (shift<0)
243       {
244          /* For shift < -2 we'd be likely to overflow, so we're capping
245                the gain here. This shouldn't happen unless the bitstream is
246                already corrupted. */
247          if (shift < -2)
248          {
249             g = 32767;
250             shift = -2;
251          }
252          do {
253             *f++ = SHL32(MULT16_16(*x++, g), -shift);
254          } while (++j<band_end);
255       } else
256 #endif
257          /* Be careful of the fixed-point "else" just above when changing this code */
258          do {
259             *f++ = SHR32(MULT16_16(*x++, g), shift);
260          } while (++j<band_end);
261    }
262    celt_assert(start <= end);
263    OPUS_CLEAR(&freq[bound], N-bound);
264 }
265
266 /* This prevents energy collapse for transients with multiple short MDCTs */
267 void anti_collapse(const CELTMode *m, celt_norm *X_, unsigned char *collapse_masks, int LM, int C, int size,
268       int start, int end, const opus_val16 *logE, const opus_val16 *prev1logE,
269       const opus_val16 *prev2logE, const int *pulses, opus_uint32 seed)
270 {
271    int c, i, j, k;
272    for (i=start;i<end;i++)
273    {
274       int N0;
275       opus_val16 thresh, sqrt_1;
276       int depth;
277 #ifdef FIXED_POINT
278       int shift;
279       opus_val32 thresh32;
280 #endif
281
282       N0 = m->eBands[i+1]-m->eBands[i];
283       /* depth in 1/8 bits */
284       celt_assert(pulses[i]>=0);
285       depth = celt_udiv(1+pulses[i], (m->eBands[i+1]-m->eBands[i])<<LM);
286
287 #ifdef FIXED_POINT
288       thresh32 = SHR32(celt_exp2(-SHL16(depth, 10-BITRES)),1);
289       thresh = MULT16_32_Q15(QCONST16(0.5f, 15), MIN32(32767,thresh32));
290       {
291          opus_val32 t;
292          t = N0<<LM;
293          shift = celt_ilog2(t)>>1;
294          t = SHL32(t, (7-shift)<<1);
295          sqrt_1 = celt_rsqrt_norm(t);
296       }
297 #else
298       thresh = .5f*celt_exp2(-.125f*depth);
299       sqrt_1 = celt_rsqrt(N0<<LM);
300 #endif
301
302       c=0; do
303       {
304          celt_norm *X;
305          opus_val16 prev1;
306          opus_val16 prev2;
307          opus_val32 Ediff;
308          opus_val16 r;
309          int renormalize=0;
310          prev1 = prev1logE[c*m->nbEBands+i];
311          prev2 = prev2logE[c*m->nbEBands+i];
312          if (C==1)
313          {
314             prev1 = MAX16(prev1,prev1logE[m->nbEBands+i]);
315             prev2 = MAX16(prev2,prev2logE[m->nbEBands+i]);
316          }
317          Ediff = EXTEND32(logE[c*m->nbEBands+i])-EXTEND32(MIN16(prev1,prev2));
318          Ediff = MAX32(0, Ediff);
319
320 #ifdef FIXED_POINT
321          if (Ediff < 16384)
322          {
323             opus_val32 r32 = SHR32(celt_exp2(-EXTRACT16(Ediff)),1);
324             r = 2*MIN16(16383,r32);
325          } else {
326             r = 0;
327          }
328          if (LM==3)
329             r = MULT16_16_Q14(23170, MIN32(23169, r));
330          r = SHR16(MIN16(thresh, r),1);
331          r = SHR32(MULT16_16_Q15(sqrt_1, r),shift);
332 #else
333          /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because
334             short blocks don't have the same energy as long */
335          r = 2.f*celt_exp2(-Ediff);
336          if (LM==3)
337             r *= 1.41421356f;
338          r = MIN16(thresh, r);
339          r = r*sqrt_1;
340 #endif
341          X = X_+c*size+(m->eBands[i]<<LM);
342          for (k=0;k<1<<LM;k++)
343          {
344             /* Detect collapse */
345             if (!(collapse_masks[i*C+c]&1<<k))
346             {
347                /* Fill with noise */
348                for (j=0;j<N0;j++)
349                {
350                   seed = celt_lcg_rand(seed);
351                   X[(j<<LM)+k] = (seed&0x8000 ? r : -r);
352                }
353                renormalize = 1;
354             }
355          }
356          /* We just added some energy, so we need to renormalise */
357          if (renormalize)
358             renormalise_vector(X, N0<<LM, Q15ONE);
359       } while (++c<C);
360    }
361 }
362
363 static void intensity_stereo(const CELTMode *m, celt_norm * OPUS_RESTRICT X, const celt_norm * OPUS_RESTRICT Y, const celt_ener *bandE, int bandID, int N)
364 {
365    int i = bandID;
366    int j;
367    opus_val16 a1, a2;
368    opus_val16 left, right;
369    opus_val16 norm;
370 #ifdef FIXED_POINT
371    int shift = celt_zlog2(MAX32(bandE[i], bandE[i+m->nbEBands]))-13;
372 #endif
373    left = VSHR32(bandE[i],shift);
374    right = VSHR32(bandE[i+m->nbEBands],shift);
375    norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
376    a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
377    a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);
378    for (j=0;j<N;j++)
379    {
380       celt_norm r, l;
381       l = X[j];
382       r = Y[j];
383       X[j] = EXTRACT16(SHR32(MAC16_16(MULT16_16(a1, l), a2, r), 14));
384       /* Side is not encoded, no need to calculate */
385    }
386 }
387
388 static void stereo_split(celt_norm * OPUS_RESTRICT X, celt_norm * OPUS_RESTRICT Y, int N)
389 {
390    int j;
391    for (j=0;j<N;j++)
392    {
393       opus_val32 r, l;
394       l = MULT16_16(QCONST16(.70710678f, 15), X[j]);
395       r = MULT16_16(QCONST16(.70710678f, 15), Y[j]);
396       X[j] = EXTRACT16(SHR32(ADD32(l, r), 15));
397       Y[j] = EXTRACT16(SHR32(SUB32(r, l), 15));
398    }
399 }
400
401 static void stereo_merge(celt_norm * OPUS_RESTRICT X, celt_norm * OPUS_RESTRICT Y, opus_val16 mid, int N)
402 {
403    int j;
404    opus_val32 xp=0, side=0;
405    opus_val32 El, Er;
406    opus_val16 mid2;
407 #ifdef FIXED_POINT
408    int kl, kr;
409 #endif
410    opus_val32 t, lgain, rgain;
411
412    /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
413    dual_inner_prod(Y, X, Y, N, &xp, &side);
414    /* Compensating for the mid normalization */
415    xp = MULT16_32_Q15(mid, xp);
416    /* mid and side are in Q15, not Q14 like X and Y */
417    mid2 = SHR32(mid, 1);
418    El = MULT16_16(mid2, mid2) + side - 2*xp;
419    Er = MULT16_16(mid2, mid2) + side + 2*xp;
420    if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28))
421    {
422       OPUS_COPY(Y, X, N);
423       return;
424    }
425
426 #ifdef FIXED_POINT
427    kl = celt_ilog2(El)>>1;
428    kr = celt_ilog2(Er)>>1;
429 #endif
430    t = VSHR32(El, (kl-7)<<1);
431    lgain = celt_rsqrt_norm(t);
432    t = VSHR32(Er, (kr-7)<<1);
433    rgain = celt_rsqrt_norm(t);
434
435 #ifdef FIXED_POINT
436    if (kl < 7)
437       kl = 7;
438    if (kr < 7)
439       kr = 7;
440 #endif
441
442    for (j=0;j<N;j++)
443    {
444       celt_norm r, l;
445       /* Apply mid scaling (side is already scaled) */
446       l = MULT16_16_P15(mid, X[j]);
447       r = Y[j];
448       X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));
449       Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));
450    }
451 }
452
453 /* Decide whether we should spread the pulses in the current frame */
454 int spreading_decision(const CELTMode *m, const celt_norm *X, int *average,
455       int last_decision, int *hf_average, int *tapset_decision, int update_hf,
456       int end, int C, int M)
457 {
458    int i, c, N0;
459    int sum = 0, nbBands=0;
460    const opus_int16 * OPUS_RESTRICT eBands = m->eBands;
461    int decision;
462    int hf_sum=0;
463
464    celt_assert(end>0);
465
466    N0 = M*m->shortMdctSize;
467
468    if (M*(eBands[end]-eBands[end-1]) <= 8)
469       return SPREAD_NONE;
470    c=0; do {
471       for (i=0;i<end;i++)
472       {
473          int j, N, tmp=0;
474          int tcount[3] = {0,0,0};
475          const celt_norm * OPUS_RESTRICT x = X+M*eBands[i]+c*N0;
476          N = M*(eBands[i+1]-eBands[i]);
477          if (N<=8)
478             continue;
479          /* Compute rough CDF of |x[j]| */
480          for (j=0;j<N;j++)
481          {
482             opus_val32 x2N; /* Q13 */
483
484             x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
485             if (x2N < QCONST16(0.25f,13))
486                tcount[0]++;
487             if (x2N < QCONST16(0.0625f,13))
488                tcount[1]++;
489             if (x2N < QCONST16(0.015625f,13))
490                tcount[2]++;
491          }
492
493          /* Only include four last bands (8 kHz and up) */
494          if (i>m->nbEBands-4)
495             hf_sum += celt_udiv(32*(tcount[1]+tcount[0]), N);
496          tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
497          sum += tmp*256;
498          nbBands++;
499       }
500    } while (++c<C);
501
502    if (update_hf)
503    {
504       if (hf_sum)
505          hf_sum = celt_udiv(hf_sum, C*(4-m->nbEBands+end));
506       *hf_average = (*hf_average+hf_sum)>>1;
507       hf_sum = *hf_average;
508       if (*tapset_decision==2)
509          hf_sum += 4;
510       else if (*tapset_decision==0)
511          hf_sum -= 4;
512       if (hf_sum > 22)
513          *tapset_decision=2;
514       else if (hf_sum > 18)
515          *tapset_decision=1;
516       else
517          *tapset_decision=0;
518    }
519    /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/
520    celt_assert(nbBands>0); /* end has to be non-zero */
521    celt_assert(sum>=0);
522    sum = celt_udiv(sum, nbBands);
523    /* Recursive averaging */
524    sum = (sum+*average)>>1;
525    *average = sum;
526    /* Hysteresis */
527    sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2;
528    if (sum < 80)
529    {
530       decision = SPREAD_AGGRESSIVE;
531    } else if (sum < 256)
532    {
533       decision = SPREAD_NORMAL;
534    } else if (sum < 384)
535    {
536       decision = SPREAD_LIGHT;
537    } else {
538       decision = SPREAD_NONE;
539    }
540 #ifdef FUZZING
541    decision = rand()&0x3;
542    *tapset_decision=rand()%3;
543 #endif
544    return decision;
545 }
546
547 /* Indexing table for converting from natural Hadamard to ordery Hadamard
548    This is essentially a bit-reversed Gray, on top of which we've added
549    an inversion of the order because we want the DC at the end rather than
550    the beginning. The lines are for N=2, 4, 8, 16 */
551 static const int ordery_table[] = {
552        1,  0,
553        3,  0,  2,  1,
554        7,  0,  4,  3,  6,  1,  5,  2,
555       15,  0,  8,  7, 12,  3, 11,  4, 14,  1,  9,  6, 13,  2, 10,  5,
556 };
557
558 static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
559 {
560    int i,j;
561    VARDECL(celt_norm, tmp);
562    int N;
563    SAVE_STACK;
564    N = N0*stride;
565    ALLOC(tmp, N, celt_norm);
566    celt_assert(stride>0);
567    if (hadamard)
568    {
569       const int *ordery = ordery_table+stride-2;
570       for (i=0;i<stride;i++)
571       {
572          for (j=0;j<N0;j++)
573             tmp[ordery[i]*N0+j] = X[j*stride+i];
574       }
575    } else {
576       for (i=0;i<stride;i++)
577          for (j=0;j<N0;j++)
578             tmp[i*N0+j] = X[j*stride+i];
579    }
580    OPUS_COPY(X, tmp, N);
581    RESTORE_STACK;
582 }
583
584 static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
585 {
586    int i,j;
587    VARDECL(celt_norm, tmp);
588    int N;
589    SAVE_STACK;
590    N = N0*stride;
591    ALLOC(tmp, N, celt_norm);
592    if (hadamard)
593    {
594       const int *ordery = ordery_table+stride-2;
595       for (i=0;i<stride;i++)
596          for (j=0;j<N0;j++)
597             tmp[j*stride+i] = X[ordery[i]*N0+j];
598    } else {
599       for (i=0;i<stride;i++)
600          for (j=0;j<N0;j++)
601             tmp[j*stride+i] = X[i*N0+j];
602    }
603    OPUS_COPY(X, tmp, N);
604    RESTORE_STACK;
605 }
606
607 void haar1(celt_norm *X, int N0, int stride)
608 {
609    int i, j;
610    N0 >>= 1;
611    for (i=0;i<stride;i++)
612       for (j=0;j<N0;j++)
613       {
614          opus_val32 tmp1, tmp2;
615          tmp1 = MULT16_16(QCONST16(.70710678f,15), X[stride*2*j+i]);
616          tmp2 = MULT16_16(QCONST16(.70710678f,15), X[stride*(2*j+1)+i]);
617          X[stride*2*j+i] = EXTRACT16(PSHR32(ADD32(tmp1, tmp2), 15));
618          X[stride*(2*j+1)+i] = EXTRACT16(PSHR32(SUB32(tmp1, tmp2), 15));
619       }
620 }
621
622 static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo)
623 {
624    static const opus_int16 exp2_table8[8] =
625       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
626    int qn, qb;
627    int N2 = 2*N-1;
628    if (stereo && N==2)
629       N2--;
630    /* The upper limit ensures that in a stereo split with itheta==16384, we'll
631        always have enough bits left over to code at least one pulse in the
632        side; otherwise it would collapse, since it doesn't get folded. */
633    qb = IMIN(b-pulse_cap-(4<<BITRES), (b+N2*offset)/N2);
634
635    qb = IMIN(8<<BITRES, qb);
636
637    if (qb<(1<<BITRES>>1)) {
638       qn = 1;
639    } else {
640       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
641       qn = (qn+1)>>1<<1;
642    }
643    celt_assert(qn <= 256);
644    return qn;
645 }
646
647 struct band_ctx {
648    int encode;
649    const CELTMode *m;
650    int i;
651    int intensity;
652    int spread;
653    int tf_change;
654    ec_ctx *ec;
655    opus_int32 remaining_bits;
656    const celt_ener *bandE;
657    opus_uint32 seed;
658 };
659
660 struct split_ctx {
661    int inv;
662    int imid;
663    int iside;
664    int delta;
665    int itheta;
666    int qalloc;
667 };
668
669 static void compute_theta(struct band_ctx *ctx, struct split_ctx *sctx,
670       celt_norm *X, celt_norm *Y, int N, int *b, int B, int B0,
671       int LM,
672       int stereo, int *fill)
673 {
674    int qn;
675    int itheta=0;
676    int delta;
677    int imid, iside;
678    int qalloc;
679    int pulse_cap;
680    int offset;
681    opus_int32 tell;
682    int inv=0;
683    int encode;
684    const CELTMode *m;
685    int i;
686    int intensity;
687    ec_ctx *ec;
688    const celt_ener *bandE;
689
690    encode = ctx->encode;
691    m = ctx->m;
692    i = ctx->i;
693    intensity = ctx->intensity;
694    ec = ctx->ec;
695    bandE = ctx->bandE;
696
697    /* Decide on the resolution to give to the split parameter theta */
698    pulse_cap = m->logN[i]+LM*(1<<BITRES);
699    offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
700    qn = compute_qn(N, *b, offset, pulse_cap, stereo);
701    if (stereo && i>=intensity)
702       qn = 1;
703    if (encode)
704    {
705       /* theta is the atan() of the ratio between the (normalized)
706          side and mid. With just that parameter, we can re-scale both
707          mid and side because we know that 1) they have unit norm and
708          2) they are orthogonal. */
709       itheta = stereo_itheta(X, Y, stereo, N);
710    }
711    tell = ec_tell_frac(ec);
712    if (qn!=1)
713    {
714       if (encode)
715          itheta = (itheta*qn+8192)>>14;
716
717       /* Entropy coding of the angle. We use a uniform pdf for the
718          time split, a step for stereo, and a triangular one for the rest. */
719       if (stereo && N>2)
720       {
721          int p0 = 3;
722          int x = itheta;
723          int x0 = qn/2;
724          int ft = p0*(x0+1) + x0;
725          /* Use a probability of p0 up to itheta=8192 and then use 1 after */
726          if (encode)
727          {
728             ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
729          } else {
730             int fs;
731             fs=ec_decode(ec,ft);
732             if (fs<(x0+1)*p0)
733                x=fs/p0;
734             else
735                x=x0+1+(fs-(x0+1)*p0);
736             ec_dec_update(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
737             itheta = x;
738          }
739       } else if (B0>1 || stereo) {
740          /* Uniform pdf */
741          if (encode)
742             ec_enc_uint(ec, itheta, qn+1);
743          else
744             itheta = ec_dec_uint(ec, qn+1);
745       } else {
746          int fs=1, ft;
747          ft = ((qn>>1)+1)*((qn>>1)+1);
748          if (encode)
749          {
750             int fl;
751
752             fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
753             fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
754              ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
755
756             ec_encode(ec, fl, fl+fs, ft);
757          } else {
758             /* Triangular pdf */
759             int fl=0;
760             int fm;
761             fm = ec_decode(ec, ft);
762
763             if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
764             {
765                itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;
766                fs = itheta + 1;
767                fl = itheta*(itheta + 1)>>1;
768             }
769             else
770             {
771                itheta = (2*(qn + 1)
772                 - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;
773                fs = qn + 1 - itheta;
774                fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
775             }
776
777             ec_dec_update(ec, fl, fl+fs, ft);
778          }
779       }
780       celt_assert(itheta>=0);
781       itheta = celt_udiv((opus_int32)itheta*16384, qn);
782       if (encode && stereo)
783       {
784          if (itheta==0)
785             intensity_stereo(m, X, Y, bandE, i, N);
786          else
787             stereo_split(X, Y, N);
788       }
789       /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.
790                Let's do that at higher complexity */
791    } else if (stereo) {
792       if (encode)
793       {
794          inv = itheta > 8192;
795          if (inv)
796          {
797             int j;
798             for (j=0;j<N;j++)
799                Y[j] = -Y[j];
800          }
801          intensity_stereo(m, X, Y, bandE, i, N);
802       }
803       if (*b>2<<BITRES && ctx->remaining_bits > 2<<BITRES)
804       {
805          if (encode)
806             ec_enc_bit_logp(ec, inv, 2);
807          else
808             inv = ec_dec_bit_logp(ec, 2);
809       } else
810          inv = 0;
811       itheta = 0;
812    }
813    qalloc = ec_tell_frac(ec) - tell;
814    *b -= qalloc;
815
816    if (itheta == 0)
817    {
818       imid = 32767;
819       iside = 0;
820       *fill &= (1<<B)-1;
821       delta = -16384;
822    } else if (itheta == 16384)
823    {
824       imid = 0;
825       iside = 32767;
826       *fill &= ((1<<B)-1)<<B;
827       delta = 16384;
828    } else {
829       imid = bitexact_cos((opus_int16)itheta);
830       iside = bitexact_cos((opus_int16)(16384-itheta));
831       /* This is the mid vs side allocation that minimizes squared error
832          in that band. */
833       delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid));
834    }
835
836    sctx->inv = inv;
837    sctx->imid = imid;
838    sctx->iside = iside;
839    sctx->delta = delta;
840    sctx->itheta = itheta;
841    sctx->qalloc = qalloc;
842 }
843 static unsigned quant_band_n1(struct band_ctx *ctx, celt_norm *X, celt_norm *Y, int b,
844       celt_norm *lowband_out)
845 {
846 #ifdef RESYNTH
847    int resynth = 1;
848 #else
849    int resynth = !ctx->encode;
850 #endif
851    int c;
852    int stereo;
853    celt_norm *x = X;
854    int encode;
855    ec_ctx *ec;
856
857    encode = ctx->encode;
858    ec = ctx->ec;
859
860    stereo = Y != NULL;
861    c=0; do {
862       int sign=0;
863       if (ctx->remaining_bits>=1<<BITRES)
864       {
865          if (encode)
866          {
867             sign = x[0]<0;
868             ec_enc_bits(ec, sign, 1);
869          } else {
870             sign = ec_dec_bits(ec, 1);
871          }
872          ctx->remaining_bits -= 1<<BITRES;
873          b-=1<<BITRES;
874       }
875       if (resynth)
876          x[0] = sign ? -NORM_SCALING : NORM_SCALING;
877       x = Y;
878    } while (++c<1+stereo);
879    if (lowband_out)
880       lowband_out[0] = SHR16(X[0],4);
881    return 1;
882 }
883
884 /* This function is responsible for encoding and decoding a mono partition.
885    It can split the band in two and transmit the energy difference with
886    the two half-bands. It can be called recursively so bands can end up being
887    split in 8 parts. */
888 static unsigned quant_partition(struct band_ctx *ctx, celt_norm *X,
889       int N, int b, int B, celt_norm *lowband,
890       int LM,
891       opus_val16 gain, int fill)
892 {
893    const unsigned char *cache;
894    int q;
895    int curr_bits;
896    int imid=0, iside=0;
897    int B0=B;
898    opus_val16 mid=0, side=0;
899    unsigned cm=0;
900 #ifdef RESYNTH
901    int resynth = 1;
902 #else
903    int resynth = !ctx->encode;
904 #endif
905    celt_norm *Y=NULL;
906    int encode;
907    const CELTMode *m;
908    int i;
909    int spread;
910    ec_ctx *ec;
911
912    encode = ctx->encode;
913    m = ctx->m;
914    i = ctx->i;
915    spread = ctx->spread;
916    ec = ctx->ec;
917
918    /* If we need 1.5 more bit than we can produce, split the band in two. */
919    cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
920    if (LM != -1 && b > cache[cache[0]]+12 && N>2)
921    {
922       int mbits, sbits, delta;
923       int itheta;
924       int qalloc;
925       struct split_ctx sctx;
926       celt_norm *next_lowband2=NULL;
927       opus_int32 rebalance;
928
929       N >>= 1;
930       Y = X+N;
931       LM -= 1;
932       if (B==1)
933          fill = (fill&1)|(fill<<1);
934       B = (B+1)>>1;
935
936       compute_theta(ctx, &sctx, X, Y, N, &b, B, B0,
937             LM, 0, &fill);
938       imid = sctx.imid;
939       iside = sctx.iside;
940       delta = sctx.delta;
941       itheta = sctx.itheta;
942       qalloc = sctx.qalloc;
943 #ifdef FIXED_POINT
944       mid = imid;
945       side = iside;
946 #else
947       mid = (1.f/32768)*imid;
948       side = (1.f/32768)*iside;
949 #endif
950
951       /* Give more bits to low-energy MDCTs than they would otherwise deserve */
952       if (B0>1 && (itheta&0x3fff))
953       {
954          if (itheta > 8192)
955             /* Rough approximation for pre-echo masking */
956             delta -= delta>>(4-LM);
957          else
958             /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
959             delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
960       }
961       mbits = IMAX(0, IMIN(b, (b-delta)/2));
962       sbits = b-mbits;
963       ctx->remaining_bits -= qalloc;
964
965       if (lowband)
966          next_lowband2 = lowband+N; /* >32-bit split case */
967
968       rebalance = ctx->remaining_bits;
969       if (mbits >= sbits)
970       {
971          cm = quant_partition(ctx, X, N, mbits, B,
972                lowband, LM,
973                MULT16_16_P15(gain,mid), fill);
974          rebalance = mbits - (rebalance-ctx->remaining_bits);
975          if (rebalance > 3<<BITRES && itheta!=0)
976             sbits += rebalance - (3<<BITRES);
977          cm |= quant_partition(ctx, Y, N, sbits, B,
978                next_lowband2, LM,
979                MULT16_16_P15(gain,side), fill>>B)<<(B0>>1);
980       } else {
981          cm = quant_partition(ctx, Y, N, sbits, B,
982                next_lowband2, LM,
983                MULT16_16_P15(gain,side), fill>>B)<<(B0>>1);
984          rebalance = sbits - (rebalance-ctx->remaining_bits);
985          if (rebalance > 3<<BITRES && itheta!=16384)
986             mbits += rebalance - (3<<BITRES);
987          cm |= quant_partition(ctx, X, N, mbits, B,
988                lowband, LM,
989                MULT16_16_P15(gain,mid), fill);
990       }
991    } else {
992       /* This is the basic no-split case */
993       q = bits2pulses(m, i, LM, b);
994       curr_bits = pulses2bits(m, i, LM, q);
995       ctx->remaining_bits -= curr_bits;
996
997       /* Ensures we can never bust the budget */
998       while (ctx->remaining_bits < 0 && q > 0)
999       {
1000          ctx->remaining_bits += curr_bits;
1001          q--;
1002          curr_bits = pulses2bits(m, i, LM, q);
1003          ctx->remaining_bits -= curr_bits;
1004       }
1005
1006       if (q!=0)
1007       {
1008          int K = get_pulses(q);
1009
1010          /* Finally do the actual quantization */
1011          if (encode)
1012          {
1013             cm = alg_quant(X, N, K, spread, B, ec
1014 #ifdef RESYNTH
1015                  , gain
1016 #endif
1017                  );
1018          } else {
1019             cm = alg_unquant(X, N, K, spread, B, ec, gain);
1020          }
1021       } else {
1022          /* If there's no pulse, fill the band anyway */
1023          int j;
1024          if (resynth)
1025          {
1026             unsigned cm_mask;
1027             /* B can be as large as 16, so this shift might overflow an int on a
1028                16-bit platform; use a long to get defined behavior.*/
1029             cm_mask = (unsigned)(1UL<<B)-1;
1030             fill &= cm_mask;
1031             if (!fill)
1032             {
1033                OPUS_CLEAR(X, N);
1034             } else {
1035                if (lowband == NULL)
1036                {
1037                   /* Noise */
1038                   for (j=0;j<N;j++)
1039                   {
1040                      ctx->seed = celt_lcg_rand(ctx->seed);
1041                      X[j] = (celt_norm)((opus_int32)ctx->seed>>20);
1042                   }
1043                   cm = cm_mask;
1044                } else {
1045                   /* Folded spectrum */
1046                   for (j=0;j<N;j++)
1047                   {
1048                      opus_val16 tmp;
1049                      ctx->seed = celt_lcg_rand(ctx->seed);
1050                      /* About 48 dB below the "normal" folding level */
1051                      tmp = QCONST16(1.0f/256, 10);
1052                      tmp = (ctx->seed)&0x8000 ? tmp : -tmp;
1053                      X[j] = lowband[j]+tmp;
1054                   }
1055                   cm = fill;
1056                }
1057                renormalise_vector(X, N, gain);
1058             }
1059          }
1060       }
1061    }
1062
1063    return cm;
1064 }
1065
1066
1067 /* This function is responsible for encoding and decoding a band for the mono case. */
1068 static unsigned quant_band(struct band_ctx *ctx, celt_norm *X,
1069       int N, int b, int B, celt_norm *lowband,
1070       int LM, celt_norm *lowband_out,
1071       opus_val16 gain, celt_norm *lowband_scratch, int fill)
1072 {
1073    int N0=N;
1074    int N_B=N;
1075    int N_B0;
1076    int B0=B;
1077    int time_divide=0;
1078    int recombine=0;
1079    int longBlocks;
1080    unsigned cm=0;
1081 #ifdef RESYNTH
1082    int resynth = 1;
1083 #else
1084    int resynth = !ctx->encode;
1085 #endif
1086    int k;
1087    int encode;
1088    int tf_change;
1089
1090    encode = ctx->encode;
1091    tf_change = ctx->tf_change;
1092
1093    longBlocks = B0==1;
1094
1095    N_B = celt_udiv(N_B, B);
1096
1097    /* Special case for one sample */
1098    if (N==1)
1099    {
1100       return quant_band_n1(ctx, X, NULL, b, lowband_out);
1101    }
1102
1103    if (tf_change>0)
1104       recombine = tf_change;
1105    /* Band recombining to increase frequency resolution */
1106
1107    if (lowband_scratch && lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
1108    {
1109       OPUS_COPY(lowband_scratch, lowband, N);
1110       lowband = lowband_scratch;
1111    }
1112
1113    for (k=0;k<recombine;k++)
1114    {
1115       static const unsigned char bit_interleave_table[16]={
1116             0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
1117       };
1118       if (encode)
1119          haar1(X, N>>k, 1<<k);
1120       if (lowband)
1121          haar1(lowband, N>>k, 1<<k);
1122       fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
1123    }
1124    B>>=recombine;
1125    N_B<<=recombine;
1126
1127    /* Increasing the time resolution */
1128    while ((N_B&1) == 0 && tf_change<0)
1129    {
1130       if (encode)
1131          haar1(X, N_B, B);
1132       if (lowband)
1133          haar1(lowband, N_B, B);
1134       fill |= fill<<B;
1135       B <<= 1;
1136       N_B >>= 1;
1137       time_divide++;
1138       tf_change++;
1139    }
1140    B0=B;
1141    N_B0 = N_B;
1142
1143    /* Reorganize the samples in time order instead of frequency order */
1144    if (B0>1)
1145    {
1146       if (encode)
1147          deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1148       if (lowband)
1149          deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
1150    }
1151
1152    cm = quant_partition(ctx, X, N, b, B, lowband,
1153          LM, gain, fill);
1154
1155    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1156    if (resynth)
1157    {
1158       /* Undo the sample reorganization going from time order to frequency order */
1159       if (B0>1)
1160          interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1161
1162       /* Undo time-freq changes that we did earlier */
1163       N_B = N_B0;
1164       B = B0;
1165       for (k=0;k<time_divide;k++)
1166       {
1167          B >>= 1;
1168          N_B <<= 1;
1169          cm |= cm>>B;
1170          haar1(X, N_B, B);
1171       }
1172
1173       for (k=0;k<recombine;k++)
1174       {
1175          static const unsigned char bit_deinterleave_table[16]={
1176                0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1177                0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1178          };
1179          cm = bit_deinterleave_table[cm];
1180          haar1(X, N0>>k, 1<<k);
1181       }
1182       B<<=recombine;
1183
1184       /* Scale output for later folding */
1185       if (lowband_out)
1186       {
1187          int j;
1188          opus_val16 n;
1189          n = celt_sqrt(SHL32(EXTEND32(N0),22));
1190          for (j=0;j<N0;j++)
1191             lowband_out[j] = MULT16_16_Q15(n,X[j]);
1192       }
1193       cm &= (1<<B)-1;
1194    }
1195    return cm;
1196 }
1197
1198
1199 /* This function is responsible for encoding and decoding a band for the stereo case. */
1200 static unsigned quant_band_stereo(struct band_ctx *ctx, celt_norm *X, celt_norm *Y,
1201       int N, int b, int B, celt_norm *lowband,
1202       int LM, celt_norm *lowband_out,
1203       celt_norm *lowband_scratch, int fill)
1204 {
1205    int imid=0, iside=0;
1206    int inv = 0;
1207    opus_val16 mid=0, side=0;
1208    unsigned cm=0;
1209 #ifdef RESYNTH
1210    int resynth = 1;
1211 #else
1212    int resynth = !ctx->encode;
1213 #endif
1214    int mbits, sbits, delta;
1215    int itheta;
1216    int qalloc;
1217    struct split_ctx sctx;
1218    int orig_fill;
1219    int encode;
1220    ec_ctx *ec;
1221
1222    encode = ctx->encode;
1223    ec = ctx->ec;
1224
1225    /* Special case for one sample */
1226    if (N==1)
1227    {
1228       return quant_band_n1(ctx, X, Y, b, lowband_out);
1229    }
1230
1231    orig_fill = fill;
1232
1233    compute_theta(ctx, &sctx, X, Y, N, &b, B, B,
1234          LM, 1, &fill);
1235    inv = sctx.inv;
1236    imid = sctx.imid;
1237    iside = sctx.iside;
1238    delta = sctx.delta;
1239    itheta = sctx.itheta;
1240    qalloc = sctx.qalloc;
1241 #ifdef FIXED_POINT
1242    mid = imid;
1243    side = iside;
1244 #else
1245    mid = (1.f/32768)*imid;
1246    side = (1.f/32768)*iside;
1247 #endif
1248
1249    /* This is a special case for N=2 that only works for stereo and takes
1250       advantage of the fact that mid and side are orthogonal to encode
1251       the side with just one bit. */
1252    if (N==2)
1253    {
1254       int c;
1255       int sign=0;
1256       celt_norm *x2, *y2;
1257       mbits = b;
1258       sbits = 0;
1259       /* Only need one bit for the side. */
1260       if (itheta != 0 && itheta != 16384)
1261          sbits = 1<<BITRES;
1262       mbits -= sbits;
1263       c = itheta > 8192;
1264       ctx->remaining_bits -= qalloc+sbits;
1265
1266       x2 = c ? Y : X;
1267       y2 = c ? X : Y;
1268       if (sbits)
1269       {
1270          if (encode)
1271          {
1272             /* Here we only need to encode a sign for the side. */
1273             sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
1274             ec_enc_bits(ec, sign, 1);
1275          } else {
1276             sign = ec_dec_bits(ec, 1);
1277          }
1278       }
1279       sign = 1-2*sign;
1280       /* We use orig_fill here because we want to fold the side, but if
1281          itheta==16384, we'll have cleared the low bits of fill. */
1282       cm = quant_band(ctx, x2, N, mbits, B, lowband,
1283             LM, lowband_out, Q15ONE, lowband_scratch, orig_fill);
1284       /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
1285          and there's no need to worry about mixing with the other channel. */
1286       y2[0] = -sign*x2[1];
1287       y2[1] = sign*x2[0];
1288       if (resynth)
1289       {
1290          celt_norm tmp;
1291          X[0] = MULT16_16_Q15(mid, X[0]);
1292          X[1] = MULT16_16_Q15(mid, X[1]);
1293          Y[0] = MULT16_16_Q15(side, Y[0]);
1294          Y[1] = MULT16_16_Q15(side, Y[1]);
1295          tmp = X[0];
1296          X[0] = SUB16(tmp,Y[0]);
1297          Y[0] = ADD16(tmp,Y[0]);
1298          tmp = X[1];
1299          X[1] = SUB16(tmp,Y[1]);
1300          Y[1] = ADD16(tmp,Y[1]);
1301       }
1302    } else {
1303       /* "Normal" split code */
1304       opus_int32 rebalance;
1305
1306       mbits = IMAX(0, IMIN(b, (b-delta)/2));
1307       sbits = b-mbits;
1308       ctx->remaining_bits -= qalloc;
1309
1310       rebalance = ctx->remaining_bits;
1311       if (mbits >= sbits)
1312       {
1313          /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1314             mid for folding later. */
1315          cm = quant_band(ctx, X, N, mbits, B,
1316                lowband, LM, lowband_out,
1317                Q15ONE, lowband_scratch, fill);
1318          rebalance = mbits - (rebalance-ctx->remaining_bits);
1319          if (rebalance > 3<<BITRES && itheta!=0)
1320             sbits += rebalance - (3<<BITRES);
1321
1322          /* For a stereo split, the high bits of fill are always zero, so no
1323             folding will be done to the side. */
1324          cm |= quant_band(ctx, Y, N, sbits, B,
1325                NULL, LM, NULL,
1326                side, NULL, fill>>B);
1327       } else {
1328          /* For a stereo split, the high bits of fill are always zero, so no
1329             folding will be done to the side. */
1330          cm = quant_band(ctx, Y, N, sbits, B,
1331                NULL, LM, NULL,
1332                side, NULL, fill>>B);
1333          rebalance = sbits - (rebalance-ctx->remaining_bits);
1334          if (rebalance > 3<<BITRES && itheta!=16384)
1335             mbits += rebalance - (3<<BITRES);
1336          /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1337             mid for folding later. */
1338          cm |= quant_band(ctx, X, N, mbits, B,
1339                lowband, LM, lowband_out,
1340                Q15ONE, lowband_scratch, fill);
1341       }
1342    }
1343
1344
1345    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1346    if (resynth)
1347    {
1348       if (N!=2)
1349          stereo_merge(X, Y, mid, N);
1350       if (inv)
1351       {
1352          int j;
1353          for (j=0;j<N;j++)
1354             Y[j] = -Y[j];
1355       }
1356    }
1357    return cm;
1358 }
1359
1360
1361 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1362       celt_norm *X_, celt_norm *Y_, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1363       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res,
1364       opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed)
1365 {
1366    int i;
1367    opus_int32 remaining_bits;
1368    const opus_int16 * OPUS_RESTRICT eBands = m->eBands;
1369    celt_norm * OPUS_RESTRICT norm, * OPUS_RESTRICT norm2;
1370    VARDECL(celt_norm, _norm);
1371    celt_norm *lowband_scratch;
1372    int B;
1373    int M;
1374    int lowband_offset;
1375    int update_lowband = 1;
1376    int C = Y_ != NULL ? 2 : 1;
1377    int norm_offset;
1378 #ifdef RESYNTH
1379    int resynth = 1;
1380 #else
1381    int resynth = !encode;
1382 #endif
1383    struct band_ctx ctx;
1384    SAVE_STACK;
1385
1386    M = 1<<LM;
1387    B = shortBlocks ? M : 1;
1388    norm_offset = M*eBands[start];
1389    /* No need to allocate norm for the last band because we don't need an
1390       output in that band. */
1391    ALLOC(_norm, C*(M*eBands[m->nbEBands-1]-norm_offset), celt_norm);
1392    norm = _norm;
1393    norm2 = norm + M*eBands[m->nbEBands-1]-norm_offset;
1394    /* We can use the last band as scratch space because we don't need that
1395       scratch space for the last band. */
1396    lowband_scratch = X_+M*eBands[m->nbEBands-1];
1397
1398    lowband_offset = 0;
1399    ctx.bandE = bandE;
1400    ctx.ec = ec;
1401    ctx.encode = encode;
1402    ctx.intensity = intensity;
1403    ctx.m = m;
1404    ctx.seed = *seed;
1405    ctx.spread = spread;
1406    for (i=start;i<end;i++)
1407    {
1408       opus_int32 tell;
1409       int b;
1410       int N;
1411       opus_int32 curr_balance;
1412       int effective_lowband=-1;
1413       celt_norm * OPUS_RESTRICT X, * OPUS_RESTRICT Y;
1414       int tf_change=0;
1415       unsigned x_cm;
1416       unsigned y_cm;
1417       int last;
1418
1419       ctx.i = i;
1420       last = (i==end-1);
1421
1422       X = X_+M*eBands[i];
1423       if (Y_!=NULL)
1424          Y = Y_+M*eBands[i];
1425       else
1426          Y = NULL;
1427       N = M*eBands[i+1]-M*eBands[i];
1428       tell = ec_tell_frac(ec);
1429
1430       /* Compute how many bits we want to allocate to this band */
1431       if (i != start)
1432          balance -= tell;
1433       remaining_bits = total_bits-tell-1;
1434       ctx.remaining_bits = remaining_bits;
1435       if (i <= codedBands-1)
1436       {
1437          curr_balance = balance / IMIN(3, codedBands-i);
1438          b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1439       } else {
1440          b = 0;
1441       }
1442
1443       if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1444             lowband_offset = i;
1445
1446       tf_change = tf_res[i];
1447       ctx.tf_change = tf_change;
1448       if (i>=m->effEBands)
1449       {
1450          X=norm;
1451          if (Y_!=NULL)
1452             Y = norm;
1453          lowband_scratch = NULL;
1454       }
1455       if (i==end-1)
1456          lowband_scratch = NULL;
1457
1458       /* Get a conservative estimate of the collapse_mask's for the bands we're
1459          going to be folding from. */
1460       if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1461       {
1462          int fold_start;
1463          int fold_end;
1464          int fold_i;
1465          /* This ensures we never repeat spectral content within one band */
1466          effective_lowband = IMAX(0, M*eBands[lowband_offset]-norm_offset-N);
1467          fold_start = lowband_offset;
1468          while(M*eBands[--fold_start] > effective_lowband+norm_offset);
1469          fold_end = lowband_offset-1;
1470          while(M*eBands[++fold_end] < effective_lowband+norm_offset+N);
1471          x_cm = y_cm = 0;
1472          fold_i = fold_start; do {
1473            x_cm |= collapse_masks[fold_i*C+0];
1474            y_cm |= collapse_masks[fold_i*C+C-1];
1475          } while (++fold_i<fold_end);
1476       }
1477       /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1478          always) be non-zero. */
1479       else
1480          x_cm = y_cm = (1<<B)-1;
1481
1482       if (dual_stereo && i==intensity)
1483       {
1484          int j;
1485
1486          /* Switch off dual stereo to do intensity. */
1487          dual_stereo = 0;
1488          if (resynth)
1489             for (j=0;j<M*eBands[i]-norm_offset;j++)
1490                norm[j] = HALF32(norm[j]+norm2[j]);
1491       }
1492       if (dual_stereo)
1493       {
1494          x_cm = quant_band(&ctx, X, N, b/2, B,
1495                effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1496                last?NULL:norm+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, x_cm);
1497          y_cm = quant_band(&ctx, Y, N, b/2, B,
1498                effective_lowband != -1 ? norm2+effective_lowband : NULL, LM,
1499                last?NULL:norm2+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, y_cm);
1500       } else {
1501          if (Y!=NULL)
1502          {
1503             x_cm = quant_band_stereo(&ctx, X, Y, N, b, B,
1504                   effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1505                         last?NULL:norm+M*eBands[i]-norm_offset, lowband_scratch, x_cm|y_cm);
1506          } else {
1507             x_cm = quant_band(&ctx, X, N, b, B,
1508                   effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1509                         last?NULL:norm+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, x_cm|y_cm);
1510          }
1511          y_cm = x_cm;
1512       }
1513       collapse_masks[i*C+0] = (unsigned char)x_cm;
1514       collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1515       balance += pulses[i] + tell;
1516
1517       /* Update the folding position only as long as we have 1 bit/sample depth. */
1518       update_lowband = b>(N<<BITRES);
1519    }
1520    *seed = ctx.seed;
1521
1522    RESTORE_STACK;
1523 }
1524