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