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