Renormalize in anti-collapse only when needed
[opus.git] / libcelt / bands.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Copyright (c) 2008-2009 Gregory Maxwell 
4    Written by Jean-Marc Valin and Gregory Maxwell */
5 /*
6    Redistribution and use in source and binary forms, with or without
7    modification, are permitted provided that the following conditions
8    are met:
9    
10    - Redistributions of source code must retain the above copyright
11    notice, this list of conditions and the following disclaimer.
12    
13    - Redistributions in binary form must reproduce the above copyright
14    notice, this list of conditions and the following disclaimer in the
15    documentation and/or other materials provided with the distribution.
16    
17    - Neither the name of the Xiph.org Foundation nor the names of its
18    contributors may be used to endorse or promote products derived from
19    this software without specific prior written permission.
20    
21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
25    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 #ifdef HAVE_CONFIG_H
35 #include "config.h"
36 #endif
37
38 #include <math.h>
39 #include "bands.h"
40 #include "modes.h"
41 #include "vq.h"
42 #include "cwrs.h"
43 #include "stack_alloc.h"
44 #include "os_support.h"
45 #include "mathops.h"
46 #include "rate.h"
47
48 celt_uint32 lcg_rand(celt_uint32 seed)
49 {
50    return 1664525 * seed + 1013904223;
51 }
52
53 /* This is a cos() approximation designed to be bit-exact on any platform. Bit exactness
54    with this approximation is important because it has an impact on the bit allocation */
55 static celt_int16 bitexact_cos(celt_int16 x)
56 {
57    celt_int32 tmp;
58    celt_int16 x2;
59    tmp = (4096+((celt_int32)(x)*(x)))>>13;
60    if (tmp > 32767)
61       tmp = 32767;
62    x2 = tmp;
63    x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));
64    if (x2 > 32766)
65       x2 = 32766;
66    return 1+x2;
67 }
68
69 static int bitexact_log2tan(int isin,int icos)
70 {
71    int lc;
72    int ls;
73    lc=EC_ILOG(icos);
74    ls=EC_ILOG(isin);
75    icos<<=15-lc;
76    isin<<=15-ls;
77    return (ls-lc<<11)
78          +FRAC_MUL16(isin, FRAC_MUL16(isin, -2597) + 7932)
79          -FRAC_MUL16(icos, FRAC_MUL16(icos, -2597) + 7932);
80 }
81
82 #ifdef FIXED_POINT
83 /* Compute the amplitude (sqrt energy) in each of the bands */
84 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
85 {
86    int i, c, N;
87    const celt_int16 *eBands = m->eBands;
88    const int C = CHANNELS(_C);
89    N = M*m->shortMdctSize;
90    c=0; do {
91       for (i=0;i<end;i++)
92       {
93          int j;
94          celt_word32 maxval=0;
95          celt_word32 sum = 0;
96          
97          j=M*eBands[i]; do {
98             maxval = MAX32(maxval, X[j+c*N]);
99             maxval = MAX32(maxval, -X[j+c*N]);
100          } while (++j<M*eBands[i+1]);
101          
102          if (maxval > 0)
103          {
104             int shift = celt_ilog2(maxval)-10;
105             j=M*eBands[i]; do {
106                sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
107                                    EXTRACT16(VSHR32(X[j+c*N],shift)));
108             } while (++j<M*eBands[i+1]);
109             /* We're adding one here to make damn sure we never end up with a pitch vector that's
110                larger than unity norm */
111             bank[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
112          } else {
113             bank[i+c*m->nbEBands] = EPSILON;
114          }
115          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
116       }
117    } while (++c<C);
118    /*printf ("\n");*/
119 }
120
121 /* Normalise each band such that the energy is one. */
122 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)
123 {
124    int i, c, N;
125    const celt_int16 *eBands = m->eBands;
126    const int C = CHANNELS(_C);
127    N = M*m->shortMdctSize;
128    c=0; do {
129       i=0; do {
130          celt_word16 g;
131          int j,shift;
132          celt_word16 E;
133          shift = celt_zlog2(bank[i+c*m->nbEBands])-13;
134          E = VSHR32(bank[i+c*m->nbEBands], shift);
135          g = EXTRACT16(celt_rcp(SHL32(E,3)));
136          j=M*eBands[i]; do {
137             X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
138          } while (++j<M*eBands[i+1]);
139       } while (++i<end);
140    } while (++c<C);
141 }
142
143 #else /* FIXED_POINT */
144 /* Compute the amplitude (sqrt energy) in each of the bands */
145 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
146 {
147    int i, c, N;
148    const celt_int16 *eBands = m->eBands;
149    const int C = CHANNELS(_C);
150    N = M*m->shortMdctSize;
151    c=0; do {
152       for (i=0;i<end;i++)
153       {
154          int j;
155          celt_word32 sum = 1e-27f;
156          for (j=M*eBands[i];j<M*eBands[i+1];j++)
157             sum += X[j+c*N]*X[j+c*N];
158          bank[i+c*m->nbEBands] = celt_sqrt(sum);
159          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
160       }
161    } while (++c<C);
162    /*printf ("\n");*/
163 }
164
165 /* Normalise each band such that the energy is one. */
166 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)
167 {
168    int i, c, N;
169    const celt_int16 *eBands = m->eBands;
170    const int C = CHANNELS(_C);
171    N = M*m->shortMdctSize;
172    c=0; do {
173       for (i=0;i<end;i++)
174       {
175          int j;
176          celt_word16 g = 1.f/(1e-27f+bank[i+c*m->nbEBands]);
177          for (j=M*eBands[i];j<M*eBands[i+1];j++)
178             X[j+c*N] = freq[j+c*N]*g;
179       }
180    } while (++c<C);
181 }
182
183 #endif /* FIXED_POINT */
184
185 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
186 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)
187 {
188    int i, c, N;
189    const celt_int16 *eBands = m->eBands;
190    const int C = CHANNELS(_C);
191    N = M*m->shortMdctSize;
192    celt_assert2(C<=2, "denormalise_bands() not implemented for >2 channels");
193    c=0; do {
194       celt_sig * restrict f;
195       const celt_norm * restrict x;
196       f = freq+c*N;
197       x = X+c*N;
198       for (i=0;i<end;i++)
199       {
200          int j, band_end;
201          celt_word32 g = SHR32(bank[i+c*m->nbEBands],1);
202          j=M*eBands[i];
203          band_end = M*eBands[i+1];
204          do {
205             *f++ = SHL32(MULT16_32_Q15(*x, g),2);
206             x++;
207          } while (++j<band_end);
208       }
209       for (i=M*eBands[m->nbEBands];i<N;i++)
210          *f++ = 0;
211    } while (++c<C);
212 }
213
214 /* This prevents energy collapse for transients with multiple short MDCTs */
215 void anti_collapse(const CELTMode *m, celt_norm *_X, unsigned char *collapse_masks, int LM, int C, int CC, int size,
216       int start, int end, celt_word16 *logE, celt_word16 *prev1logE,
217       celt_word16 *prev2logE, int *pulses, celt_uint32 seed)
218 {
219    int c, i, j, k;
220    for (i=start;i<end;i++)
221    {
222       int N0;
223       celt_word16 thresh, sqrt_1;
224       int depth;
225 #ifdef FIXED_POINT
226       int shift;
227 #endif
228
229       N0 = m->eBands[i+1]-m->eBands[i];
230       /* depth in 1/8 bits */
231       depth = (1+pulses[i])/(m->eBands[i+1]-m->eBands[i]<<LM);
232
233 #ifdef FIXED_POINT
234       thresh = MULT16_32_Q15(QCONST16(0.5f, 15), MIN32(32767,SHR32(celt_exp2(-SHL16(depth, 10-BITRES)),1) ));
235       {
236          celt_word32 t;
237          t = N0<<LM;
238          shift = celt_ilog2(t)>>1;
239          t = SHL32(t, (7-shift)<<1);
240          sqrt_1 = celt_rsqrt_norm(t);
241       }
242 #else
243       thresh = .5f*celt_exp2(-.125f*depth);
244       sqrt_1 = celt_rsqrt(N0<<LM);
245 #endif
246
247       c=0; do
248       {
249          celt_norm *X;
250          celt_word16 prev1;
251          celt_word16 prev2;
252          celt_word16 Ediff;
253          celt_word16 r;
254          int renormalize=0;
255          prev1 = prev1logE[c*m->nbEBands+i];
256          prev2 = prev2logE[c*m->nbEBands+i];
257          if (C<CC)
258          {
259             prev1 = MAX16(prev1,prev1logE[m->nbEBands+i]);
260             prev2 = MAX16(prev2,prev2logE[m->nbEBands+i]);
261          }
262          Ediff = logE[c*m->nbEBands+i]-MIN16(prev1,prev2);
263          Ediff = MAX16(0, Ediff);
264
265 #ifdef FIXED_POINT
266          if (Ediff < 16384)
267             r = 2*MIN16(16383,SHR32(celt_exp2(-Ediff),1));
268          else
269             r = 0;
270          if (LM==3)
271             r = MULT16_16_Q14(23170, MIN32(23169, r));
272          r = SHR16(MIN16(thresh, r),1);
273          r = SHR32(MULT16_16_Q15(sqrt_1, r),shift);
274 #else
275          /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because
276             short blocks don't have the same energy as long */
277          r = 2.f*celt_exp2(-Ediff);
278          if (LM==3)
279             r *= 1.41421356f;
280          r = MIN16(thresh, r);
281          r = r*sqrt_1;
282 #endif
283          X = _X+c*size+(m->eBands[i]<<LM);
284          for (k=0;k<1<<LM;k++)
285          {
286             /* Detect collapse */
287             if (!(collapse_masks[i*C+c]&1<<k))
288             {
289                /* Fill with noise */
290                for (j=0;j<N0;j++)
291                {
292                   seed = lcg_rand(seed);
293                   X[(j<<LM)+k] = (seed&0x8000 ? r : -r);
294                }
295                renormalize = 1;
296             }
297          }
298          /* We just added some energy, so we need to renormalise */
299          if (renormalize)
300             renormalise_vector(X, N0<<LM, Q15ONE);
301       } while (++c<C);
302    }
303
304 }
305
306
307 static void intensity_stereo(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bank, int bandID, int N)
308 {
309    int i = bandID;
310    int j;
311    celt_word16 a1, a2;
312    celt_word16 left, right;
313    celt_word16 norm;
314 #ifdef FIXED_POINT
315    int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
316 #endif
317    left = VSHR32(bank[i],shift);
318    right = VSHR32(bank[i+m->nbEBands],shift);
319    norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
320    a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
321    a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);
322    for (j=0;j<N;j++)
323    {
324       celt_norm r, l;
325       l = X[j];
326       r = Y[j];
327       X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
328       /* Side is not encoded, no need to calculate */
329    }
330 }
331
332 static void stereo_split(celt_norm *X, celt_norm *Y, int N)
333 {
334    int j;
335    for (j=0;j<N;j++)
336    {
337       celt_norm r, l;
338       l = MULT16_16_Q15(QCONST16(.70710678f,15), X[j]);
339       r = MULT16_16_Q15(QCONST16(.70710678f,15), Y[j]);
340       X[j] = l+r;
341       Y[j] = r-l;
342    }
343 }
344
345 static void stereo_merge(celt_norm *X, celt_norm *Y, celt_word16 mid, int N)
346 {
347    int j;
348    celt_word32 xp=0, side=0;
349    celt_word32 El, Er;
350    celt_word16 mid2;
351 #ifdef FIXED_POINT
352    int kl, kr;
353 #endif
354    celt_word32 t, lgain, rgain;
355
356    /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
357    for (j=0;j<N;j++)
358    {
359       xp = MAC16_16(xp, X[j], Y[j]);
360       side = MAC16_16(side, Y[j], Y[j]);
361    }
362    /* Compensating for the mid normalization */
363    xp = MULT16_32_Q15(mid, xp);
364    /* mid and side are in Q15, not Q14 like X and Y */
365    mid2 = SHR32(mid, 1);
366    El = MULT16_16(mid2, mid2) + side - 2*xp;
367    Er = MULT16_16(mid2, mid2) + side + 2*xp;
368    if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28))
369    {
370       for (j=0;j<N;j++)
371          Y[j] = X[j];
372       return;
373    }
374
375 #ifdef FIXED_POINT
376    kl = celt_ilog2(El)>>1;
377    kr = celt_ilog2(Er)>>1;
378 #endif
379    t = VSHR32(El, (kl-7)<<1);
380    lgain = celt_rsqrt_norm(t);
381    t = VSHR32(Er, (kr-7)<<1);
382    rgain = celt_rsqrt_norm(t);
383
384 #ifdef FIXED_POINT
385    if (kl < 7)
386       kl = 7;
387    if (kr < 7)
388       kr = 7;
389 #endif
390
391    for (j=0;j<N;j++)
392    {
393       celt_norm r, l;
394       /* Apply mid scaling (side is already scaled) */
395       l = MULT16_16_Q15(mid, X[j]);
396       r = Y[j];
397       X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));
398       Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));
399    }
400 }
401
402 /* Decide whether we should spread the pulses in the current frame */
403 int spreading_decision(const CELTMode *m, celt_norm *X, int *average,
404       int last_decision, int *hf_average, int *tapset_decision, int update_hf,
405       int end, int _C, int M)
406 {
407    int i, c, N0;
408    int sum = 0, nbBands=0;
409    const int C = CHANNELS(_C);
410    const celt_int16 * restrict eBands = m->eBands;
411    int decision;
412    int hf_sum=0;
413    
414    N0 = M*m->shortMdctSize;
415
416    if (M*(eBands[end]-eBands[end-1]) <= 8)
417       return SPREAD_NONE;
418    c=0; do {
419       for (i=0;i<end;i++)
420       {
421          int j, N, tmp=0;
422          int tcount[3] = {0};
423          celt_norm * restrict x = X+M*eBands[i]+c*N0;
424          N = M*(eBands[i+1]-eBands[i]);
425          if (N<=8)
426             continue;
427          /* Compute rough CDF of |x[j]| */
428          for (j=0;j<N;j++)
429          {
430             celt_word32 x2N; /* Q13 */
431
432             x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
433             if (x2N < QCONST16(0.25f,13))
434                tcount[0]++;
435             if (x2N < QCONST16(0.0625f,13))
436                tcount[1]++;
437             if (x2N < QCONST16(0.015625f,13))
438                tcount[2]++;
439          }
440
441          /* Only include four last bands (8 kHz and up) */
442          if (i>m->nbEBands-4)
443             hf_sum += 32*(tcount[1]+tcount[0])/N;
444          tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
445          sum += tmp*256;
446          nbBands++;
447       }
448    } while (++c<C);
449
450    if (update_hf)
451    {
452       if (hf_sum)
453          hf_sum /= C*(4-m->nbEBands+end);
454       *hf_average = (*hf_average+hf_sum)>>1;
455       hf_sum = *hf_average;
456       if (*tapset_decision==2)
457          hf_sum += 4;
458       else if (*tapset_decision==0)
459          hf_sum -= 4;
460       if (hf_sum > 22)
461          *tapset_decision=2;
462       else if (hf_sum > 18)
463          *tapset_decision=1;
464       else
465          *tapset_decision=0;
466    }
467    /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/
468    sum /= nbBands;
469    /* Recursive averaging */
470    sum = (sum+*average)>>1;
471    *average = sum;
472    /* Hysteresis */
473    sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2;
474    if (sum < 80)
475    {
476       decision = SPREAD_AGGRESSIVE;
477    } else if (sum < 256)
478    {
479       decision = SPREAD_NORMAL;
480    } else if (sum < 384)
481    {
482       decision = SPREAD_LIGHT;
483    } else {
484       decision = SPREAD_NONE;
485    }
486    return decision;
487 }
488
489 #ifdef MEASURE_NORM_MSE
490
491 float MSE[30] = {0};
492 int nbMSEBands = 0;
493 int MSECount[30] = {0};
494
495 void dump_norm_mse(void)
496 {
497    int i;
498    for (i=0;i<nbMSEBands;i++)
499    {
500       printf ("%g ", MSE[i]/MSECount[i]);
501    }
502    printf ("\n");
503 }
504
505 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
506 {
507    static int init = 0;
508    int i;
509    if (!init)
510    {
511       atexit(dump_norm_mse);
512       init = 1;
513    }
514    for (i=0;i<m->nbEBands;i++)
515    {
516       int j;
517       int c;
518       float g;
519       if (bandE0[i]<10 || (C==2 && bandE0[i+m->nbEBands]<1))
520          continue;
521       c=0; do {
522          g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
523          for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
524             MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
525       } while (++c<C);
526       MSECount[i]+=C;
527    }
528    nbMSEBands = m->nbEBands;
529 }
530
531 #endif
532
533 /* Indexing table for converting from natural Hadamard to ordery Hadamard
534    This is essentially a bit-reversed Gray, on top of which we've added
535    an inversion of the order because we want the DC at the end rather than
536    the beginning. The lines are for N=2, 4, 8, 16 */
537 static const int ordery_table[] = {
538        1,  0,
539        3,  0,  2,  1,
540        7,  0,  4,  3,  6,  1,  5,  2,
541       15,  0,  8,  7, 12,  3, 11,  4, 14,  1,  9,  6, 13,  2, 10,  5,
542 };
543
544 static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
545 {
546    int i,j;
547    VARDECL(celt_norm, tmp);
548    int N;
549    SAVE_STACK;
550    N = N0*stride;
551    ALLOC(tmp, N, celt_norm);
552    if (hadamard)
553    {
554       const int *ordery = ordery_table+stride-2;
555       for (i=0;i<stride;i++)
556       {
557          for (j=0;j<N0;j++)
558             tmp[ordery[i]*N0+j] = X[j*stride+i];
559       }
560    } else {
561       for (i=0;i<stride;i++)
562          for (j=0;j<N0;j++)
563             tmp[i*N0+j] = X[j*stride+i];
564    }
565    for (j=0;j<N;j++)
566       X[j] = tmp[j];
567    RESTORE_STACK;
568 }
569
570 static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
571 {
572    int i,j;
573    VARDECL(celt_norm, tmp);
574    int N;
575    SAVE_STACK;
576    N = N0*stride;
577    ALLOC(tmp, N, celt_norm);
578    if (hadamard)
579    {
580       const int *ordery = ordery_table+stride-2;
581       for (i=0;i<stride;i++)
582          for (j=0;j<N0;j++)
583             tmp[j*stride+i] = X[ordery[i]*N0+j];
584    } else {
585       for (i=0;i<stride;i++)
586          for (j=0;j<N0;j++)
587             tmp[j*stride+i] = X[i*N0+j];
588    }
589    for (j=0;j<N;j++)
590       X[j] = tmp[j];
591    RESTORE_STACK;
592 }
593
594 void haar1(celt_norm *X, int N0, int stride)
595 {
596    int i, j;
597    N0 >>= 1;
598    for (i=0;i<stride;i++)
599       for (j=0;j<N0;j++)
600       {
601          celt_norm tmp1, tmp2;
602          tmp1 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*2*j+i]);
603          tmp2 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*(2*j+1)+i]);
604          X[stride*2*j+i] = tmp1 + tmp2;
605          X[stride*(2*j+1)+i] = tmp1 - tmp2;
606       }
607 }
608
609 static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo)
610 {
611    static const celt_int16 exp2_table8[8] =
612       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
613    int qn, qb;
614    int N2 = 2*N-1;
615    if (stereo && N==2)
616       N2--;
617    /* The upper limit ensures that in a stereo split with itheta==16384, we'll
618        always have enough bits left over to code at least one pulse in the
619        side; otherwise it would collapse, since it doesn't get folded. */
620    qb = IMIN(b-pulse_cap-(4<<BITRES), (b+N2*offset)/N2);
621
622    qb = IMIN(8<<BITRES, qb);
623
624    if (qb<(1<<BITRES>>1)) {
625       qn = 1;
626    } else {
627       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
628       qn = (qn+1)>>1<<1;
629    }
630    celt_assert(qn <= 256);
631    return qn;
632 }
633
634 /* This function is responsible for encoding and decoding a band for both
635    the mono and stereo case. Even in the mono case, it can split the band
636    in two and transmit the energy difference with the two half-bands. It
637    can be called recursively so bands can end up being split in 8 parts. */
638 static unsigned quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
639       int N, int b, int spread, int B, int intensity, int tf_change, celt_norm *lowband, int resynth, void *ec,
640       celt_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level,
641       celt_uint32 *seed, celt_word16 gain, celt_norm *lowband_scratch, int fill)
642 {
643    const unsigned char *cache;
644    int q;
645    int curr_bits;
646    int stereo, split;
647    int imid=0, iside=0;
648    int N0=N;
649    int N_B=N;
650    int N_B0;
651    int B0=B;
652    int time_divide=0;
653    int recombine=0;
654    int inv = 0;
655    celt_word16 mid=0, side=0;
656    int longBlocks;
657    unsigned cm=0;
658
659    longBlocks = B0==1;
660
661    N_B /= B;
662    N_B0 = N_B;
663
664    split = stereo = Y != NULL;
665
666    /* Special case for one sample */
667    if (N==1)
668    {
669       int c;
670       celt_norm *x = X;
671       c=0; do {
672          int sign=0;
673          if (*remaining_bits>=1<<BITRES)
674          {
675             if (encode)
676             {
677                sign = x[0]<0;
678                ec_enc_bits((ec_enc*)ec, sign, 1);
679             } else {
680                sign = ec_dec_bits((ec_dec*)ec, 1);
681             }
682             *remaining_bits -= 1<<BITRES;
683             b-=1<<BITRES;
684          }
685          if (resynth)
686             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
687          x = Y;
688       } while (++c<1+stereo);
689       if (lowband_out)
690          lowband_out[0] = SHR16(X[0],4);
691       return 1;
692    }
693
694    if (!stereo && level == 0)
695    {
696       int k;
697       if (tf_change>0)
698          recombine = tf_change;
699       /* Band recombining to increase frequency resolution */
700
701       if (lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
702       {
703          int j;
704          for (j=0;j<N;j++)
705             lowband_scratch[j] = lowband[j];
706          lowband = lowband_scratch;
707       }
708
709       for (k=0;k<recombine;k++)
710       {
711          static const unsigned char bit_interleave_table[16]={
712            0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
713          };
714          if (encode)
715             haar1(X, N>>k, 1<<k);
716          if (lowband)
717             haar1(lowband, N>>k, 1<<k);
718          fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
719       }
720       B>>=recombine;
721       N_B<<=recombine;
722
723       /* Increasing the time resolution */
724       while ((N_B&1) == 0 && tf_change<0)
725       {
726          if (encode)
727             haar1(X, N_B, B);
728          if (lowband)
729             haar1(lowband, N_B, B);
730          fill |= fill<<B;
731          B <<= 1;
732          N_B >>= 1;
733          time_divide++;
734          tf_change++;
735       }
736       B0=B;
737       N_B0 = N_B;
738
739       /* Reorganize the samples in time order instead of frequency order */
740       if (B0>1)
741       {
742          if (encode)
743             deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
744          if (lowband)
745             deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
746       }
747    }
748
749    /* If we need 1.5 more bit than we can produce, split the band in two. */
750    cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
751    if (!stereo && LM != -1 && b > cache[cache[0]]+12 && N>2)
752    {
753       if (LM>0 || (N&1)==0)
754       {
755          N >>= 1;
756          Y = X+N;
757          split = 1;
758          LM -= 1;
759          if (B==1)
760             fill = fill&1|fill<<1;
761          B = (B+1)>>1;
762       }
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       celt_int32 tell;
775
776       /* Decide on the resolution to give to the split parameter theta */
777       pulse_cap = m->logN[i]+(LM<<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 = encode ? ec_enc_tell(ec, BITRES) : ec_dec_tell(ec, BITRES);
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_enc*)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_enc*)ec, itheta, qn+1);
822             else
823                itheta = ec_dec_uint((ec_dec*)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_enc*)ec, fl, fl+fs, ft);
836             } else {
837                /* Triangular pdf */
838                int fl=0;
839                int fm;
840                fm = ec_decode((ec_dec*)ec, ft);
841
842                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
843                {
844                   itheta = (isqrt32(8*(celt_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*(celt_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_dec*)ec, fl, fl+fs, ft);
857             }
858          }
859          itheta = (celt_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          /* TODO: 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 = (encode ? ec_enc_tell(ec, BITRES) : ec_dec_tell(ec, BITRES))
892                - 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_enc*)ec, sign, 1);
950             } else {
951                sign = ec_dec_bits((ec_dec*)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, resynth, 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          celt_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, resynth, 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, resynth, 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, resynth, 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, resynth, ec, remaining_bits, LM, next_lowband_out1,
1037                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1038          }
1039       }
1040
1041    } else {
1042       /* This is the basic no-split case */
1043       q = bits2pulses(m, i, LM, b);
1044       curr_bits = pulses2bits(m, i, LM, q);
1045       *remaining_bits -= curr_bits;
1046
1047       /* Ensures we can never bust the budget */
1048       while (*remaining_bits < 0 && q > 0)
1049       {
1050          *remaining_bits += curr_bits;
1051          q--;
1052          curr_bits = pulses2bits(m, i, LM, q);
1053          *remaining_bits -= curr_bits;
1054       }
1055
1056       if (q!=0)
1057       {
1058          int K = get_pulses(q);
1059
1060          /* Finally do the actual quantization */
1061          if (encode)
1062             cm = alg_quant(X, N, K, spread, B, resynth, (ec_enc*)ec, gain);
1063          else
1064             cm = alg_unquant(X, N, K, spread, B, (ec_dec*)ec, gain);
1065       } else {
1066          /* If there's no pulse, fill the band anyway */
1067          int j;
1068          if (resynth)
1069          {
1070             unsigned cm_mask;
1071             /*B can be as large as 16, so this shift might overflow an int on a
1072                16-bit platform; use a long to get defined behavior.*/
1073             cm_mask = (unsigned)(1UL<<B)-1;
1074             fill &= cm_mask;
1075             if (!fill)
1076             {
1077                for (j=0;j<N;j++)
1078                   X[j] = 0;
1079             } else {
1080                if (lowband == NULL)
1081                {
1082                   /* Noise */
1083                   for (j=0;j<N;j++)
1084                   {
1085                      *seed = lcg_rand(*seed);
1086                      X[j] = (celt_int32)(*seed)>>20;
1087                   }
1088                   cm = cm_mask;
1089                } else {
1090                   /* Folded spectrum */
1091                   for (j=0;j<N;j++)
1092                      X[j] = lowband[j];
1093                   cm = fill;
1094                }
1095                renormalise_vector(X, N, gain);
1096             }
1097          }
1098       }
1099    }
1100
1101    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1102    if (resynth)
1103    {
1104       if (stereo)
1105       {
1106          if (N!=2)
1107             stereo_merge(X, Y, mid, N);
1108          if (inv)
1109          {
1110             int j;
1111             for (j=0;j<N;j++)
1112                Y[j] = -Y[j];
1113          }
1114       } else if (level == 0)
1115       {
1116          int k;
1117
1118          /* Undo the sample reorganization going from time order to frequency order */
1119          if (B0>1)
1120             interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1121
1122          /* Undo time-freq changes that we did earlier */
1123          N_B = N_B0;
1124          B = B0;
1125          for (k=0;k<time_divide;k++)
1126          {
1127             B >>= 1;
1128             N_B <<= 1;
1129             cm |= cm>>B;
1130             haar1(X, N_B, B);
1131          }
1132
1133          for (k=0;k<recombine;k++)
1134          {
1135             static const unsigned char bit_deinterleave_table[16]={
1136               0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1137               0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1138             };
1139             cm = bit_deinterleave_table[cm];
1140             haar1(X, N0>>k, 1<<k);
1141          }
1142          B<<=recombine;
1143          N_B>>=recombine;
1144
1145          /* Scale output for later folding */
1146          if (lowband_out)
1147          {
1148             int j;
1149             celt_word16 n;
1150             n = celt_sqrt(SHL32(EXTEND32(N0),22));
1151             for (j=0;j<N0;j++)
1152                lowband_out[j] = MULT16_16_Q15(n,X[j]);
1153          }
1154          cm &= (1<<B)-1;
1155       }
1156    }
1157    return cm;
1158 }
1159
1160 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1161       celt_norm *_X, celt_norm *_Y, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1162       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res, int resynth,
1163       celt_int32 total_bits, celt_int32 balance, void *ec, int LM, int codedBands, ec_uint32 *seed)
1164 {
1165    int i;
1166    celt_int32 remaining_bits;
1167    const celt_int16 * restrict eBands = m->eBands;
1168    celt_norm * restrict norm, * restrict norm2;
1169    VARDECL(celt_norm, _norm);
1170    VARDECL(celt_norm, lowband_scratch);
1171    int B;
1172    int M;
1173    int lowband_offset;
1174    int update_lowband = 1;
1175    int C = _Y != NULL ? 2 : 1;
1176    SAVE_STACK;
1177
1178    M = 1<<LM;
1179    B = shortBlocks ? M : 1;
1180    ALLOC(_norm, C*M*eBands[m->nbEBands], celt_norm);
1181    ALLOC(lowband_scratch, M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]), celt_norm);
1182    norm = _norm;
1183    norm2 = norm + M*eBands[m->nbEBands];
1184
1185    lowband_offset = 0;
1186    for (i=start;i<end;i++)
1187    {
1188       celt_int32 tell;
1189       int b;
1190       int N;
1191       celt_int32 curr_balance;
1192       int effective_lowband=-1;
1193       celt_norm * restrict X, * restrict Y;
1194       int tf_change=0;
1195       unsigned x_cm;
1196       unsigned y_cm;
1197
1198       X = _X+M*eBands[i];
1199       if (_Y!=NULL)
1200          Y = _Y+M*eBands[i];
1201       else
1202          Y = NULL;
1203       N = M*eBands[i+1]-M*eBands[i];
1204       if (encode)
1205          tell = ec_enc_tell((ec_enc*)ec, BITRES);
1206       else
1207          tell = ec_dec_tell((ec_dec*)ec, BITRES);
1208
1209       /* Compute how many bits we want to allocate to this band */
1210       if (i != start)
1211          balance -= tell;
1212       remaining_bits = total_bits-tell-1;
1213       if (i <= codedBands-1)
1214       {
1215          curr_balance = balance / IMIN(3, codedBands-i);
1216          b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1217       } else {
1218          b = 0;
1219       }
1220
1221       if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1222             lowband_offset = i;
1223
1224       tf_change = tf_res[i];
1225       if (i>=m->effEBands)
1226       {
1227          X=norm;
1228          if (_Y!=NULL)
1229             Y = norm;
1230       }
1231
1232       /* Get a conservative estimate of the collapse_mask's for the bands we're
1233           going to be folding from. */
1234       if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1235       {
1236          int fold_start;
1237          int fold_end;
1238          int fold_i;
1239          /* This ensures we never repeat spectral content within one band */
1240          effective_lowband = IMAX(M*eBands[start], M*eBands[lowband_offset]-N);
1241          fold_start = lowband_offset;
1242          while(M*eBands[--fold_start] > effective_lowband);
1243          fold_end = lowband_offset-1;
1244          while(M*eBands[++fold_end] < effective_lowband+N);
1245          x_cm = y_cm = 0;
1246          fold_i = fold_start; do {
1247            x_cm |= collapse_masks[fold_i*C+0];
1248            y_cm |= collapse_masks[fold_i*C+C-1];
1249          } while (++fold_i<fold_end);
1250       }
1251       /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1252           always) be non-zero.*/
1253       else
1254          x_cm = y_cm = (1<<B)-1;
1255
1256       if (dual_stereo && i==intensity)
1257       {
1258          int j;
1259
1260          /* Switch off dual stereo to do intensity */
1261          dual_stereo = 0;
1262          for (j=M*eBands[start];j<M*eBands[i];j++)
1263             norm[j] = HALF32(norm[j]+norm2[j]);
1264       }
1265       if (dual_stereo)
1266       {
1267          x_cm = quant_band(encode, m, i, X, NULL, N, b/2, spread, B, intensity, tf_change,
1268                effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1269                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm);
1270          y_cm = quant_band(encode, m, i, Y, NULL, N, b/2, spread, B, intensity, tf_change,
1271                effective_lowband != -1 ? norm2+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1272                norm2+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, y_cm);
1273       } else {
1274          x_cm = quant_band(encode, m, i, X, Y, N, b, spread, B, intensity, tf_change,
1275                effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1276                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm|y_cm);
1277          y_cm = x_cm;
1278       }
1279       collapse_masks[i*C+0] = (unsigned char)x_cm;
1280       collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1281       balance += pulses[i] + tell;
1282
1283       /* Update the folding position only as long as we have 1 bit/sample depth */
1284       update_lowband = b>(N<<BITRES);
1285    }
1286    RESTORE_STACK;
1287 }
1288