Renamed some funciton that were likely to clash with other (non-Opus) code
[opus.git] / libcelt / 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)<<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       if (LM>0 || (N&1)==0)
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
766    if (split)
767    {
768       int qn;
769       int itheta=0;
770       int mbits, sbits, delta;
771       int qalloc;
772       int pulse_cap;
773       int offset;
774       int orig_fill;
775       opus_int32 tell;
776
777       /* Decide on the resolution to give to the split parameter theta */
778       pulse_cap = m->logN[i]+(LM<<BITRES);
779       offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
780       qn = compute_qn(N, b, offset, pulse_cap, stereo);
781       if (stereo && i>=intensity)
782          qn = 1;
783       if (encode)
784       {
785          /* theta is the atan() of the ratio between the (normalized)
786             side and mid. With just that parameter, we can re-scale both
787             mid and side because we know that 1) they have unit norm and
788             2) they are orthogonal. */
789          itheta = stereo_itheta(X, Y, stereo, N);
790       }
791       tell = ec_tell_frac(ec);
792       if (qn!=1)
793       {
794          if (encode)
795             itheta = (itheta*qn+8192)>>14;
796
797          /* Entropy coding of the angle. We use a uniform pdf for the
798             time split, a step for stereo, and a triangular one for the rest. */
799          if (stereo && N>2)
800          {
801             int p0 = 3;
802             int x = itheta;
803             int x0 = qn/2;
804             int ft = p0*(x0+1) + x0;
805             /* Use a probability of p0 up to itheta=8192 and then use 1 after */
806             if (encode)
807             {
808                ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
809             } else {
810                int fs;
811                fs=ec_decode(ec,ft);
812                if (fs<(x0+1)*p0)
813                   x=fs/p0;
814                else
815                   x=x0+1+(fs-(x0+1)*p0);
816                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);
817                itheta = x;
818             }
819          } else if (B0>1 || stereo) {
820             /* Uniform pdf */
821             if (encode)
822                ec_enc_uint(ec, itheta, qn+1);
823             else
824                itheta = ec_dec_uint(ec, qn+1);
825          } else {
826             int fs=1, ft;
827             ft = ((qn>>1)+1)*((qn>>1)+1);
828             if (encode)
829             {
830                int fl;
831
832                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
833                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
834                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
835
836                ec_encode(ec, fl, fl+fs, ft);
837             } else {
838                /* Triangular pdf */
839                int fl=0;
840                int fm;
841                fm = ec_decode(ec, ft);
842
843                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
844                {
845                   itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;
846                   fs = itheta + 1;
847                   fl = itheta*(itheta + 1)>>1;
848                }
849                else
850                {
851                   itheta = (2*(qn + 1)
852                    - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;
853                   fs = qn + 1 - itheta;
854                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
855                }
856
857                ec_dec_update(ec, fl, fl+fs, ft);
858             }
859          }
860          itheta = (opus_int32)itheta*16384/qn;
861          if (encode && stereo)
862          {
863             if (itheta==0)
864                intensity_stereo(m, X, Y, bandE, i, N);
865             else
866                stereo_split(X, Y, N);
867          }
868          /* TODO: Renormalising X and Y *may* help fixed-point a bit at very high rate.
869                   Let's do that at higher complexity */
870       } else if (stereo) {
871          if (encode)
872          {
873             inv = itheta > 8192;
874             if (inv)
875             {
876                int j;
877                for (j=0;j<N;j++)
878                   Y[j] = -Y[j];
879             }
880             intensity_stereo(m, X, Y, bandE, i, N);
881          }
882          if (b>2<<BITRES && *remaining_bits > 2<<BITRES)
883          {
884             if (encode)
885                ec_enc_bit_logp(ec, inv, 2);
886             else
887                inv = ec_dec_bit_logp(ec, 2);
888          } else
889             inv = 0;
890          itheta = 0;
891       }
892       qalloc = ec_tell_frac(ec) - tell;
893       b -= qalloc;
894
895       orig_fill = fill;
896       if (itheta == 0)
897       {
898          imid = 32767;
899          iside = 0;
900          fill &= (1<<B)-1;
901          delta = -16384;
902       } else if (itheta == 16384)
903       {
904          imid = 0;
905          iside = 32767;
906          fill &= ((1<<B)-1)<<B;
907          delta = 16384;
908       } else {
909          imid = bitexact_cos(itheta);
910          iside = bitexact_cos(16384-itheta);
911          /* This is the mid vs side allocation that minimizes squared error
912             in that band. */
913          delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid));
914       }
915
916 #ifdef FIXED_POINT
917       mid = imid;
918       side = iside;
919 #else
920       mid = (1.f/32768)*imid;
921       side = (1.f/32768)*iside;
922 #endif
923
924       /* This is a special case for N=2 that only works for stereo and takes
925          advantage of the fact that mid and side are orthogonal to encode
926          the side with just one bit. */
927       if (N==2 && stereo)
928       {
929          int c;
930          int sign=0;
931          celt_norm *x2, *y2;
932          mbits = b;
933          sbits = 0;
934          /* Only need one bit for the side */
935          if (itheta != 0 && itheta != 16384)
936             sbits = 1<<BITRES;
937          mbits -= sbits;
938          c = itheta > 8192;
939          *remaining_bits -= qalloc+sbits;
940
941          x2 = c ? Y : X;
942          y2 = c ? X : Y;
943          if (sbits)
944          {
945             if (encode)
946             {
947                /* Here we only need to encode a sign for the side */
948                sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
949                ec_enc_bits(ec, sign, 1);
950             } else {
951                sign = ec_dec_bits(ec, 1);
952             }
953          }
954          sign = 1-2*sign;
955          /* We use orig_fill here because we want to fold the side, but if
956              itheta==16384, we'll have cleared the low bits of fill. */
957          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);
958          /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
959              and there's no need to worry about mixing with the other channel. */
960          y2[0] = -sign*x2[1];
961          y2[1] = sign*x2[0];
962          if (resynth)
963          {
964             celt_norm tmp;
965             X[0] = MULT16_16_Q15(mid, X[0]);
966             X[1] = MULT16_16_Q15(mid, X[1]);
967             Y[0] = MULT16_16_Q15(side, Y[0]);
968             Y[1] = MULT16_16_Q15(side, Y[1]);
969             tmp = X[0];
970             X[0] = SUB16(tmp,Y[0]);
971             Y[0] = ADD16(tmp,Y[0]);
972             tmp = X[1];
973             X[1] = SUB16(tmp,Y[1]);
974             Y[1] = ADD16(tmp,Y[1]);
975          }
976       } else {
977          /* "Normal" split code */
978          celt_norm *next_lowband2=NULL;
979          celt_norm *next_lowband_out1=NULL;
980          int next_level=0;
981          opus_int32 rebalance;
982
983          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
984          if (B0>1 && !stereo && (itheta&0x3fff))
985          {
986             if (itheta > 8192)
987                /* Rough approximation for pre-echo masking */
988                delta -= delta>>(4-LM);
989             else
990                /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
991                delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
992          }
993          mbits = IMAX(0, IMIN(b, (b-delta)/2));
994          sbits = b-mbits;
995          *remaining_bits -= qalloc;
996
997          if (lowband && !stereo)
998             next_lowband2 = lowband+N; /* >32-bit split case */
999
1000          /* Only stereo needs to pass on lowband_out. Otherwise, it's
1001             handled at the end */
1002          if (stereo)
1003             next_lowband_out1 = lowband_out;
1004          else
1005             next_level = level+1;
1006
1007          rebalance = *remaining_bits;
1008          if (mbits >= sbits)
1009          {
1010             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1011                mid for folding later */
1012             cm = quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1013                   lowband, ec, remaining_bits, LM, next_lowband_out1,
1014                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1015             rebalance = mbits - (rebalance-*remaining_bits);
1016             if (rebalance > 3<<BITRES && itheta!=0)
1017                sbits += rebalance - (3<<BITRES);
1018
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          } else {
1025             /* For a stereo split, the high bits of fill are always zero, so no
1026                folding will be done to the side. */
1027             cm = quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1028                   next_lowband2, ec, remaining_bits, LM, NULL,
1029                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<((B0>>1)&(stereo-1));
1030             rebalance = sbits - (rebalance-*remaining_bits);
1031             if (rebalance > 3<<BITRES && itheta!=16384)
1032                mbits += rebalance - (3<<BITRES);
1033             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1034                mid for folding later */
1035             cm |= quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1036                   lowband, ec, remaining_bits, LM, next_lowband_out1,
1037                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1038          }
1039       }
1040
1041    } else {
1042       /* This is the basic no-split case */
1043       q = bits2pulses(m, i, LM, b);
1044       curr_bits = pulses2bits(m, i, LM, q);
1045       *remaining_bits -= curr_bits;
1046
1047       /* Ensures we can never bust the budget */
1048       while (*remaining_bits < 0 && q > 0)
1049       {
1050          *remaining_bits += curr_bits;
1051          q--;
1052          curr_bits = pulses2bits(m, i, LM, q);
1053          *remaining_bits -= curr_bits;
1054       }
1055
1056       if (q!=0)
1057       {
1058          int K = get_pulses(q);
1059
1060          /* Finally do the actual quantization */
1061          if (encode)
1062             cm = alg_quant(X, N, K, spread, B, ec, gain);
1063          else
1064             cm = alg_unquant(X, N, K, spread, B, ec, gain);
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          for (j=M*eBands[start];j<M*eBands[i];j++)
1271             norm[j] = HALF32(norm[j]+norm2[j]);
1272       }
1273       if (dual_stereo)
1274       {
1275          x_cm = quant_band(encode, m, i, X, NULL, N, b/2, spread, B, intensity, tf_change,
1276                effective_lowband != -1 ? norm+effective_lowband : NULL, ec, &remaining_bits, LM,
1277                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm);
1278          y_cm = quant_band(encode, m, i, Y, NULL, N, b/2, spread, B, intensity, tf_change,
1279                effective_lowband != -1 ? norm2+effective_lowband : NULL, ec, &remaining_bits, LM,
1280                norm2+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, y_cm);
1281       } else {
1282          x_cm = quant_band(encode, m, i, X, Y, N, b, spread, B, intensity, tf_change,
1283                effective_lowband != -1 ? norm+effective_lowband : NULL, ec, &remaining_bits, LM,
1284                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm|y_cm);
1285          y_cm = x_cm;
1286       }
1287       collapse_masks[i*C+0] = (unsigned char)x_cm;
1288       collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1289       balance += pulses[i] + tell;
1290
1291       /* Update the folding position only as long as we have 1 bit/sample depth */
1292       update_lowband = b>(N<<BITRES);
1293    }
1294    RESTORE_STACK;
1295 }
1296