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