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