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