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