As of commit 44203907 all celt bands N>1 are even, including custom modes. In light...
[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    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    sum /= nbBands;
461    /* Recursive averaging */
462    sum = (sum+*average)>>1;
463    *average = sum;
464    /* Hysteresis */
465    sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2;
466    if (sum < 80)
467    {
468       decision = SPREAD_AGGRESSIVE;
469    } else if (sum < 256)
470    {
471       decision = SPREAD_NORMAL;
472    } else if (sum < 384)
473    {
474       decision = SPREAD_LIGHT;
475    } else {
476       decision = SPREAD_NONE;
477    }
478 #ifdef FUZZING
479    decision = rand()&0x3;
480    *tapset_decision=rand()%3;
481 #endif
482    return decision;
483 }
484
485 #ifdef MEASURE_NORM_MSE
486
487 float MSE[30] = {0};
488 int nbMSEBands = 0;
489 int MSECount[30] = {0};
490
491 void dump_norm_mse(void)
492 {
493    int i;
494    for (i=0;i<nbMSEBands;i++)
495    {
496       printf ("%g ", MSE[i]/MSECount[i]);
497    }
498    printf ("\n");
499 }
500
501 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
502 {
503    static int init = 0;
504    int i;
505    if (!init)
506    {
507       atexit(dump_norm_mse);
508       init = 1;
509    }
510    for (i=0;i<m->nbEBands;i++)
511    {
512       int j;
513       int c;
514       float g;
515       if (bandE0[i]<10 || (C==2 && bandE0[i+m->nbEBands]<1))
516          continue;
517       c=0; do {
518          g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
519          for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
520             MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
521       } while (++c<C);
522       MSECount[i]+=C;
523    }
524    nbMSEBands = m->nbEBands;
525 }
526
527 #endif
528
529 /* Indexing table for converting from natural Hadamard to ordery Hadamard
530    This is essentially a bit-reversed Gray, on top of which we've added
531    an inversion of the order because we want the DC at the end rather than
532    the beginning. The lines are for N=2, 4, 8, 16 */
533 static const int ordery_table[] = {
534        1,  0,
535        3,  0,  2,  1,
536        7,  0,  4,  3,  6,  1,  5,  2,
537       15,  0,  8,  7, 12,  3, 11,  4, 14,  1,  9,  6, 13,  2, 10,  5,
538 };
539
540 static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
541 {
542    int i,j;
543    VARDECL(celt_norm, tmp);
544    int N;
545    SAVE_STACK;
546    N = N0*stride;
547    ALLOC(tmp, N, celt_norm);
548    if (hadamard)
549    {
550       const int *ordery = ordery_table+stride-2;
551       for (i=0;i<stride;i++)
552       {
553          for (j=0;j<N0;j++)
554             tmp[ordery[i]*N0+j] = X[j*stride+i];
555       }
556    } else {
557       for (i=0;i<stride;i++)
558          for (j=0;j<N0;j++)
559             tmp[i*N0+j] = X[j*stride+i];
560    }
561    for (j=0;j<N;j++)
562       X[j] = tmp[j];
563    RESTORE_STACK;
564 }
565
566 static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
567 {
568    int i,j;
569    VARDECL(celt_norm, tmp);
570    int N;
571    SAVE_STACK;
572    N = N0*stride;
573    ALLOC(tmp, N, celt_norm);
574    if (hadamard)
575    {
576       const int *ordery = ordery_table+stride-2;
577       for (i=0;i<stride;i++)
578          for (j=0;j<N0;j++)
579             tmp[j*stride+i] = X[ordery[i]*N0+j];
580    } else {
581       for (i=0;i<stride;i++)
582          for (j=0;j<N0;j++)
583             tmp[j*stride+i] = X[i*N0+j];
584    }
585    for (j=0;j<N;j++)
586       X[j] = tmp[j];
587    RESTORE_STACK;
588 }
589
590 void haar1(celt_norm *X, int N0, int stride)
591 {
592    int i, j;
593    N0 >>= 1;
594    for (i=0;i<stride;i++)
595       for (j=0;j<N0;j++)
596       {
597          celt_norm tmp1, tmp2;
598          tmp1 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*2*j+i]);
599          tmp2 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*(2*j+1)+i]);
600          X[stride*2*j+i] = tmp1 + tmp2;
601          X[stride*(2*j+1)+i] = tmp1 - tmp2;
602       }
603 }
604
605 static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo)
606 {
607    static const opus_int16 exp2_table8[8] =
608       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
609    int qn, qb;
610    int N2 = 2*N-1;
611    if (stereo && N==2)
612       N2--;
613    /* The upper limit ensures that in a stereo split with itheta==16384, we'll
614        always have enough bits left over to code at least one pulse in the
615        side; otherwise it would collapse, since it doesn't get folded. */
616    qb = IMIN(b-pulse_cap-(4<<BITRES), (b+N2*offset)/N2);
617
618    qb = IMIN(8<<BITRES, qb);
619
620    if (qb<(1<<BITRES>>1)) {
621       qn = 1;
622    } else {
623       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
624       qn = (qn+1)>>1<<1;
625    }
626    celt_assert(qn <= 256);
627    return qn;
628 }
629
630 /* This function is responsible for encoding and decoding a band for both
631    the mono and stereo case. Even in the mono case, it can split the band
632    in two and transmit the energy difference with the two half-bands. It
633    can be called recursively so bands can end up being split in 8 parts. */
634 static unsigned quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
635       int N, int b, int spread, int B, int intensity, int tf_change, celt_norm *lowband, ec_ctx *ec,
636       opus_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level,
637       opus_uint32 *seed, opus_val16 gain, celt_norm *lowband_scratch, int fill)
638 {
639    const unsigned char *cache;
640    int q;
641    int curr_bits;
642    int stereo, split;
643    int imid=0, iside=0;
644    int N0=N;
645    int N_B=N;
646    int N_B0;
647    int B0=B;
648    int time_divide=0;
649    int recombine=0;
650    int inv = 0;
651    opus_val16 mid=0, side=0;
652    int longBlocks;
653    unsigned cm=0;
654 #ifdef RESYNTH
655    int resynth = 1;
656 #else
657    int resynth = !encode;
658 #endif
659
660    longBlocks = B0==1;
661
662    N_B /= B;
663    N_B0 = N_B;
664
665    split = stereo = Y != NULL;
666
667    /* Special case for one sample */
668    if (N==1)
669    {
670       int c;
671       celt_norm *x = X;
672       c=0; do {
673          int sign=0;
674          if (*remaining_bits>=1<<BITRES)
675          {
676             if (encode)
677             {
678                sign = x[0]<0;
679                ec_enc_bits(ec, sign, 1);
680             } else {
681                sign = ec_dec_bits(ec, 1);
682             }
683             *remaining_bits -= 1<<BITRES;
684             b-=1<<BITRES;
685          }
686          if (resynth)
687             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
688          x = Y;
689       } while (++c<1+stereo);
690       if (lowband_out)
691          lowband_out[0] = SHR16(X[0],4);
692       return 1;
693    }
694
695    if (!stereo && level == 0)
696    {
697       int k;
698       if (tf_change>0)
699          recombine = tf_change;
700       /* Band recombining to increase frequency resolution */
701
702       if (lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
703       {
704          int j;
705          for (j=0;j<N;j++)
706             lowband_scratch[j] = lowband[j];
707          lowband = lowband_scratch;
708       }
709
710       for (k=0;k<recombine;k++)
711       {
712          static const unsigned char bit_interleave_table[16]={
713            0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
714          };
715          if (encode)
716             haar1(X, N>>k, 1<<k);
717          if (lowband)
718             haar1(lowband, N>>k, 1<<k);
719          fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
720       }
721       B>>=recombine;
722       N_B<<=recombine;
723
724       /* Increasing the time resolution */
725       while ((N_B&1) == 0 && tf_change<0)
726       {
727          if (encode)
728             haar1(X, N_B, B);
729          if (lowband)
730             haar1(lowband, N_B, B);
731          fill |= fill<<B;
732          B <<= 1;
733          N_B >>= 1;
734          time_divide++;
735          tf_change++;
736       }
737       B0=B;
738       N_B0 = N_B;
739
740       /* Reorganize the samples in time order instead of frequency order */
741       if (B0>1)
742       {
743          if (encode)
744             deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
745          if (lowband)
746             deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
747       }
748    }
749
750    /* If we need 1.5 more bit than we can produce, split the band in two. */
751    cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
752    if (!stereo && LM != -1 && b > cache[cache[0]]+12 && N>2)
753    {
754       N >>= 1;
755       Y = X+N;
756       split = 1;
757       LM -= 1;
758       if (B==1)
759          fill = (fill&1)|(fill<<1);
760       B = (B+1)>>1;
761    }
762
763    if (split)
764    {
765       int qn;
766       int itheta=0;
767       int mbits, sbits, delta;
768       int qalloc;
769       int pulse_cap;
770       int offset;
771       int orig_fill;
772       opus_int32 tell;
773
774       /* Decide on the resolution to give to the split parameter theta */
775       pulse_cap = m->logN[i]+LM*(1<<BITRES);
776       offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
777       qn = compute_qn(N, b, offset, pulse_cap, stereo);
778       if (stereo && i>=intensity)
779          qn = 1;
780       if (encode)
781       {
782          /* theta is the atan() of the ratio between the (normalized)
783             side and mid. With just that parameter, we can re-scale both
784             mid and side because we know that 1) they have unit norm and
785             2) they are orthogonal. */
786          itheta = stereo_itheta(X, Y, stereo, N);
787       }
788       tell = ec_tell_frac(ec);
789       if (qn!=1)
790       {
791          if (encode)
792             itheta = (itheta*qn+8192)>>14;
793
794          /* Entropy coding of the angle. We use a uniform pdf for the
795             time split, a step for stereo, and a triangular one for the rest. */
796          if (stereo && N>2)
797          {
798             int p0 = 3;
799             int x = itheta;
800             int x0 = qn/2;
801             int ft = p0*(x0+1) + x0;
802             /* Use a probability of p0 up to itheta=8192 and then use 1 after */
803             if (encode)
804             {
805                ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
806             } else {
807                int fs;
808                fs=ec_decode(ec,ft);
809                if (fs<(x0+1)*p0)
810                   x=fs/p0;
811                else
812                   x=x0+1+(fs-(x0+1)*p0);
813                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);
814                itheta = x;
815             }
816          } else if (B0>1 || stereo) {
817             /* Uniform pdf */
818             if (encode)
819                ec_enc_uint(ec, itheta, qn+1);
820             else
821                itheta = ec_dec_uint(ec, qn+1);
822          } else {
823             int fs=1, ft;
824             ft = ((qn>>1)+1)*((qn>>1)+1);
825             if (encode)
826             {
827                int fl;
828
829                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
830                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
831                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
832
833                ec_encode(ec, fl, fl+fs, ft);
834             } else {
835                /* Triangular pdf */
836                int fl=0;
837                int fm;
838                fm = ec_decode(ec, ft);
839
840                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
841                {
842                   itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;
843                   fs = itheta + 1;
844                   fl = itheta*(itheta + 1)>>1;
845                }
846                else
847                {
848                   itheta = (2*(qn + 1)
849                    - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;
850                   fs = qn + 1 - itheta;
851                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
852                }
853
854                ec_dec_update(ec, fl, fl+fs, ft);
855             }
856          }
857          itheta = (opus_int32)itheta*16384/qn;
858          if (encode && stereo)
859          {
860             if (itheta==0)
861                intensity_stereo(m, X, Y, bandE, i, N);
862             else
863                stereo_split(X, Y, N);
864          }
865          /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.
866                   Let's do that at higher complexity */
867       } else if (stereo) {
868          if (encode)
869          {
870             inv = itheta > 8192;
871             if (inv)
872             {
873                int j;
874                for (j=0;j<N;j++)
875                   Y[j] = -Y[j];
876             }
877             intensity_stereo(m, X, Y, bandE, i, N);
878          }
879          if (b>2<<BITRES && *remaining_bits > 2<<BITRES)
880          {
881             if (encode)
882                ec_enc_bit_logp(ec, inv, 2);
883             else
884                inv = ec_dec_bit_logp(ec, 2);
885          } else
886             inv = 0;
887          itheta = 0;
888       }
889       qalloc = ec_tell_frac(ec) - tell;
890       b -= qalloc;
891
892       orig_fill = fill;
893       if (itheta == 0)
894       {
895          imid = 32767;
896          iside = 0;
897          fill &= (1<<B)-1;
898          delta = -16384;
899       } else if (itheta == 16384)
900       {
901          imid = 0;
902          iside = 32767;
903          fill &= ((1<<B)-1)<<B;
904          delta = 16384;
905       } else {
906          imid = bitexact_cos(itheta);
907          iside = bitexact_cos(16384-itheta);
908          /* This is the mid vs side allocation that minimizes squared error
909             in that band. */
910          delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid));
911       }
912
913 #ifdef FIXED_POINT
914       mid = imid;
915       side = iside;
916 #else
917       mid = (1.f/32768)*imid;
918       side = (1.f/32768)*iside;
919 #endif
920
921       /* This is a special case for N=2 that only works for stereo and takes
922          advantage of the fact that mid and side are orthogonal to encode
923          the side with just one bit. */
924       if (N==2 && stereo)
925       {
926          int c;
927          int sign=0;
928          celt_norm *x2, *y2;
929          mbits = b;
930          sbits = 0;
931          /* Only need one bit for the side */
932          if (itheta != 0 && itheta != 16384)
933             sbits = 1<<BITRES;
934          mbits -= sbits;
935          c = itheta > 8192;
936          *remaining_bits -= qalloc+sbits;
937
938          x2 = c ? Y : X;
939          y2 = c ? X : Y;
940          if (sbits)
941          {
942             if (encode)
943             {
944                /* Here we only need to encode a sign for the side */
945                sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
946                ec_enc_bits(ec, sign, 1);
947             } else {
948                sign = ec_dec_bits(ec, 1);
949             }
950          }
951          sign = 1-2*sign;
952          /* We use orig_fill here because we want to fold the side, but if
953              itheta==16384, we'll have cleared the low bits of fill. */
954          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);
955          /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
956              and there's no need to worry about mixing with the other channel. */
957          y2[0] = -sign*x2[1];
958          y2[1] = sign*x2[0];
959          if (resynth)
960          {
961             celt_norm tmp;
962             X[0] = MULT16_16_Q15(mid, X[0]);
963             X[1] = MULT16_16_Q15(mid, X[1]);
964             Y[0] = MULT16_16_Q15(side, Y[0]);
965             Y[1] = MULT16_16_Q15(side, Y[1]);
966             tmp = X[0];
967             X[0] = SUB16(tmp,Y[0]);
968             Y[0] = ADD16(tmp,Y[0]);
969             tmp = X[1];
970             X[1] = SUB16(tmp,Y[1]);
971             Y[1] = ADD16(tmp,Y[1]);
972          }
973       } else {
974          /* "Normal" split code */
975          celt_norm *next_lowband2=NULL;
976          celt_norm *next_lowband_out1=NULL;
977          int next_level=0;
978          opus_int32 rebalance;
979
980          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
981          if (B0>1 && !stereo && (itheta&0x3fff))
982          {
983             if (itheta > 8192)
984                /* Rough approximation for pre-echo masking */
985                delta -= delta>>(4-LM);
986             else
987                /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
988                delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
989          }
990          mbits = IMAX(0, IMIN(b, (b-delta)/2));
991          sbits = b-mbits;
992          *remaining_bits -= qalloc;
993
994          if (lowband && !stereo)
995             next_lowband2 = lowband+N; /* >32-bit split case */
996
997          /* Only stereo needs to pass on lowband_out. Otherwise, it's
998             handled at the end */
999          if (stereo)
1000             next_lowband_out1 = lowband_out;
1001          else
1002             next_level = level+1;
1003
1004          rebalance = *remaining_bits;
1005          if (mbits >= sbits)
1006          {
1007             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1008                mid for folding later */
1009             cm = quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1010                   lowband, ec, remaining_bits, LM, next_lowband_out1,
1011                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1012             rebalance = mbits - (rebalance-*remaining_bits);
1013             if (rebalance > 3<<BITRES && itheta!=0)
1014                sbits += rebalance - (3<<BITRES);
1015
1016             /* For a stereo split, the high bits of fill are always zero, so no
1017                folding will be done to the side. */
1018             cm |= quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1019                   next_lowband2, ec, remaining_bits, LM, NULL,
1020                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<((B0>>1)&(stereo-1));
1021          } else {
1022             /* For a stereo split, the high bits of fill are always zero, so no
1023                folding will be done to the side. */
1024             cm = quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1025                   next_lowband2, ec, remaining_bits, LM, NULL,
1026                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<((B0>>1)&(stereo-1));
1027             rebalance = sbits - (rebalance-*remaining_bits);
1028             if (rebalance > 3<<BITRES && itheta!=16384)
1029                mbits += rebalance - (3<<BITRES);
1030             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1031                mid for folding later */
1032             cm |= quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1033                   lowband, ec, remaining_bits, LM, next_lowband_out1,
1034                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1035          }
1036       }
1037
1038    } else {
1039       /* This is the basic no-split case */
1040       q = bits2pulses(m, i, LM, b);
1041       curr_bits = pulses2bits(m, i, LM, q);
1042       *remaining_bits -= curr_bits;
1043
1044       /* Ensures we can never bust the budget */
1045       while (*remaining_bits < 0 && q > 0)
1046       {
1047          *remaining_bits += curr_bits;
1048          q--;
1049          curr_bits = pulses2bits(m, i, LM, q);
1050          *remaining_bits -= curr_bits;
1051       }
1052
1053       if (q!=0)
1054       {
1055          int K = get_pulses(q);
1056
1057          /* Finally do the actual quantization */
1058          if (encode)
1059          {
1060             cm = alg_quant(X, N, K, spread, B, ec
1061 #ifdef RESYNTH
1062                  , gain
1063 #endif
1064                  );
1065          } else {
1066             cm = alg_unquant(X, N, K, spread, B, ec, gain);
1067          }
1068       } else {
1069          /* If there's no pulse, fill the band anyway */
1070          int j;
1071          if (resynth)
1072          {
1073             unsigned cm_mask;
1074             /*B can be as large as 16, so this shift might overflow an int on a
1075                16-bit platform; use a long to get defined behavior.*/
1076             cm_mask = (unsigned)(1UL<<B)-1;
1077             fill &= cm_mask;
1078             if (!fill)
1079             {
1080                for (j=0;j<N;j++)
1081                   X[j] = 0;
1082             } else {
1083                if (lowband == NULL)
1084                {
1085                   /* Noise */
1086                   for (j=0;j<N;j++)
1087                   {
1088                      *seed = celt_lcg_rand(*seed);
1089                      X[j] = (celt_norm)((opus_int32)*seed>>20);
1090                   }
1091                   cm = cm_mask;
1092                } else {
1093                   /* Folded spectrum */
1094                   for (j=0;j<N;j++)
1095                   {
1096                      opus_val16 tmp;
1097                      *seed = celt_lcg_rand(*seed);
1098                      /* About 48 dB below the "normal" folding level */
1099                      tmp = QCONST16(1.0f/256, 10);
1100                      tmp = (*seed)&0x8000 ? tmp : -tmp;
1101                      X[j] = lowband[j]+tmp;
1102                   }
1103                   cm = fill;
1104                }
1105                renormalise_vector(X, N, gain);
1106             }
1107          }
1108       }
1109    }
1110
1111    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1112    if (resynth)
1113    {
1114       if (stereo)
1115       {
1116          if (N!=2)
1117             stereo_merge(X, Y, mid, N);
1118          if (inv)
1119          {
1120             int j;
1121             for (j=0;j<N;j++)
1122                Y[j] = -Y[j];
1123          }
1124       } else if (level == 0)
1125       {
1126          int k;
1127
1128          /* Undo the sample reorganization going from time order to frequency order */
1129          if (B0>1)
1130             interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1131
1132          /* Undo time-freq changes that we did earlier */
1133          N_B = N_B0;
1134          B = B0;
1135          for (k=0;k<time_divide;k++)
1136          {
1137             B >>= 1;
1138             N_B <<= 1;
1139             cm |= cm>>B;
1140             haar1(X, N_B, B);
1141          }
1142
1143          for (k=0;k<recombine;k++)
1144          {
1145             static const unsigned char bit_deinterleave_table[16]={
1146               0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1147               0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1148             };
1149             cm = bit_deinterleave_table[cm];
1150             haar1(X, N0>>k, 1<<k);
1151          }
1152          B<<=recombine;
1153
1154          /* Scale output for later folding */
1155          if (lowband_out)
1156          {
1157             int j;
1158             opus_val16 n;
1159             n = celt_sqrt(SHL32(EXTEND32(N0),22));
1160             for (j=0;j<N0;j++)
1161                lowband_out[j] = MULT16_16_Q15(n,X[j]);
1162          }
1163          cm &= (1<<B)-1;
1164       }
1165    }
1166    return cm;
1167 }
1168
1169 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1170       celt_norm *_X, celt_norm *_Y, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1171       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res,
1172       opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed)
1173 {
1174    int i;
1175    opus_int32 remaining_bits;
1176    const opus_int16 * restrict eBands = m->eBands;
1177    celt_norm * restrict norm, * restrict norm2;
1178    VARDECL(celt_norm, _norm);
1179    VARDECL(celt_norm, lowband_scratch);
1180    int B;
1181    int M;
1182    int lowband_offset;
1183    int update_lowband = 1;
1184    int C = _Y != NULL ? 2 : 1;
1185 #ifdef RESYNTH
1186    int resynth = 1;
1187 #else
1188    int resynth = !encode;
1189 #endif
1190    SAVE_STACK;
1191
1192    M = 1<<LM;
1193    B = shortBlocks ? M : 1;
1194    ALLOC(_norm, C*M*eBands[m->nbEBands], celt_norm);
1195    ALLOC(lowband_scratch, M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]), celt_norm);
1196    norm = _norm;
1197    norm2 = norm + M*eBands[m->nbEBands];
1198
1199    lowband_offset = 0;
1200    for (i=start;i<end;i++)
1201    {
1202       opus_int32 tell;
1203       int b;
1204       int N;
1205       opus_int32 curr_balance;
1206       int effective_lowband=-1;
1207       celt_norm * restrict X, * restrict Y;
1208       int tf_change=0;
1209       unsigned x_cm;
1210       unsigned y_cm;
1211
1212       X = _X+M*eBands[i];
1213       if (_Y!=NULL)
1214          Y = _Y+M*eBands[i];
1215       else
1216          Y = NULL;
1217       N = M*eBands[i+1]-M*eBands[i];
1218       tell = ec_tell_frac(ec);
1219
1220       /* Compute how many bits we want to allocate to this band */
1221       if (i != start)
1222          balance -= tell;
1223       remaining_bits = total_bits-tell-1;
1224       if (i <= codedBands-1)
1225       {
1226          curr_balance = balance / IMIN(3, codedBands-i);
1227          b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1228       } else {
1229          b = 0;
1230       }
1231
1232       if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1233             lowband_offset = i;
1234
1235       tf_change = tf_res[i];
1236       if (i>=m->effEBands)
1237       {
1238          X=norm;
1239          if (_Y!=NULL)
1240             Y = norm;
1241       }
1242
1243       /* Get a conservative estimate of the collapse_mask's for the bands we're
1244           going to be folding from. */
1245       if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1246       {
1247          int fold_start;
1248          int fold_end;
1249          int fold_i;
1250          /* This ensures we never repeat spectral content within one band */
1251          effective_lowband = IMAX(M*eBands[start], M*eBands[lowband_offset]-N);
1252          fold_start = lowband_offset;
1253          while(M*eBands[--fold_start] > effective_lowband);
1254          fold_end = lowband_offset-1;
1255          while(M*eBands[++fold_end] < effective_lowband+N);
1256          x_cm = y_cm = 0;
1257          fold_i = fold_start; do {
1258            x_cm |= collapse_masks[fold_i*C+0];
1259            y_cm |= collapse_masks[fold_i*C+C-1];
1260          } while (++fold_i<fold_end);
1261       }
1262       /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1263           always) be non-zero.*/
1264       else
1265          x_cm = y_cm = (1<<B)-1;
1266
1267       if (dual_stereo && i==intensity)
1268       {
1269          int j;
1270
1271          /* Switch off dual stereo to do intensity */
1272          dual_stereo = 0;
1273          for (j=M*eBands[start];j<M*eBands[i];j++)
1274             norm[j] = HALF32(norm[j]+norm2[j]);
1275       }
1276       if (dual_stereo)
1277       {
1278          x_cm = quant_band(encode, m, i, X, NULL, N, b/2, spread, B, intensity, tf_change,
1279                effective_lowband != -1 ? norm+effective_lowband : NULL, ec, &remaining_bits, LM,
1280                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm);
1281          y_cm = quant_band(encode, m, i, Y, NULL, N, b/2, spread, B, intensity, tf_change,
1282                effective_lowband != -1 ? norm2+effective_lowband : NULL, ec, &remaining_bits, LM,
1283                norm2+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, y_cm);
1284       } else {
1285          x_cm = quant_band(encode, m, i, X, Y, N, b, spread, B, intensity, tf_change,
1286                effective_lowband != -1 ? norm+effective_lowband : NULL, ec, &remaining_bits, LM,
1287                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm|y_cm);
1288          y_cm = x_cm;
1289       }
1290       collapse_masks[i*C+0] = (unsigned char)x_cm;
1291       collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1292       balance += pulses[i] + tell;
1293
1294       /* Update the folding position only as long as we have 1 bit/sample depth */
1295       update_lowband = b>(N<<BITRES);
1296    }
1297    RESTORE_STACK;
1298 }
1299