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