renames the libcelt/ directory to celt/
[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)<<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          /* NOTE: 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          {
1063             cm = alg_quant(X, N, K, spread, B, ec
1064 #ifdef RESYNTH
1065                  , gain
1066 #endif
1067                  );
1068          } else {
1069             cm = alg_unquant(X, N, K, spread, B, ec, gain);
1070          }
1071       } else {
1072          /* If there's no pulse, fill the band anyway */
1073          int j;
1074          if (resynth)
1075          {
1076             unsigned cm_mask;
1077             /*B can be as large as 16, so this shift might overflow an int on a
1078                16-bit platform; use a long to get defined behavior.*/
1079             cm_mask = (unsigned)(1UL<<B)-1;
1080             fill &= cm_mask;
1081             if (!fill)
1082             {
1083                for (j=0;j<N;j++)
1084                   X[j] = 0;
1085             } else {
1086                if (lowband == NULL)
1087                {
1088                   /* Noise */
1089                   for (j=0;j<N;j++)
1090                   {
1091                      *seed = celt_lcg_rand(*seed);
1092                      X[j] = (celt_norm)((opus_int32)*seed>>20);
1093                   }
1094                   cm = cm_mask;
1095                } else {
1096                   /* Folded spectrum */
1097                   for (j=0;j<N;j++)
1098                   {
1099                      opus_val16 tmp;
1100                      *seed = celt_lcg_rand(*seed);
1101                      /* About 48 dB below the "normal" folding level */
1102                      tmp = QCONST16(1.0f/256, 10);
1103                      tmp = (*seed)&0x8000 ? tmp : -tmp;
1104                      X[j] = lowband[j]+tmp;
1105                   }
1106                   cm = fill;
1107                }
1108                renormalise_vector(X, N, gain);
1109             }
1110          }
1111       }
1112    }
1113
1114    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1115    if (resynth)
1116    {
1117       if (stereo)
1118       {
1119          if (N!=2)
1120             stereo_merge(X, Y, mid, N);
1121          if (inv)
1122          {
1123             int j;
1124             for (j=0;j<N;j++)
1125                Y[j] = -Y[j];
1126          }
1127       } else if (level == 0)
1128       {
1129          int k;
1130
1131          /* Undo the sample reorganization going from time order to frequency order */
1132          if (B0>1)
1133             interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1134
1135          /* Undo time-freq changes that we did earlier */
1136          N_B = N_B0;
1137          B = B0;
1138          for (k=0;k<time_divide;k++)
1139          {
1140             B >>= 1;
1141             N_B <<= 1;
1142             cm |= cm>>B;
1143             haar1(X, N_B, B);
1144          }
1145
1146          for (k=0;k<recombine;k++)
1147          {
1148             static const unsigned char bit_deinterleave_table[16]={
1149               0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1150               0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1151             };
1152             cm = bit_deinterleave_table[cm];
1153             haar1(X, N0>>k, 1<<k);
1154          }
1155          B<<=recombine;
1156
1157          /* Scale output for later folding */
1158          if (lowband_out)
1159          {
1160             int j;
1161             opus_val16 n;
1162             n = celt_sqrt(SHL32(EXTEND32(N0),22));
1163             for (j=0;j<N0;j++)
1164                lowband_out[j] = MULT16_16_Q15(n,X[j]);
1165          }
1166          cm &= (1<<B)-1;
1167       }
1168    }
1169    return cm;
1170 }
1171
1172 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1173       celt_norm *_X, celt_norm *_Y, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1174       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res,
1175       opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed)
1176 {
1177    int i;
1178    opus_int32 remaining_bits;
1179    const opus_int16 * restrict eBands = m->eBands;
1180    celt_norm * restrict norm, * restrict norm2;
1181    VARDECL(celt_norm, _norm);
1182    VARDECL(celt_norm, lowband_scratch);
1183    int B;
1184    int M;
1185    int lowband_offset;
1186    int update_lowband = 1;
1187    int C = _Y != NULL ? 2 : 1;
1188 #ifdef RESYNTH
1189    int resynth = 1;
1190 #else
1191    int resynth = !encode;
1192 #endif
1193    SAVE_STACK;
1194
1195    M = 1<<LM;
1196    B = shortBlocks ? M : 1;
1197    ALLOC(_norm, C*M*eBands[m->nbEBands], celt_norm);
1198    ALLOC(lowband_scratch, M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]), celt_norm);
1199    norm = _norm;
1200    norm2 = norm + M*eBands[m->nbEBands];
1201
1202    lowband_offset = 0;
1203    for (i=start;i<end;i++)
1204    {
1205       opus_int32 tell;
1206       int b;
1207       int N;
1208       opus_int32 curr_balance;
1209       int effective_lowband=-1;
1210       celt_norm * restrict X, * restrict Y;
1211       int tf_change=0;
1212       unsigned x_cm;
1213       unsigned y_cm;
1214
1215       X = _X+M*eBands[i];
1216       if (_Y!=NULL)
1217          Y = _Y+M*eBands[i];
1218       else
1219          Y = NULL;
1220       N = M*eBands[i+1]-M*eBands[i];
1221       tell = ec_tell_frac(ec);
1222
1223       /* Compute how many bits we want to allocate to this band */
1224       if (i != start)
1225          balance -= tell;
1226       remaining_bits = total_bits-tell-1;
1227       if (i <= codedBands-1)
1228       {
1229          curr_balance = balance / IMIN(3, codedBands-i);
1230          b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1231       } else {
1232          b = 0;
1233       }
1234
1235       if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1236             lowband_offset = i;
1237
1238       tf_change = tf_res[i];
1239       if (i>=m->effEBands)
1240       {
1241          X=norm;
1242          if (_Y!=NULL)
1243             Y = norm;
1244       }
1245
1246       /* Get a conservative estimate of the collapse_mask's for the bands we're
1247           going to be folding from. */
1248       if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1249       {
1250          int fold_start;
1251          int fold_end;
1252          int fold_i;
1253          /* This ensures we never repeat spectral content within one band */
1254          effective_lowband = IMAX(M*eBands[start], M*eBands[lowband_offset]-N);
1255          fold_start = lowband_offset;
1256          while(M*eBands[--fold_start] > effective_lowband);
1257          fold_end = lowband_offset-1;
1258          while(M*eBands[++fold_end] < effective_lowband+N);
1259          x_cm = y_cm = 0;
1260          fold_i = fold_start; do {
1261            x_cm |= collapse_masks[fold_i*C+0];
1262            y_cm |= collapse_masks[fold_i*C+C-1];
1263          } while (++fold_i<fold_end);
1264       }
1265       /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1266           always) be non-zero.*/
1267       else
1268          x_cm = y_cm = (1<<B)-1;
1269
1270       if (dual_stereo && i==intensity)
1271       {
1272          int j;
1273
1274          /* Switch off dual stereo to do intensity */
1275          dual_stereo = 0;
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