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