Cleaning up leftovers of "freq" in celt_decode_with_ec()
[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       depth = (1+pulses[i])/((m->eBands[i+1]-m->eBands[i])<<LM);
285
286 #ifdef FIXED_POINT
287       thresh32 = SHR32(celt_exp2(-SHL16(depth, 10-BITRES)),1);
288       thresh = MULT16_32_Q15(QCONST16(0.5f, 15), MIN32(32767,thresh32));
289       {
290          opus_val32 t;
291          t = N0<<LM;
292          shift = celt_ilog2(t)>>1;
293          t = SHL32(t, (7-shift)<<1);
294          sqrt_1 = celt_rsqrt_norm(t);
295       }
296 #else
297       thresh = .5f*celt_exp2(-.125f*depth);
298       sqrt_1 = celt_rsqrt(N0<<LM);
299 #endif
300
301       c=0; do
302       {
303          celt_norm *X;
304          opus_val16 prev1;
305          opus_val16 prev2;
306          opus_val32 Ediff;
307          opus_val16 r;
308          int renormalize=0;
309          prev1 = prev1logE[c*m->nbEBands+i];
310          prev2 = prev2logE[c*m->nbEBands+i];
311          if (C==1)
312          {
313             prev1 = MAX16(prev1,prev1logE[m->nbEBands+i]);
314             prev2 = MAX16(prev2,prev2logE[m->nbEBands+i]);
315          }
316          Ediff = EXTEND32(logE[c*m->nbEBands+i])-EXTEND32(MIN16(prev1,prev2));
317          Ediff = MAX32(0, Ediff);
318
319 #ifdef FIXED_POINT
320          if (Ediff < 16384)
321          {
322             opus_val32 r32 = SHR32(celt_exp2(-EXTRACT16(Ediff)),1);
323             r = 2*MIN16(16383,r32);
324          } else {
325             r = 0;
326          }
327          if (LM==3)
328             r = MULT16_16_Q14(23170, MIN32(23169, r));
329          r = SHR16(MIN16(thresh, r),1);
330          r = SHR32(MULT16_16_Q15(sqrt_1, r),shift);
331 #else
332          /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because
333             short blocks don't have the same energy as long */
334          r = 2.f*celt_exp2(-Ediff);
335          if (LM==3)
336             r *= 1.41421356f;
337          r = MIN16(thresh, r);
338          r = r*sqrt_1;
339 #endif
340          X = X_+c*size+(m->eBands[i]<<LM);
341          for (k=0;k<1<<LM;k++)
342          {
343             /* Detect collapse */
344             if (!(collapse_masks[i*C+c]&1<<k))
345             {
346                /* Fill with noise */
347                for (j=0;j<N0;j++)
348                {
349                   seed = celt_lcg_rand(seed);
350                   X[(j<<LM)+k] = (seed&0x8000 ? r : -r);
351                }
352                renormalize = 1;
353             }
354          }
355          /* We just added some energy, so we need to renormalise */
356          if (renormalize)
357             renormalise_vector(X, N0<<LM, Q15ONE);
358       } while (++c<C);
359    }
360 }
361
362 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)
363 {
364    int i = bandID;
365    int j;
366    opus_val16 a1, a2;
367    opus_val16 left, right;
368    opus_val16 norm;
369 #ifdef FIXED_POINT
370    int shift = celt_zlog2(MAX32(bandE[i], bandE[i+m->nbEBands]))-13;
371 #endif
372    left = VSHR32(bandE[i],shift);
373    right = VSHR32(bandE[i+m->nbEBands],shift);
374    norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
375    a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
376    a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);
377    for (j=0;j<N;j++)
378    {
379       celt_norm r, l;
380       l = X[j];
381       r = Y[j];
382       X[j] = EXTRACT16(SHR32(MAC16_16(MULT16_16(a1, l), a2, r), 14));
383       /* Side is not encoded, no need to calculate */
384    }
385 }
386
387 static void stereo_split(celt_norm * OPUS_RESTRICT X, celt_norm * OPUS_RESTRICT Y, int N)
388 {
389    int j;
390    for (j=0;j<N;j++)
391    {
392       opus_val32 r, l;
393       l = MULT16_16(QCONST16(.70710678f, 15), X[j]);
394       r = MULT16_16(QCONST16(.70710678f, 15), Y[j]);
395       X[j] = EXTRACT16(SHR32(ADD32(l, r), 15));
396       Y[j] = EXTRACT16(SHR32(SUB32(r, l), 15));
397    }
398 }
399
400 static void stereo_merge(celt_norm * OPUS_RESTRICT X, celt_norm * OPUS_RESTRICT Y, opus_val16 mid, int N)
401 {
402    int j;
403    opus_val32 xp=0, side=0;
404    opus_val32 El, Er;
405    opus_val16 mid2;
406 #ifdef FIXED_POINT
407    int kl, kr;
408 #endif
409    opus_val32 t, lgain, rgain;
410
411    /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
412    dual_inner_prod(Y, X, Y, N, &xp, &side);
413    /* Compensating for the mid normalization */
414    xp = MULT16_32_Q15(mid, xp);
415    /* mid and side are in Q15, not Q14 like X and Y */
416    mid2 = SHR32(mid, 1);
417    El = MULT16_16(mid2, mid2) + side - 2*xp;
418    Er = MULT16_16(mid2, mid2) + side + 2*xp;
419    if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28))
420    {
421       OPUS_COPY(Y, X, N);
422       return;
423    }
424
425 #ifdef FIXED_POINT
426    kl = celt_ilog2(El)>>1;
427    kr = celt_ilog2(Er)>>1;
428 #endif
429    t = VSHR32(El, (kl-7)<<1);
430    lgain = celt_rsqrt_norm(t);
431    t = VSHR32(Er, (kr-7)<<1);
432    rgain = celt_rsqrt_norm(t);
433
434 #ifdef FIXED_POINT
435    if (kl < 7)
436       kl = 7;
437    if (kr < 7)
438       kr = 7;
439 #endif
440
441    for (j=0;j<N;j++)
442    {
443       celt_norm r, l;
444       /* Apply mid scaling (side is already scaled) */
445       l = MULT16_16_P15(mid, X[j]);
446       r = Y[j];
447       X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));
448       Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));
449    }
450 }
451
452 /* Decide whether we should spread the pulses in the current frame */
453 int spreading_decision(const CELTMode *m, const celt_norm *X, int *average,
454       int last_decision, int *hf_average, int *tapset_decision, int update_hf,
455       int end, int C, int M)
456 {
457    int i, c, N0;
458    int sum = 0, nbBands=0;
459    const opus_int16 * OPUS_RESTRICT eBands = m->eBands;
460    int decision;
461    int hf_sum=0;
462
463    celt_assert(end>0);
464
465    N0 = M*m->shortMdctSize;
466
467    if (M*(eBands[end]-eBands[end-1]) <= 8)
468       return SPREAD_NONE;
469    c=0; do {
470       for (i=0;i<end;i++)
471       {
472          int j, N, tmp=0;
473          int tcount[3] = {0,0,0};
474          const celt_norm * OPUS_RESTRICT x = X+M*eBands[i]+c*N0;
475          N = M*(eBands[i+1]-eBands[i]);
476          if (N<=8)
477             continue;
478          /* Compute rough CDF of |x[j]| */
479          for (j=0;j<N;j++)
480          {
481             opus_val32 x2N; /* Q13 */
482
483             x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
484             if (x2N < QCONST16(0.25f,13))
485                tcount[0]++;
486             if (x2N < QCONST16(0.0625f,13))
487                tcount[1]++;
488             if (x2N < QCONST16(0.015625f,13))
489                tcount[2]++;
490          }
491
492          /* Only include four last bands (8 kHz and up) */
493          if (i>m->nbEBands-4)
494             hf_sum += 32*(tcount[1]+tcount[0])/N;
495          tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
496          sum += tmp*256;
497          nbBands++;
498       }
499    } while (++c<C);
500
501    if (update_hf)
502    {
503       if (hf_sum)
504          hf_sum /= C*(4-m->nbEBands+end);
505       *hf_average = (*hf_average+hf_sum)>>1;
506       hf_sum = *hf_average;
507       if (*tapset_decision==2)
508          hf_sum += 4;
509       else if (*tapset_decision==0)
510          hf_sum -= 4;
511       if (hf_sum > 22)
512          *tapset_decision=2;
513       else if (hf_sum > 18)
514          *tapset_decision=1;
515       else
516          *tapset_decision=0;
517    }
518    /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/
519    celt_assert(nbBands>0); /* end has to be non-zero */
520    sum /= nbBands;
521    /* Recursive averaging */
522    sum = (sum+*average)>>1;
523    *average = sum;
524    /* Hysteresis */
525    sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2;
526    if (sum < 80)
527    {
528       decision = SPREAD_AGGRESSIVE;
529    } else if (sum < 256)
530    {
531       decision = SPREAD_NORMAL;
532    } else if (sum < 384)
533    {
534       decision = SPREAD_LIGHT;
535    } else {
536       decision = SPREAD_NONE;
537    }
538 #ifdef FUZZING
539    decision = rand()&0x3;
540    *tapset_decision=rand()%3;
541 #endif
542    return decision;
543 }
544
545 /* Indexing table for converting from natural Hadamard to ordery Hadamard
546    This is essentially a bit-reversed Gray, on top of which we've added
547    an inversion of the order because we want the DC at the end rather than
548    the beginning. The lines are for N=2, 4, 8, 16 */
549 static const int ordery_table[] = {
550        1,  0,
551        3,  0,  2,  1,
552        7,  0,  4,  3,  6,  1,  5,  2,
553       15,  0,  8,  7, 12,  3, 11,  4, 14,  1,  9,  6, 13,  2, 10,  5,
554 };
555
556 static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
557 {
558    int i,j;
559    VARDECL(celt_norm, tmp);
560    int N;
561    SAVE_STACK;
562    N = N0*stride;
563    ALLOC(tmp, N, celt_norm);
564    celt_assert(stride>0);
565    if (hadamard)
566    {
567       const int *ordery = ordery_table+stride-2;
568       for (i=0;i<stride;i++)
569       {
570          for (j=0;j<N0;j++)
571             tmp[ordery[i]*N0+j] = X[j*stride+i];
572       }
573    } else {
574       for (i=0;i<stride;i++)
575          for (j=0;j<N0;j++)
576             tmp[i*N0+j] = X[j*stride+i];
577    }
578    OPUS_COPY(X, tmp, N);
579    RESTORE_STACK;
580 }
581
582 static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
583 {
584    int i,j;
585    VARDECL(celt_norm, tmp);
586    int N;
587    SAVE_STACK;
588    N = N0*stride;
589    ALLOC(tmp, N, celt_norm);
590    if (hadamard)
591    {
592       const int *ordery = ordery_table+stride-2;
593       for (i=0;i<stride;i++)
594          for (j=0;j<N0;j++)
595             tmp[j*stride+i] = X[ordery[i]*N0+j];
596    } else {
597       for (i=0;i<stride;i++)
598          for (j=0;j<N0;j++)
599             tmp[j*stride+i] = X[i*N0+j];
600    }
601    OPUS_COPY(X, tmp, N);
602    RESTORE_STACK;
603 }
604
605 void haar1(celt_norm *X, int N0, int stride)
606 {
607    int i, j;
608    N0 >>= 1;
609    for (i=0;i<stride;i++)
610       for (j=0;j<N0;j++)
611       {
612          opus_val32 tmp1, tmp2;
613          tmp1 = MULT16_16(QCONST16(.70710678f,15), X[stride*2*j+i]);
614          tmp2 = MULT16_16(QCONST16(.70710678f,15), X[stride*(2*j+1)+i]);
615          X[stride*2*j+i] = EXTRACT16(PSHR32(ADD32(tmp1, tmp2), 15));
616          X[stride*(2*j+1)+i] = EXTRACT16(PSHR32(SUB32(tmp1, tmp2), 15));
617       }
618 }
619
620 static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo)
621 {
622    static const opus_int16 exp2_table8[8] =
623       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
624    int qn, qb;
625    int N2 = 2*N-1;
626    if (stereo && N==2)
627       N2--;
628    /* The upper limit ensures that in a stereo split with itheta==16384, we'll
629        always have enough bits left over to code at least one pulse in the
630        side; otherwise it would collapse, since it doesn't get folded. */
631    qb = IMIN(b-pulse_cap-(4<<BITRES), (b+N2*offset)/N2);
632
633    qb = IMIN(8<<BITRES, qb);
634
635    if (qb<(1<<BITRES>>1)) {
636       qn = 1;
637    } else {
638       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
639       qn = (qn+1)>>1<<1;
640    }
641    celt_assert(qn <= 256);
642    return qn;
643 }
644
645 struct band_ctx {
646    int encode;
647    const CELTMode *m;
648    int i;
649    int intensity;
650    int spread;
651    int tf_change;
652    ec_ctx *ec;
653    opus_int32 remaining_bits;
654    const celt_ener *bandE;
655    opus_uint32 seed;
656 };
657
658 struct split_ctx {
659    int inv;
660    int imid;
661    int iside;
662    int delta;
663    int itheta;
664    int qalloc;
665 };
666
667 static void compute_theta(struct band_ctx *ctx, struct split_ctx *sctx,
668       celt_norm *X, celt_norm *Y, int N, int *b, int B, int B0,
669       int LM,
670       int stereo, int *fill)
671 {
672    int qn;
673    int itheta=0;
674    int delta;
675    int imid, iside;
676    int qalloc;
677    int pulse_cap;
678    int offset;
679    opus_int32 tell;
680    int inv=0;
681    int encode;
682    const CELTMode *m;
683    int i;
684    int intensity;
685    ec_ctx *ec;
686    const celt_ener *bandE;
687
688    encode = ctx->encode;
689    m = ctx->m;
690    i = ctx->i;
691    intensity = ctx->intensity;
692    ec = ctx->ec;
693    bandE = ctx->bandE;
694
695    /* Decide on the resolution to give to the split parameter theta */
696    pulse_cap = m->logN[i]+LM*(1<<BITRES);
697    offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
698    qn = compute_qn(N, *b, offset, pulse_cap, stereo);
699    if (stereo && i>=intensity)
700       qn = 1;
701    if (encode)
702    {
703       /* theta is the atan() of the ratio between the (normalized)
704          side and mid. With just that parameter, we can re-scale both
705          mid and side because we know that 1) they have unit norm and
706          2) they are orthogonal. */
707       itheta = stereo_itheta(X, Y, stereo, N);
708    }
709    tell = ec_tell_frac(ec);
710    if (qn!=1)
711    {
712       if (encode)
713          itheta = (itheta*qn+8192)>>14;
714
715       /* Entropy coding of the angle. We use a uniform pdf for the
716          time split, a step for stereo, and a triangular one for the rest. */
717       if (stereo && N>2)
718       {
719          int p0 = 3;
720          int x = itheta;
721          int x0 = qn/2;
722          int ft = p0*(x0+1) + x0;
723          /* Use a probability of p0 up to itheta=8192 and then use 1 after */
724          if (encode)
725          {
726             ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
727          } else {
728             int fs;
729             fs=ec_decode(ec,ft);
730             if (fs<(x0+1)*p0)
731                x=fs/p0;
732             else
733                x=x0+1+(fs-(x0+1)*p0);
734             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);
735             itheta = x;
736          }
737       } else if (B0>1 || stereo) {
738          /* Uniform pdf */
739          if (encode)
740             ec_enc_uint(ec, itheta, qn+1);
741          else
742             itheta = ec_dec_uint(ec, qn+1);
743       } else {
744          int fs=1, ft;
745          ft = ((qn>>1)+1)*((qn>>1)+1);
746          if (encode)
747          {
748             int fl;
749
750             fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
751             fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
752              ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
753
754             ec_encode(ec, fl, fl+fs, ft);
755          } else {
756             /* Triangular pdf */
757             int fl=0;
758             int fm;
759             fm = ec_decode(ec, ft);
760
761             if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
762             {
763                itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;
764                fs = itheta + 1;
765                fl = itheta*(itheta + 1)>>1;
766             }
767             else
768             {
769                itheta = (2*(qn + 1)
770                 - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;
771                fs = qn + 1 - itheta;
772                fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
773             }
774
775             ec_dec_update(ec, fl, fl+fs, ft);
776          }
777       }
778       itheta = (opus_int32)itheta*16384/qn;
779       if (encode && stereo)
780       {
781          if (itheta==0)
782             intensity_stereo(m, X, Y, bandE, i, N);
783          else
784             stereo_split(X, Y, N);
785       }
786       /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.
787                Let's do that at higher complexity */
788    } else if (stereo) {
789       if (encode)
790       {
791          inv = itheta > 8192;
792          if (inv)
793          {
794             int j;
795             for (j=0;j<N;j++)
796                Y[j] = -Y[j];
797          }
798          intensity_stereo(m, X, Y, bandE, i, N);
799       }
800       if (*b>2<<BITRES && ctx->remaining_bits > 2<<BITRES)
801       {
802          if (encode)
803             ec_enc_bit_logp(ec, inv, 2);
804          else
805             inv = ec_dec_bit_logp(ec, 2);
806       } else
807          inv = 0;
808       itheta = 0;
809    }
810    qalloc = ec_tell_frac(ec) - tell;
811    *b -= qalloc;
812
813    if (itheta == 0)
814    {
815       imid = 32767;
816       iside = 0;
817       *fill &= (1<<B)-1;
818       delta = -16384;
819    } else if (itheta == 16384)
820    {
821       imid = 0;
822       iside = 32767;
823       *fill &= ((1<<B)-1)<<B;
824       delta = 16384;
825    } else {
826       imid = bitexact_cos((opus_int16)itheta);
827       iside = bitexact_cos((opus_int16)(16384-itheta));
828       /* This is the mid vs side allocation that minimizes squared error
829          in that band. */
830       delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid));
831    }
832
833    sctx->inv = inv;
834    sctx->imid = imid;
835    sctx->iside = iside;
836    sctx->delta = delta;
837    sctx->itheta = itheta;
838    sctx->qalloc = qalloc;
839 }
840 static unsigned quant_band_n1(struct band_ctx *ctx, celt_norm *X, celt_norm *Y, int b,
841       celt_norm *lowband_out)
842 {
843 #ifdef RESYNTH
844    int resynth = 1;
845 #else
846    int resynth = !ctx->encode;
847 #endif
848    int c;
849    int stereo;
850    celt_norm *x = X;
851    int encode;
852    ec_ctx *ec;
853
854    encode = ctx->encode;
855    ec = ctx->ec;
856
857    stereo = Y != NULL;
858    c=0; do {
859       int sign=0;
860       if (ctx->remaining_bits>=1<<BITRES)
861       {
862          if (encode)
863          {
864             sign = x[0]<0;
865             ec_enc_bits(ec, sign, 1);
866          } else {
867             sign = ec_dec_bits(ec, 1);
868          }
869          ctx->remaining_bits -= 1<<BITRES;
870          b-=1<<BITRES;
871       }
872       if (resynth)
873          x[0] = sign ? -NORM_SCALING : NORM_SCALING;
874       x = Y;
875    } while (++c<1+stereo);
876    if (lowband_out)
877       lowband_out[0] = SHR16(X[0],4);
878    return 1;
879 }
880
881 /* This function is responsible for encoding and decoding a mono partition.
882    It can split the band in two and transmit the energy difference with
883    the two half-bands. It can be called recursively so bands can end up being
884    split in 8 parts. */
885 static unsigned quant_partition(struct band_ctx *ctx, celt_norm *X,
886       int N, int b, int B, celt_norm *lowband,
887       int LM,
888       opus_val16 gain, int fill)
889 {
890    const unsigned char *cache;
891    int q;
892    int curr_bits;
893    int imid=0, iside=0;
894    int B0=B;
895    opus_val16 mid=0, side=0;
896    unsigned cm=0;
897 #ifdef RESYNTH
898    int resynth = 1;
899 #else
900    int resynth = !ctx->encode;
901 #endif
902    celt_norm *Y=NULL;
903    int encode;
904    const CELTMode *m;
905    int i;
906    int spread;
907    ec_ctx *ec;
908
909    encode = ctx->encode;
910    m = ctx->m;
911    i = ctx->i;
912    spread = ctx->spread;
913    ec = ctx->ec;
914
915    /* If we need 1.5 more bit than we can produce, split the band in two. */
916    cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
917    if (LM != -1 && b > cache[cache[0]]+12 && N>2)
918    {
919       int mbits, sbits, delta;
920       int itheta;
921       int qalloc;
922       struct split_ctx sctx;
923       celt_norm *next_lowband2=NULL;
924       opus_int32 rebalance;
925
926       N >>= 1;
927       Y = X+N;
928       LM -= 1;
929       if (B==1)
930          fill = (fill&1)|(fill<<1);
931       B = (B+1)>>1;
932
933       compute_theta(ctx, &sctx, X, Y, N, &b, B, B0,
934             LM, 0, &fill);
935       imid = sctx.imid;
936       iside = sctx.iside;
937       delta = sctx.delta;
938       itheta = sctx.itheta;
939       qalloc = sctx.qalloc;
940 #ifdef FIXED_POINT
941       mid = imid;
942       side = iside;
943 #else
944       mid = (1.f/32768)*imid;
945       side = (1.f/32768)*iside;
946 #endif
947
948       /* Give more bits to low-energy MDCTs than they would otherwise deserve */
949       if (B0>1 && (itheta&0x3fff))
950       {
951          if (itheta > 8192)
952             /* Rough approximation for pre-echo masking */
953             delta -= delta>>(4-LM);
954          else
955             /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
956             delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
957       }
958       mbits = IMAX(0, IMIN(b, (b-delta)/2));
959       sbits = b-mbits;
960       ctx->remaining_bits -= qalloc;
961
962       if (lowband)
963          next_lowband2 = lowband+N; /* >32-bit split case */
964
965       rebalance = ctx->remaining_bits;
966       if (mbits >= sbits)
967       {
968          cm = quant_partition(ctx, X, N, mbits, B,
969                lowband, LM,
970                MULT16_16_P15(gain,mid), fill);
971          rebalance = mbits - (rebalance-ctx->remaining_bits);
972          if (rebalance > 3<<BITRES && itheta!=0)
973             sbits += rebalance - (3<<BITRES);
974          cm |= quant_partition(ctx, Y, N, sbits, B,
975                next_lowband2, LM,
976                MULT16_16_P15(gain,side), fill>>B)<<(B0>>1);
977       } else {
978          cm = quant_partition(ctx, Y, N, sbits, B,
979                next_lowband2, LM,
980                MULT16_16_P15(gain,side), fill>>B)<<(B0>>1);
981          rebalance = sbits - (rebalance-ctx->remaining_bits);
982          if (rebalance > 3<<BITRES && itheta!=16384)
983             mbits += rebalance - (3<<BITRES);
984          cm |= quant_partition(ctx, X, N, mbits, B,
985                lowband, LM,
986                MULT16_16_P15(gain,mid), fill);
987       }
988    } else {
989       /* This is the basic no-split case */
990       q = bits2pulses(m, i, LM, b);
991       curr_bits = pulses2bits(m, i, LM, q);
992       ctx->remaining_bits -= curr_bits;
993
994       /* Ensures we can never bust the budget */
995       while (ctx->remaining_bits < 0 && q > 0)
996       {
997          ctx->remaining_bits += curr_bits;
998          q--;
999          curr_bits = pulses2bits(m, i, LM, q);
1000          ctx->remaining_bits -= curr_bits;
1001       }
1002
1003       if (q!=0)
1004       {
1005          int K = get_pulses(q);
1006
1007          /* Finally do the actual quantization */
1008          if (encode)
1009          {
1010             cm = alg_quant(X, N, K, spread, B, ec
1011 #ifdef RESYNTH
1012                  , gain
1013 #endif
1014                  );
1015          } else {
1016             cm = alg_unquant(X, N, K, spread, B, ec, gain);
1017          }
1018       } else {
1019          /* If there's no pulse, fill the band anyway */
1020          int j;
1021          if (resynth)
1022          {
1023             unsigned cm_mask;
1024             /* B can be as large as 16, so this shift might overflow an int on a
1025                16-bit platform; use a long to get defined behavior.*/
1026             cm_mask = (unsigned)(1UL<<B)-1;
1027             fill &= cm_mask;
1028             if (!fill)
1029             {
1030                OPUS_CLEAR(X, N);
1031             } else {
1032                if (lowband == NULL)
1033                {
1034                   /* Noise */
1035                   for (j=0;j<N;j++)
1036                   {
1037                      ctx->seed = celt_lcg_rand(ctx->seed);
1038                      X[j] = (celt_norm)((opus_int32)ctx->seed>>20);
1039                   }
1040                   cm = cm_mask;
1041                } else {
1042                   /* Folded spectrum */
1043                   for (j=0;j<N;j++)
1044                   {
1045                      opus_val16 tmp;
1046                      ctx->seed = celt_lcg_rand(ctx->seed);
1047                      /* About 48 dB below the "normal" folding level */
1048                      tmp = QCONST16(1.0f/256, 10);
1049                      tmp = (ctx->seed)&0x8000 ? tmp : -tmp;
1050                      X[j] = lowband[j]+tmp;
1051                   }
1052                   cm = fill;
1053                }
1054                renormalise_vector(X, N, gain);
1055             }
1056          }
1057       }
1058    }
1059
1060    return cm;
1061 }
1062
1063
1064 /* This function is responsible for encoding and decoding a band for the mono case. */
1065 static unsigned quant_band(struct band_ctx *ctx, celt_norm *X,
1066       int N, int b, int B, celt_norm *lowband,
1067       int LM, celt_norm *lowband_out,
1068       opus_val16 gain, celt_norm *lowband_scratch, int fill)
1069 {
1070    int N0=N;
1071    int N_B=N;
1072    int N_B0;
1073    int B0=B;
1074    int time_divide=0;
1075    int recombine=0;
1076    int longBlocks;
1077    unsigned cm=0;
1078 #ifdef RESYNTH
1079    int resynth = 1;
1080 #else
1081    int resynth = !ctx->encode;
1082 #endif
1083    int k;
1084    int encode;
1085    int tf_change;
1086
1087    encode = ctx->encode;
1088    tf_change = ctx->tf_change;
1089
1090    longBlocks = B0==1;
1091
1092    N_B /= B;
1093
1094    /* Special case for one sample */
1095    if (N==1)
1096    {
1097       return quant_band_n1(ctx, X, NULL, b, lowband_out);
1098    }
1099
1100    if (tf_change>0)
1101       recombine = tf_change;
1102    /* Band recombining to increase frequency resolution */
1103
1104    if (lowband_scratch && lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
1105    {
1106       OPUS_COPY(lowband_scratch, lowband, N);
1107       lowband = lowband_scratch;
1108    }
1109
1110    for (k=0;k<recombine;k++)
1111    {
1112       static const unsigned char bit_interleave_table[16]={
1113             0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
1114       };
1115       if (encode)
1116          haar1(X, N>>k, 1<<k);
1117       if (lowband)
1118          haar1(lowband, N>>k, 1<<k);
1119       fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
1120    }
1121    B>>=recombine;
1122    N_B<<=recombine;
1123
1124    /* Increasing the time resolution */
1125    while ((N_B&1) == 0 && tf_change<0)
1126    {
1127       if (encode)
1128          haar1(X, N_B, B);
1129       if (lowband)
1130          haar1(lowband, N_B, B);
1131       fill |= fill<<B;
1132       B <<= 1;
1133       N_B >>= 1;
1134       time_divide++;
1135       tf_change++;
1136    }
1137    B0=B;
1138    N_B0 = N_B;
1139
1140    /* Reorganize the samples in time order instead of frequency order */
1141    if (B0>1)
1142    {
1143       if (encode)
1144          deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1145       if (lowband)
1146          deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
1147    }
1148
1149    cm = quant_partition(ctx, X, N, b, B, lowband,
1150          LM, gain, fill);
1151
1152    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1153    if (resynth)
1154    {
1155       /* Undo the sample reorganization going from time order to frequency order */
1156       if (B0>1)
1157          interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1158
1159       /* Undo time-freq changes that we did earlier */
1160       N_B = N_B0;
1161       B = B0;
1162       for (k=0;k<time_divide;k++)
1163       {
1164          B >>= 1;
1165          N_B <<= 1;
1166          cm |= cm>>B;
1167          haar1(X, N_B, B);
1168       }
1169
1170       for (k=0;k<recombine;k++)
1171       {
1172          static const unsigned char bit_deinterleave_table[16]={
1173                0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1174                0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1175          };
1176          cm = bit_deinterleave_table[cm];
1177          haar1(X, N0>>k, 1<<k);
1178       }
1179       B<<=recombine;
1180
1181       /* Scale output for later folding */
1182       if (lowband_out)
1183       {
1184          int j;
1185          opus_val16 n;
1186          n = celt_sqrt(SHL32(EXTEND32(N0),22));
1187          for (j=0;j<N0;j++)
1188             lowband_out[j] = MULT16_16_Q15(n,X[j]);
1189       }
1190       cm &= (1<<B)-1;
1191    }
1192    return cm;
1193 }
1194
1195
1196 /* This function is responsible for encoding and decoding a band for the stereo case. */
1197 static unsigned quant_band_stereo(struct band_ctx *ctx, celt_norm *X, celt_norm *Y,
1198       int N, int b, int B, celt_norm *lowband,
1199       int LM, celt_norm *lowband_out,
1200       celt_norm *lowband_scratch, int fill)
1201 {
1202    int imid=0, iside=0;
1203    int inv = 0;
1204    opus_val16 mid=0, side=0;
1205    unsigned cm=0;
1206 #ifdef RESYNTH
1207    int resynth = 1;
1208 #else
1209    int resynth = !ctx->encode;
1210 #endif
1211    int mbits, sbits, delta;
1212    int itheta;
1213    int qalloc;
1214    struct split_ctx sctx;
1215    int orig_fill;
1216    int encode;
1217    ec_ctx *ec;
1218
1219    encode = ctx->encode;
1220    ec = ctx->ec;
1221
1222    /* Special case for one sample */
1223    if (N==1)
1224    {
1225       return quant_band_n1(ctx, X, Y, b, lowband_out);
1226    }
1227
1228    orig_fill = fill;
1229
1230    compute_theta(ctx, &sctx, X, Y, N, &b, B, B,
1231          LM, 1, &fill);
1232    inv = sctx.inv;
1233    imid = sctx.imid;
1234    iside = sctx.iside;
1235    delta = sctx.delta;
1236    itheta = sctx.itheta;
1237    qalloc = sctx.qalloc;
1238 #ifdef FIXED_POINT
1239    mid = imid;
1240    side = iside;
1241 #else
1242    mid = (1.f/32768)*imid;
1243    side = (1.f/32768)*iside;
1244 #endif
1245
1246    /* This is a special case for N=2 that only works for stereo and takes
1247       advantage of the fact that mid and side are orthogonal to encode
1248       the side with just one bit. */
1249    if (N==2)
1250    {
1251       int c;
1252       int sign=0;
1253       celt_norm *x2, *y2;
1254       mbits = b;
1255       sbits = 0;
1256       /* Only need one bit for the side. */
1257       if (itheta != 0 && itheta != 16384)
1258          sbits = 1<<BITRES;
1259       mbits -= sbits;
1260       c = itheta > 8192;
1261       ctx->remaining_bits -= qalloc+sbits;
1262
1263       x2 = c ? Y : X;
1264       y2 = c ? X : Y;
1265       if (sbits)
1266       {
1267          if (encode)
1268          {
1269             /* Here we only need to encode a sign for the side. */
1270             sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
1271             ec_enc_bits(ec, sign, 1);
1272          } else {
1273             sign = ec_dec_bits(ec, 1);
1274          }
1275       }
1276       sign = 1-2*sign;
1277       /* We use orig_fill here because we want to fold the side, but if
1278          itheta==16384, we'll have cleared the low bits of fill. */
1279       cm = quant_band(ctx, x2, N, mbits, B, lowband,
1280             LM, lowband_out, Q15ONE, lowband_scratch, orig_fill);
1281       /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
1282          and there's no need to worry about mixing with the other channel. */
1283       y2[0] = -sign*x2[1];
1284       y2[1] = sign*x2[0];
1285       if (resynth)
1286       {
1287          celt_norm tmp;
1288          X[0] = MULT16_16_Q15(mid, X[0]);
1289          X[1] = MULT16_16_Q15(mid, X[1]);
1290          Y[0] = MULT16_16_Q15(side, Y[0]);
1291          Y[1] = MULT16_16_Q15(side, Y[1]);
1292          tmp = X[0];
1293          X[0] = SUB16(tmp,Y[0]);
1294          Y[0] = ADD16(tmp,Y[0]);
1295          tmp = X[1];
1296          X[1] = SUB16(tmp,Y[1]);
1297          Y[1] = ADD16(tmp,Y[1]);
1298       }
1299    } else {
1300       /* "Normal" split code */
1301       opus_int32 rebalance;
1302
1303       mbits = IMAX(0, IMIN(b, (b-delta)/2));
1304       sbits = b-mbits;
1305       ctx->remaining_bits -= qalloc;
1306
1307       rebalance = ctx->remaining_bits;
1308       if (mbits >= sbits)
1309       {
1310          /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1311             mid for folding later. */
1312          cm = quant_band(ctx, X, N, mbits, B,
1313                lowband, LM, lowband_out,
1314                Q15ONE, lowband_scratch, fill);
1315          rebalance = mbits - (rebalance-ctx->remaining_bits);
1316          if (rebalance > 3<<BITRES && itheta!=0)
1317             sbits += rebalance - (3<<BITRES);
1318
1319          /* For a stereo split, the high bits of fill are always zero, so no
1320             folding will be done to the side. */
1321          cm |= quant_band(ctx, Y, N, sbits, B,
1322                NULL, LM, NULL,
1323                side, NULL, fill>>B);
1324       } else {
1325          /* For a stereo split, the high bits of fill are always zero, so no
1326             folding will be done to the side. */
1327          cm = quant_band(ctx, Y, N, sbits, B,
1328                NULL, LM, NULL,
1329                side, NULL, fill>>B);
1330          rebalance = sbits - (rebalance-ctx->remaining_bits);
1331          if (rebalance > 3<<BITRES && itheta!=16384)
1332             mbits += rebalance - (3<<BITRES);
1333          /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1334             mid for folding later. */
1335          cm |= quant_band(ctx, X, N, mbits, B,
1336                lowband, LM, lowband_out,
1337                Q15ONE, lowband_scratch, fill);
1338       }
1339    }
1340
1341
1342    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1343    if (resynth)
1344    {
1345       if (N!=2)
1346          stereo_merge(X, Y, mid, N);
1347       if (inv)
1348       {
1349          int j;
1350          for (j=0;j<N;j++)
1351             Y[j] = -Y[j];
1352       }
1353    }
1354    return cm;
1355 }
1356
1357
1358 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1359       celt_norm *X_, celt_norm *Y_, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1360       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res,
1361       opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed)
1362 {
1363    int i;
1364    opus_int32 remaining_bits;
1365    const opus_int16 * OPUS_RESTRICT eBands = m->eBands;
1366    celt_norm * OPUS_RESTRICT norm, * OPUS_RESTRICT norm2;
1367    VARDECL(celt_norm, _norm);
1368    celt_norm *lowband_scratch;
1369    int B;
1370    int M;
1371    int lowband_offset;
1372    int update_lowband = 1;
1373    int C = Y_ != NULL ? 2 : 1;
1374    int norm_offset;
1375 #ifdef RESYNTH
1376    int resynth = 1;
1377 #else
1378    int resynth = !encode;
1379 #endif
1380    struct band_ctx ctx;
1381    SAVE_STACK;
1382
1383    M = 1<<LM;
1384    B = shortBlocks ? M : 1;
1385    norm_offset = M*eBands[start];
1386    /* No need to allocate norm for the last band because we don't need an
1387       output in that band. */
1388    ALLOC(_norm, C*(M*eBands[m->nbEBands-1]-norm_offset), celt_norm);
1389    norm = _norm;
1390    norm2 = norm + M*eBands[m->nbEBands-1]-norm_offset;
1391    /* We can use the last band as scratch space because we don't need that
1392       scratch space for the last band. */
1393    lowband_scratch = X_+M*eBands[m->nbEBands-1];
1394
1395    lowband_offset = 0;
1396    ctx.bandE = bandE;
1397    ctx.ec = ec;
1398    ctx.encode = encode;
1399    ctx.intensity = intensity;
1400    ctx.m = m;
1401    ctx.seed = *seed;
1402    ctx.spread = spread;
1403    for (i=start;i<end;i++)
1404    {
1405       opus_int32 tell;
1406       int b;
1407       int N;
1408       opus_int32 curr_balance;
1409       int effective_lowband=-1;
1410       celt_norm * OPUS_RESTRICT X, * OPUS_RESTRICT Y;
1411       int tf_change=0;
1412       unsigned x_cm;
1413       unsigned y_cm;
1414       int last;
1415
1416       ctx.i = i;
1417       last = (i==end-1);
1418
1419       X = X_+M*eBands[i];
1420       if (Y_!=NULL)
1421          Y = Y_+M*eBands[i];
1422       else
1423          Y = NULL;
1424       N = M*eBands[i+1]-M*eBands[i];
1425       tell = ec_tell_frac(ec);
1426
1427       /* Compute how many bits we want to allocate to this band */
1428       if (i != start)
1429          balance -= tell;
1430       remaining_bits = total_bits-tell-1;
1431       ctx.remaining_bits = remaining_bits;
1432       if (i <= codedBands-1)
1433       {
1434          curr_balance = balance / IMIN(3, codedBands-i);
1435          b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1436       } else {
1437          b = 0;
1438       }
1439
1440       if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1441             lowband_offset = i;
1442
1443       tf_change = tf_res[i];
1444       ctx.tf_change = tf_change;
1445       if (i>=m->effEBands)
1446       {
1447          X=norm;
1448          if (Y_!=NULL)
1449             Y = norm;
1450          lowband_scratch = NULL;
1451       }
1452       if (i==end-1)
1453          lowband_scratch = NULL;
1454
1455       /* Get a conservative estimate of the collapse_mask's for the bands we're
1456          going to be folding from. */
1457       if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1458       {
1459          int fold_start;
1460          int fold_end;
1461          int fold_i;
1462          /* This ensures we never repeat spectral content within one band */
1463          effective_lowband = IMAX(0, M*eBands[lowband_offset]-norm_offset-N);
1464          fold_start = lowband_offset;
1465          while(M*eBands[--fold_start] > effective_lowband+norm_offset);
1466          fold_end = lowband_offset-1;
1467          while(M*eBands[++fold_end] < effective_lowband+norm_offset+N);
1468          x_cm = y_cm = 0;
1469          fold_i = fold_start; do {
1470            x_cm |= collapse_masks[fold_i*C+0];
1471            y_cm |= collapse_masks[fold_i*C+C-1];
1472          } while (++fold_i<fold_end);
1473       }
1474       /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1475          always) be non-zero. */
1476       else
1477          x_cm = y_cm = (1<<B)-1;
1478
1479       if (dual_stereo && i==intensity)
1480       {
1481          int j;
1482
1483          /* Switch off dual stereo to do intensity. */
1484          dual_stereo = 0;
1485          if (resynth)
1486             for (j=0;j<M*eBands[i]-norm_offset;j++)
1487                norm[j] = HALF32(norm[j]+norm2[j]);
1488       }
1489       if (dual_stereo)
1490       {
1491          x_cm = quant_band(&ctx, X, N, b/2, B,
1492                effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1493                last?NULL:norm+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, x_cm);
1494          y_cm = quant_band(&ctx, Y, N, b/2, B,
1495                effective_lowband != -1 ? norm2+effective_lowband : NULL, LM,
1496                last?NULL:norm2+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, y_cm);
1497       } else {
1498          if (Y!=NULL)
1499          {
1500             x_cm = quant_band_stereo(&ctx, X, Y, N, b, B,
1501                   effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1502                         last?NULL:norm+M*eBands[i]-norm_offset, lowband_scratch, x_cm|y_cm);
1503          } else {
1504             x_cm = quant_band(&ctx, X, N, b, B,
1505                   effective_lowband != -1 ? norm+effective_lowband : NULL, LM,
1506                         last?NULL:norm+M*eBands[i]-norm_offset, Q15ONE, lowband_scratch, x_cm|y_cm);
1507          }
1508          y_cm = x_cm;
1509       }
1510       collapse_masks[i*C+0] = (unsigned char)x_cm;
1511       collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1512       balance += pulses[i] + tell;
1513
1514       /* Update the folding position only as long as we have 1 bit/sample depth. */
1515       update_lowband = b>(N<<BITRES);
1516    }
1517    *seed = ctx.seed;
1518
1519    RESTORE_STACK;
1520 }
1521