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