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