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