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