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