Optimizing divisions with a signed numerator
[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 = celt_sudiv(b+N2*offset, N2);
634    qb = IMIN(b-pulse_cap-(4<<BITRES), qb);
635
636    qb = IMIN(8<<BITRES, qb);
637
638    if (qb<(1<<BITRES>>1)) {
639       qn = 1;
640    } else {
641       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
642       qn = (qn+1)>>1<<1;
643    }
644    celt_assert(qn <= 256);
645    return qn;
646 }
647
648 struct band_ctx {
649    int encode;
650    const CELTMode *m;
651    int i;
652    int intensity;
653    int spread;
654    int tf_change;
655    ec_ctx *ec;
656    opus_int32 remaining_bits;
657    const celt_ener *bandE;
658    opus_uint32 seed;
659 };
660
661 struct split_ctx {
662    int inv;
663    int imid;
664    int iside;
665    int delta;
666    int itheta;
667    int qalloc;
668 };
669
670 static void compute_theta(struct band_ctx *ctx, struct split_ctx *sctx,
671       celt_norm *X, celt_norm *Y, int N, int *b, int B, int B0,
672       int LM,
673       int stereo, int *fill)
674 {
675    int qn;
676    int itheta=0;
677    int delta;
678    int imid, iside;
679    int qalloc;
680    int pulse_cap;
681    int offset;
682    opus_int32 tell;
683    int inv=0;
684    int encode;
685    const CELTMode *m;
686    int i;
687    int intensity;
688    ec_ctx *ec;
689    const celt_ener *bandE;
690
691    encode = ctx->encode;
692    m = ctx->m;
693    i = ctx->i;
694    intensity = ctx->intensity;
695    ec = ctx->ec;
696    bandE = ctx->bandE;
697
698    /* Decide on the resolution to give to the split parameter theta */
699    pulse_cap = m->logN[i]+LM*(1<<BITRES);
700    offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
701    qn = compute_qn(N, *b, offset, pulse_cap, stereo);
702    if (stereo && i>=intensity)
703       qn = 1;
704    if (encode)
705    {
706       /* theta is the atan() of the ratio between the (normalized)
707          side and mid. With just that parameter, we can re-scale both
708          mid and side because we know that 1) they have unit norm and
709          2) they are orthogonal. */
710       itheta = stereo_itheta(X, Y, stereo, N);
711    }
712    tell = ec_tell_frac(ec);
713    if (qn!=1)
714    {
715       if (encode)
716          itheta = (itheta*qn+8192)>>14;
717
718       /* Entropy coding of the angle. We use a uniform pdf for the
719          time split, a step for stereo, and a triangular one for the rest. */
720       if (stereo && N>2)
721       {
722          int p0 = 3;
723          int x = itheta;
724          int x0 = qn/2;
725          int ft = p0*(x0+1) + x0;
726          /* Use a probability of p0 up to itheta=8192 and then use 1 after */
727          if (encode)
728          {
729             ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
730          } else {
731             int fs;
732             fs=ec_decode(ec,ft);
733             if (fs<(x0+1)*p0)
734                x=fs/p0;
735             else
736                x=x0+1+(fs-(x0+1)*p0);
737             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);
738             itheta = x;
739          }
740       } else if (B0>1 || stereo) {
741          /* Uniform pdf */
742          if (encode)
743             ec_enc_uint(ec, itheta, qn+1);
744          else
745             itheta = ec_dec_uint(ec, qn+1);
746       } else {
747          int fs=1, ft;
748          ft = ((qn>>1)+1)*((qn>>1)+1);
749          if (encode)
750          {
751             int fl;
752
753             fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
754             fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
755              ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
756
757             ec_encode(ec, fl, fl+fs, ft);
758          } else {
759             /* Triangular pdf */
760             int fl=0;
761             int fm;
762             fm = ec_decode(ec, ft);
763
764             if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
765             {
766                itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;
767                fs = itheta + 1;
768                fl = itheta*(itheta + 1)>>1;
769             }
770             else
771             {
772                itheta = (2*(qn + 1)
773                 - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;
774                fs = qn + 1 - itheta;
775                fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
776             }
777
778             ec_dec_update(ec, fl, fl+fs, ft);
779          }
780       }
781       celt_assert(itheta>=0);
782       itheta = celt_udiv((opus_int32)itheta*16384, qn);
783       if (encode && stereo)
784       {
785          if (itheta==0)
786             intensity_stereo(m, X, Y, bandE, i, N);
787          else
788             stereo_split(X, Y, N);
789       }
790       /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.
791                Let's do that at higher complexity */
792    } else if (stereo) {
793       if (encode)
794       {
795          inv = itheta > 8192;
796          if (inv)
797          {
798             int j;
799             for (j=0;j<N;j++)
800                Y[j] = -Y[j];
801          }
802          intensity_stereo(m, X, Y, bandE, i, N);
803       }
804       if (*b>2<<BITRES && ctx->remaining_bits > 2<<BITRES)
805       {
806          if (encode)
807             ec_enc_bit_logp(ec, inv, 2);
808          else
809             inv = ec_dec_bit_logp(ec, 2);
810       } else
811          inv = 0;
812       itheta = 0;
813    }
814    qalloc = ec_tell_frac(ec) - tell;
815    *b -= qalloc;
816
817    if (itheta == 0)
818    {
819       imid = 32767;
820       iside = 0;
821       *fill &= (1<<B)-1;
822       delta = -16384;
823    } else if (itheta == 16384)
824    {
825       imid = 0;
826       iside = 32767;
827       *fill &= ((1<<B)-1)<<B;
828       delta = 16384;
829    } else {
830       imid = bitexact_cos((opus_int16)itheta);
831       iside = bitexact_cos((opus_int16)(16384-itheta));
832       /* This is the mid vs side allocation that minimizes squared error
833          in that band. */
834       delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid));
835    }
836
837    sctx->inv = inv;
838    sctx->imid = imid;
839    sctx->iside = iside;
840    sctx->delta = delta;
841    sctx->itheta = itheta;
842    sctx->qalloc = qalloc;
843 }
844 static unsigned quant_band_n1(struct band_ctx *ctx, celt_norm *X, celt_norm *Y, int b,
845       celt_norm *lowband_out)
846 {
847 #ifdef RESYNTH
848    int resynth = 1;
849 #else
850    int resynth = !ctx->encode;
851 #endif
852    int c;
853    int stereo;
854    celt_norm *x = X;
855    int encode;
856    ec_ctx *ec;
857
858    encode = ctx->encode;
859    ec = ctx->ec;
860
861    stereo = Y != NULL;
862    c=0; do {
863       int sign=0;
864       if (ctx->remaining_bits>=1<<BITRES)
865       {
866          if (encode)
867          {
868             sign = x[0]<0;
869             ec_enc_bits(ec, sign, 1);
870          } else {
871             sign = ec_dec_bits(ec, 1);
872          }
873          ctx->remaining_bits -= 1<<BITRES;
874          b-=1<<BITRES;
875       }
876       if (resynth)
877          x[0] = sign ? -NORM_SCALING : NORM_SCALING;
878       x = Y;
879    } while (++c<1+stereo);
880    if (lowband_out)
881       lowband_out[0] = SHR16(X[0],4);
882    return 1;
883 }
884
885 /* This function is responsible for encoding and decoding a mono partition.
886    It can split the band in two and transmit the energy difference with
887    the two half-bands. It can be called recursively so bands can end up being
888    split in 8 parts. */
889 static unsigned quant_partition(struct band_ctx *ctx, celt_norm *X,
890       int N, int b, int B, celt_norm *lowband,
891       int LM,
892       opus_val16 gain, int fill)
893 {
894    const unsigned char *cache;
895    int q;
896    int curr_bits;
897    int imid=0, iside=0;
898    int B0=B;
899    opus_val16 mid=0, side=0;
900    unsigned cm=0;
901 #ifdef RESYNTH
902    int resynth = 1;
903 #else
904    int resynth = !ctx->encode;
905 #endif
906    celt_norm *Y=NULL;
907    int encode;
908    const CELTMode *m;
909    int i;
910    int spread;
911    ec_ctx *ec;
912
913    encode = ctx->encode;
914    m = ctx->m;
915    i = ctx->i;
916    spread = ctx->spread;
917    ec = ctx->ec;
918
919    /* If we need 1.5 more bit than we can produce, split the band in two. */
920    cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
921    if (LM != -1 && b > cache[cache[0]]+12 && N>2)
922    {
923       int mbits, sbits, delta;
924       int itheta;
925       int qalloc;
926       struct split_ctx sctx;
927       celt_norm *next_lowband2=NULL;
928       opus_int32 rebalance;
929
930       N >>= 1;
931       Y = X+N;
932       LM -= 1;
933       if (B==1)
934          fill = (fill&1)|(fill<<1);
935       B = (B+1)>>1;
936
937       compute_theta(ctx, &sctx, X, Y, N, &b, B, B0,
938             LM, 0, &fill);
939       imid = sctx.imid;
940       iside = sctx.iside;
941       delta = sctx.delta;
942       itheta = sctx.itheta;
943       qalloc = sctx.qalloc;
944 #ifdef FIXED_POINT
945       mid = imid;
946       side = iside;
947 #else
948       mid = (1.f/32768)*imid;
949       side = (1.f/32768)*iside;
950 #endif
951
952       /* Give more bits to low-energy MDCTs than they would otherwise deserve */
953       if (B0>1 && (itheta&0x3fff))
954       {
955          if (itheta > 8192)
956             /* Rough approximation for pre-echo masking */
957             delta -= delta>>(4-LM);
958          else
959             /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
960             delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
961       }
962       mbits = IMAX(0, IMIN(b, (b-delta)/2));
963       sbits = b-mbits;
964       ctx->remaining_bits -= qalloc;
965
966       if (lowband)
967          next_lowband2 = lowband+N; /* >32-bit split case */
968
969       rebalance = ctx->remaining_bits;
970       if (mbits >= sbits)
971       {
972          cm = quant_partition(ctx, X, N, mbits, B,
973                lowband, LM,
974                MULT16_16_P15(gain,mid), fill);
975          rebalance = mbits - (rebalance-ctx->remaining_bits);
976          if (rebalance > 3<<BITRES && itheta!=0)
977             sbits += rebalance - (3<<BITRES);
978          cm |= quant_partition(ctx, Y, N, sbits, B,
979                next_lowband2, LM,
980                MULT16_16_P15(gain,side), fill>>B)<<(B0>>1);
981       } else {
982          cm = quant_partition(ctx, Y, N, sbits, B,
983                next_lowband2, LM,
984                MULT16_16_P15(gain,side), fill>>B)<<(B0>>1);
985          rebalance = sbits - (rebalance-ctx->remaining_bits);
986          if (rebalance > 3<<BITRES && itheta!=16384)
987             mbits += rebalance - (3<<BITRES);
988          cm |= quant_partition(ctx, X, N, mbits, B,
989                lowband, LM,
990                MULT16_16_P15(gain,mid), fill);
991       }
992    } else {
993       /* This is the basic no-split case */
994       q = bits2pulses(m, i, LM, b);
995       curr_bits = pulses2bits(m, i, LM, q);
996       ctx->remaining_bits -= curr_bits;
997
998       /* Ensures we can never bust the budget */
999       while (ctx->remaining_bits < 0 && q > 0)
1000       {
1001          ctx->remaining_bits += curr_bits;
1002          q--;
1003          curr_bits = pulses2bits(m, i, LM, q);
1004          ctx->remaining_bits -= curr_bits;
1005       }
1006
1007       if (q!=0)
1008       {
1009          int K = get_pulses(q);
1010
1011          /* Finally do the actual quantization */
1012          if (encode)
1013          {
1014             cm = alg_quant(X, N, K, spread, B, ec
1015 #ifdef RESYNTH
1016                  , gain
1017 #endif
1018                  );
1019          } else {
1020             cm = alg_unquant(X, N, K, spread, B, ec, gain);
1021          }
1022       } else {
1023          /* If there's no pulse, fill the band anyway */
1024          int j;
1025          if (resynth)
1026          {
1027             unsigned cm_mask;
1028             /* B can be as large as 16, so this shift might overflow an int on a
1029                16-bit platform; use a long to get defined behavior.*/
1030             cm_mask = (unsigned)(1UL<<B)-1;
1031             fill &= cm_mask;
1032             if (!fill)
1033             {
1034                OPUS_CLEAR(X, N);
1035             } else {
1036                if (lowband == NULL)
1037                {
1038                   /* Noise */
1039                   for (j=0;j<N;j++)
1040                   {
1041                      ctx->seed = celt_lcg_rand(ctx->seed);
1042                      X[j] = (celt_norm)((opus_int32)ctx->seed>>20);
1043                   }
1044                   cm = cm_mask;
1045                } else {
1046                   /* Folded spectrum */
1047                   for (j=0;j<N;j++)
1048                   {
1049                      opus_val16 tmp;
1050                      ctx->seed = celt_lcg_rand(ctx->seed);
1051                      /* About 48 dB below the "normal" folding level */
1052                      tmp = QCONST16(1.0f/256, 10);
1053                      tmp = (ctx->seed)&0x8000 ? tmp : -tmp;
1054                      X[j] = lowband[j]+tmp;
1055                   }
1056                   cm = fill;
1057                }
1058                renormalise_vector(X, N, gain);
1059             }
1060          }
1061       }
1062    }
1063
1064    return cm;
1065 }
1066
1067
1068 /* This function is responsible for encoding and decoding a band for the mono case. */
1069 static unsigned quant_band(struct band_ctx *ctx, celt_norm *X,
1070       int N, int b, int B, celt_norm *lowband,
1071       int LM, celt_norm *lowband_out,
1072       opus_val16 gain, celt_norm *lowband_scratch, int fill)
1073 {
1074    int N0=N;
1075    int N_B=N;
1076    int N_B0;
1077    int B0=B;
1078    int time_divide=0;
1079    int recombine=0;
1080    int longBlocks;
1081    unsigned cm=0;
1082 #ifdef RESYNTH
1083    int resynth = 1;
1084 #else
1085    int resynth = !ctx->encode;
1086 #endif
1087    int k;
1088    int encode;
1089    int tf_change;
1090
1091    encode = ctx->encode;
1092    tf_change = ctx->tf_change;
1093
1094    longBlocks = B0==1;
1095
1096    N_B = celt_udiv(N_B, B);
1097
1098    /* Special case for one sample */
1099    if (N==1)
1100    {
1101       return quant_band_n1(ctx, X, NULL, b, lowband_out);
1102    }
1103
1104    if (tf_change>0)
1105       recombine = tf_change;
1106    /* Band recombining to increase frequency resolution */
1107
1108    if (lowband_scratch && lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
1109    {
1110       OPUS_COPY(lowband_scratch, lowband, N);
1111       lowband = lowband_scratch;
1112    }
1113
1114    for (k=0;k<recombine;k++)
1115    {
1116       static const unsigned char bit_interleave_table[16]={
1117             0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
1118       };
1119       if (encode)
1120          haar1(X, N>>k, 1<<k);
1121       if (lowband)
1122          haar1(lowband, N>>k, 1<<k);
1123       fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
1124    }
1125    B>>=recombine;
1126    N_B<<=recombine;
1127
1128    /* Increasing the time resolution */
1129    while ((N_B&1) == 0 && tf_change<0)
1130    {
1131       if (encode)
1132          haar1(X, N_B, B);
1133       if (lowband)
1134          haar1(lowband, N_B, B);
1135       fill |= fill<<B;
1136       B <<= 1;
1137       N_B >>= 1;
1138       time_divide++;
1139       tf_change++;
1140    }
1141    B0=B;
1142    N_B0 = N_B;
1143
1144    /* Reorganize the samples in time order instead of frequency order */
1145    if (B0>1)
1146    {
1147       if (encode)
1148          deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1149       if (lowband)
1150          deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
1151    }
1152
1153    cm = quant_partition(ctx, X, N, b, B, lowband,
1154          LM, gain, fill);
1155
1156    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1157    if (resynth)
1158    {
1159       /* Undo the sample reorganization going from time order to frequency order */
1160       if (B0>1)
1161          interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1162
1163       /* Undo time-freq changes that we did earlier */
1164       N_B = N_B0;
1165       B = B0;
1166       for (k=0;k<time_divide;k++)
1167       {
1168          B >>= 1;
1169          N_B <<= 1;
1170          cm |= cm>>B;
1171          haar1(X, N_B, B);
1172       }
1173
1174       for (k=0;k<recombine;k++)
1175       {
1176          static const unsigned char bit_deinterleave_table[16]={
1177                0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1178                0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1179          };
1180          cm = bit_deinterleave_table[cm];
1181          haar1(X, N0>>k, 1<<k);
1182       }
1183       B<<=recombine;
1184
1185       /* Scale output for later folding */
1186       if (lowband_out)
1187       {
1188          int j;
1189          opus_val16 n;
1190          n = celt_sqrt(SHL32(EXTEND32(N0),22));
1191          for (j=0;j<N0;j++)
1192             lowband_out[j] = MULT16_16_Q15(n,X[j]);
1193       }
1194       cm &= (1<<B)-1;
1195    }
1196    return cm;
1197 }
1198
1199
1200 /* This function is responsible for encoding and decoding a band for the stereo case. */
1201 static unsigned quant_band_stereo(struct band_ctx *ctx, celt_norm *X, celt_norm *Y,
1202       int N, int b, int B, celt_norm *lowband,
1203       int LM, celt_norm *lowband_out,
1204       celt_norm *lowband_scratch, int fill)
1205 {
1206    int imid=0, iside=0;
1207    int inv = 0;
1208    opus_val16 mid=0, side=0;
1209    unsigned cm=0;
1210 #ifdef RESYNTH
1211    int resynth = 1;
1212 #else
1213    int resynth = !ctx->encode;
1214 #endif
1215    int mbits, sbits, delta;
1216    int itheta;
1217    int qalloc;
1218    struct split_ctx sctx;
1219    int orig_fill;
1220    int encode;
1221    ec_ctx *ec;
1222
1223    encode = ctx->encode;
1224    ec = ctx->ec;
1225
1226    /* Special case for one sample */
1227    if (N==1)
1228    {
1229       return quant_band_n1(ctx, X, Y, b, lowband_out);
1230    }
1231
1232    orig_fill = fill;
1233
1234    compute_theta(ctx, &sctx, X, Y, N, &b, B, B,
1235          LM, 1, &fill);
1236    inv = sctx.inv;
1237    imid = sctx.imid;
1238    iside = sctx.iside;
1239    delta = sctx.delta;
1240    itheta = sctx.itheta;
1241    qalloc = sctx.qalloc;
1242 #ifdef FIXED_POINT
1243    mid = imid;
1244    side = iside;
1245 #else
1246    mid = (1.f/32768)*imid;
1247    side = (1.f/32768)*iside;
1248 #endif
1249
1250    /* This is a special case for N=2 that only works for stereo and takes
1251       advantage of the fact that mid and side are orthogonal to encode
1252       the side with just one bit. */
1253    if (N==2)
1254    {
1255       int c;
1256       int sign=0;
1257       celt_norm *x2, *y2;
1258       mbits = b;
1259       sbits = 0;
1260       /* Only need one bit for the side. */
1261       if (itheta != 0 && itheta != 16384)
1262          sbits = 1<<BITRES;
1263       mbits -= sbits;
1264       c = itheta > 8192;
1265       ctx->remaining_bits -= qalloc+sbits;
1266
1267       x2 = c ? Y : X;
1268       y2 = c ? X : Y;
1269       if (sbits)
1270       {
1271          if (encode)
1272          {
1273             /* Here we only need to encode a sign for the side. */
1274             sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
1275             ec_enc_bits(ec, sign, 1);
1276          } else {
1277             sign = ec_dec_bits(ec, 1);
1278          }
1279       }
1280       sign = 1-2*sign;
1281       /* We use orig_fill here because we want to fold the side, but if
1282          itheta==16384, we'll have cleared the low bits of fill. */
1283       cm = quant_band(ctx, x2, N, mbits, B, lowband,
1284             LM, lowband_out, Q15ONE, lowband_scratch, orig_fill);
1285       /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
1286          and there's no need to worry about mixing with the other channel. */
1287       y2[0] = -sign*x2[1];
1288       y2[1] = sign*x2[0];
1289       if (resynth)
1290       {
1291          celt_norm tmp;
1292          X[0] = MULT16_16_Q15(mid, X[0]);
1293          X[1] = MULT16_16_Q15(mid, X[1]);
1294          Y[0] = MULT16_16_Q15(side, Y[0]);
1295          Y[1] = MULT16_16_Q15(side, Y[1]);
1296          tmp = X[0];
1297          X[0] = SUB16(tmp,Y[0]);
1298          Y[0] = ADD16(tmp,Y[0]);
1299          tmp = X[1];
1300          X[1] = SUB16(tmp,Y[1]);
1301          Y[1] = ADD16(tmp,Y[1]);
1302       }
1303    } else {
1304       /* "Normal" split code */
1305       opus_int32 rebalance;
1306
1307       mbits = IMAX(0, IMIN(b, (b-delta)/2));
1308       sbits = b-mbits;
1309       ctx->remaining_bits -= qalloc;
1310
1311       rebalance = ctx->remaining_bits;
1312       if (mbits >= sbits)
1313       {
1314          /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1315             mid for folding later. */
1316          cm = quant_band(ctx, X, N, mbits, B,
1317                lowband, LM, lowband_out,
1318                Q15ONE, lowband_scratch, fill);
1319          rebalance = mbits - (rebalance-ctx->remaining_bits);
1320          if (rebalance > 3<<BITRES && itheta!=0)
1321             sbits += rebalance - (3<<BITRES);
1322
1323          /* For a stereo split, the high bits of fill are always zero, so no
1324             folding will be done to the side. */
1325          cm |= quant_band(ctx, Y, N, sbits, B,
1326                NULL, LM, NULL,
1327                side, NULL, fill>>B);
1328       } else {
1329          /* For a stereo split, the high bits of fill are always zero, so no
1330             folding will be done to the side. */
1331          cm = quant_band(ctx, Y, N, sbits, B,
1332                NULL, LM, NULL,
1333                side, NULL, fill>>B);
1334          rebalance = sbits - (rebalance-ctx->remaining_bits);
1335          if (rebalance > 3<<BITRES && itheta!=16384)
1336             mbits += rebalance - (3<<BITRES);
1337          /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1338             mid for folding later. */
1339          cm |= quant_band(ctx, X, N, mbits, B,
1340                lowband, LM, lowband_out,
1341                Q15ONE, lowband_scratch, fill);
1342       }
1343    }
1344
1345
1346    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1347    if (resynth)
1348    {
1349       if (N!=2)
1350          stereo_merge(X, Y, mid, N);
1351       if (inv)
1352       {
1353          int j;
1354          for (j=0;j<N;j++)
1355             Y[j] = -Y[j];
1356       }
1357    }
1358    return cm;
1359 }
1360
1361
1362 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1363       celt_norm *X_, celt_norm *Y_, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1364       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res,
1365       opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed)
1366 {
1367    int i;
1368    opus_int32 remaining_bits;
1369    const opus_int16 * OPUS_RESTRICT eBands = m->eBands;
1370    celt_norm * OPUS_RESTRICT norm, * OPUS_RESTRICT norm2;
1371    VARDECL(celt_norm, _norm);
1372    celt_norm *lowband_scratch;
1373    int B;
1374    int M;
1375    int lowband_offset;
1376    int update_lowband = 1;
1377    int C = Y_ != NULL ? 2 : 1;
1378    int norm_offset;
1379 #ifdef RESYNTH
1380    int resynth = 1;
1381 #else
1382    int resynth = !encode;
1383 #endif
1384    struct band_ctx ctx;
1385    SAVE_STACK;
1386
1387    M = 1<<LM;
1388    B = shortBlocks ? M : 1;
1389    norm_offset = M*eBands[start];
1390    /* No need to allocate norm for the last band because we don't need an
1391       output in that band. */
1392    ALLOC(_norm, C*(M*eBands[m->nbEBands-1]-norm_offset), celt_norm);
1393    norm = _norm;
1394    norm2 = norm + M*eBands[m->nbEBands-1]-norm_offset;
1395    /* We can use the last band as scratch space because we don't need that
1396       scratch space for the last band. */
1397    lowband_scratch = X_+M*eBands[m->nbEBands-1];
1398
1399    lowband_offset = 0;
1400    ctx.bandE = bandE;
1401    ctx.ec = ec;
1402    ctx.encode = encode;
1403    ctx.intensity = intensity;
1404    ctx.m = m;
1405    ctx.seed = *seed;
1406    ctx.spread = spread;
1407    for (i=start;i<end;i++)
1408    {
1409       opus_int32 tell;
1410       int b;
1411       int N;
1412       opus_int32 curr_balance;
1413       int effective_lowband=-1;
1414       celt_norm * OPUS_RESTRICT X, * OPUS_RESTRICT Y;
1415       int tf_change=0;
1416       unsigned x_cm;
1417       unsigned y_cm;
1418       int last;
1419
1420       ctx.i = i;
1421       last = (i==end-1);
1422
1423       X = X_+M*eBands[i];
1424       if (Y_!=NULL)
1425          Y = Y_+M*eBands[i];
1426       else
1427          Y = NULL;
1428       N = M*eBands[i+1]-M*eBands[i];
1429       tell = ec_tell_frac(ec);
1430
1431       /* Compute how many bits we want to allocate to this band */
1432       if (i != start)
1433          balance -= tell;
1434       remaining_bits = total_bits-tell-1;
1435       ctx.remaining_bits = remaining_bits;
1436       if (i <= codedBands-1)
1437       {
1438          curr_balance = celt_sudiv(balance, IMIN(3, codedBands-i));
1439          b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1440       } else {
1441          b = 0;
1442       }
1443
1444       if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1445             lowband_offset = i;
1446
1447       tf_change = tf_res[i];
1448       ctx.tf_change = tf_change;
1449       if (i>=m->effEBands)
1450       {
1451          X=norm;
1452          if (Y_!=NULL)
1453             Y = norm;
1454          lowband_scratch = NULL;
1455       }
1456       if (i==end-1)
1457          lowband_scratch = NULL;
1458
1459       /* Get a conservative estimate of the collapse_mask's for the bands we're
1460          going to be folding from. */
1461       if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1462       {
1463          int fold_start;
1464          int fold_end;
1465          int fold_i;
1466          /* This ensures we never repeat spectral content within one band */
1467          effective_lowband = IMAX(0, M*eBands[lowband_offset]-norm_offset-N);
1468          fold_start = lowband_offset;
1469          while(M*eBands[--fold_start] > effective_lowband+norm_offset);
1470          fold_end = lowband_offset-1;
1471          while(M*eBands[++fold_end] < effective_lowband+norm_offset+N);
1472          x_cm = y_cm = 0;
1473          fold_i = fold_start; do {
1474            x_cm |= collapse_masks[fold_i*C+0];
1475            y_cm |= collapse_masks[fold_i*C+C-1];
1476          } while (++fold_i<fold_end);
1477       }
1478       /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1479          always) be non-zero. */
1480       else
1481          x_cm = y_cm = (1<<B)-1;
1482
1483       if (dual_stereo && i==intensity)
1484       {
1485          int j;
1486
1487          /* Switch off dual stereo to do intensity. */
1488          dual_stereo = 0;
1489          if (resynth)
1490             for (j=0;j<M*eBands[i]-norm_offset;j++)
1491                norm[j] = HALF32(norm[j]+norm2[j]);
1492       }
1493       if (dual_stereo)
1494       {
1495          x_cm = quant_band(&ctx, X, N, b/2, B,
1496                effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1497                last?NULL:norm+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, x_cm);
1498          y_cm = quant_band(&ctx, Y, N, b/2, B,
1499                effective_lowband != -1 ? norm2+effective_lowband : NULL, LM,
1500                last?NULL:norm2+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, y_cm);
1501       } else {
1502          if (Y!=NULL)
1503          {
1504             x_cm = quant_band_stereo(&ctx, X, Y, N, b, B,
1505                   effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1506                         last?NULL:norm+M*eBands[i]-norm_offset, lowband_scratch, x_cm|y_cm);
1507          } else {
1508             x_cm = quant_band(&ctx, X, N, b, B,
1509                   effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1510                         last?NULL:norm+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, x_cm|y_cm);
1511          }
1512          y_cm = x_cm;
1513       }
1514       collapse_masks[i*C+0] = (unsigned char)x_cm;
1515       collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1516       balance += pulses[i] + tell;
1517
1518       /* Update the folding position only as long as we have 1 bit/sample depth. */
1519       update_lowband = b>(N<<BITRES);
1520    }
1521    *seed = ctx.seed;
1522
1523    RESTORE_STACK;
1524 }
1525