Adds a fuzzing mode that causes the encoder to make random decisions
[opus.git] / libcelt / bands.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Copyright (c) 2008-2009 Gregory Maxwell
4    Written by Jean-Marc Valin and Gregory Maxwell */
5 /*
6    Redistribution and use in source and binary forms, with or without
7    modification, are permitted provided that the following conditions
8    are met:
9
10    - Redistributions of source code must retain the above copyright
11    notice, this list of conditions and the following disclaimer.
12
13    - Redistributions in binary form must reproduce the above copyright
14    notice, this list of conditions and the following disclaimer in the
15    documentation and/or other materials provided with the distribution.
16
17    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
18    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
19    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
20    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
21    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
22    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
23    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
24    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
25    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
26    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
27    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30 #ifdef HAVE_CONFIG_H
31 #include "config.h"
32 #endif
33
34 #include <math.h>
35 #include "bands.h"
36 #include "modes.h"
37 #include "vq.h"
38 #include "cwrs.h"
39 #include "stack_alloc.h"
40 #include "os_support.h"
41 #include "mathops.h"
42 #include "rate.h"
43
44 opus_uint32 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    if (tmp > 32767)
57       tmp = 32767;
58    x2 = tmp;
59    x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));
60    if (x2 > 32766)
61       x2 = 32766;
62    return 1+x2;
63 }
64
65 static int bitexact_log2tan(int isin,int icos)
66 {
67    int lc;
68    int ls;
69    lc=EC_ILOG(icos);
70    ls=EC_ILOG(isin);
71    icos<<=15-lc;
72    isin<<=15-ls;
73    return (ls-lc<<11)
74          +FRAC_MUL16(isin, FRAC_MUL16(isin, -2597) + 7932)
75          -FRAC_MUL16(icos, FRAC_MUL16(icos, -2597) + 7932);
76 }
77
78 #ifdef FIXED_POINT
79 /* Compute the amplitude (sqrt energy) in each of the bands */
80 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
81 {
82    int i, c, N;
83    const opus_int16 *eBands = m->eBands;
84    const int C = CHANNELS(_C);
85    N = M*m->shortMdctSize;
86    c=0; do {
87       for (i=0;i<end;i++)
88       {
89          int j;
90          opus_val32 maxval=0;
91          opus_val32 sum = 0;
92
93          j=M*eBands[i]; do {
94             maxval = MAX32(maxval, X[j+c*N]);
95             maxval = MAX32(maxval, -X[j+c*N]);
96          } while (++j<M*eBands[i+1]);
97
98          if (maxval > 0)
99          {
100             int shift = celt_ilog2(maxval)-10;
101             j=M*eBands[i]; do {
102                sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
103                                    EXTRACT16(VSHR32(X[j+c*N],shift)));
104             } while (++j<M*eBands[i+1]);
105             /* We're adding one here to make damn sure we never end up with a pitch vector that's
106                larger than unity norm */
107             bank[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
108          } else {
109             bank[i+c*m->nbEBands] = EPSILON;
110          }
111          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
112       }
113    } while (++c<C);
114    /*printf ("\n");*/
115 }
116
117 /* Normalise each band such that the energy is one. */
118 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)
119 {
120    int i, c, N;
121    const opus_int16 *eBands = m->eBands;
122    const int C = CHANNELS(_C);
123    N = M*m->shortMdctSize;
124    c=0; do {
125       i=0; do {
126          opus_val16 g;
127          int j,shift;
128          opus_val16 E;
129          shift = celt_zlog2(bank[i+c*m->nbEBands])-13;
130          E = VSHR32(bank[i+c*m->nbEBands], shift);
131          g = EXTRACT16(celt_rcp(SHL32(E,3)));
132          j=M*eBands[i]; do {
133             X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
134          } while (++j<M*eBands[i+1]);
135       } while (++i<end);
136    } while (++c<C);
137 }
138
139 #else /* FIXED_POINT */
140 /* Compute the amplitude (sqrt energy) in each of the bands */
141 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
142 {
143    int i, c, N;
144    const opus_int16 *eBands = m->eBands;
145    const int C = CHANNELS(_C);
146    N = M*m->shortMdctSize;
147    c=0; do {
148       for (i=0;i<end;i++)
149       {
150          int j;
151          opus_val32 sum = 1e-27f;
152          for (j=M*eBands[i];j<M*eBands[i+1];j++)
153             sum += X[j+c*N]*X[j+c*N];
154          bank[i+c*m->nbEBands] = celt_sqrt(sum);
155          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
156       }
157    } while (++c<C);
158    /*printf ("\n");*/
159 }
160
161 /* Normalise each band such that the energy is one. */
162 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)
163 {
164    int i, c, N;
165    const opus_int16 *eBands = m->eBands;
166    const int C = CHANNELS(_C);
167    N = M*m->shortMdctSize;
168    c=0; do {
169       for (i=0;i<end;i++)
170       {
171          int j;
172          opus_val16 g = 1.f/(1e-27f+bank[i+c*m->nbEBands]);
173          for (j=M*eBands[i];j<M*eBands[i+1];j++)
174             X[j+c*N] = freq[j+c*N]*g;
175       }
176    } while (++c<C);
177 }
178
179 #endif /* FIXED_POINT */
180
181 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
182 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)
183 {
184    int i, c, N;
185    const opus_int16 *eBands = m->eBands;
186    const int C = CHANNELS(_C);
187    N = M*m->shortMdctSize;
188    celt_assert2(C<=2, "denormalise_bands() not implemented for >2 channels");
189    c=0; do {
190       celt_sig * restrict f;
191       const celt_norm * restrict x;
192       f = freq+c*N;
193       x = X+c*N;
194       for (i=0;i<end;i++)
195       {
196          int j, band_end;
197          opus_val32 g = SHR32(bank[i+c*m->nbEBands],1);
198          j=M*eBands[i];
199          band_end = M*eBands[i+1];
200          do {
201             *f++ = SHL32(MULT16_32_Q15(*x, g),2);
202             x++;
203          } while (++j<band_end);
204       }
205       for (i=M*eBands[end];i<N;i++)
206          *f++ = 0;
207    } while (++c<C);
208 }
209
210 /* This prevents energy collapse for transients with multiple short MDCTs */
211 void anti_collapse(const CELTMode *m, celt_norm *_X, unsigned char *collapse_masks, int LM, int C, int CC, int size,
212       int start, int end, opus_val16 *logE, opus_val16 *prev1logE,
213       opus_val16 *prev2logE, int *pulses, opus_uint32 seed)
214 {
215    int c, i, j, k;
216    for (i=start;i<end;i++)
217    {
218       int N0;
219       opus_val16 thresh, sqrt_1;
220       int depth;
221 #ifdef FIXED_POINT
222       int shift;
223 #endif
224
225       N0 = m->eBands[i+1]-m->eBands[i];
226       /* depth in 1/8 bits */
227       depth = (1+pulses[i])/(m->eBands[i+1]-m->eBands[i]<<LM);
228
229 #ifdef FIXED_POINT
230       thresh = MULT16_32_Q15(QCONST16(0.5f, 15), MIN32(32767,SHR32(celt_exp2(-SHL16(depth, 10-BITRES)),1) ));
231       {
232          opus_val32 t;
233          t = N0<<LM;
234          shift = celt_ilog2(t)>>1;
235          t = SHL32(t, (7-shift)<<1);
236          sqrt_1 = celt_rsqrt_norm(t);
237       }
238 #else
239       thresh = .5f*celt_exp2(-.125f*depth);
240       sqrt_1 = celt_rsqrt(N0<<LM);
241 #endif
242
243       c=0; do
244       {
245          celt_norm *X;
246          opus_val16 prev1;
247          opus_val16 prev2;
248          opus_val16 Ediff;
249          opus_val16 r;
250          int renormalize=0;
251          prev1 = prev1logE[c*m->nbEBands+i];
252          prev2 = prev2logE[c*m->nbEBands+i];
253          if (C<CC)
254          {
255             prev1 = MAX16(prev1,prev1logE[m->nbEBands+i]);
256             prev2 = MAX16(prev2,prev2logE[m->nbEBands+i]);
257          }
258          Ediff = logE[c*m->nbEBands+i]-MIN16(prev1,prev2);
259          Ediff = MAX16(0, Ediff);
260
261 #ifdef FIXED_POINT
262          if (Ediff < 16384)
263             r = 2*MIN16(16383,SHR32(celt_exp2(-Ediff),1));
264          else
265             r = 0;
266          if (LM==3)
267             r = MULT16_16_Q14(23170, MIN32(23169, r));
268          r = SHR16(MIN16(thresh, r),1);
269          r = SHR32(MULT16_16_Q15(sqrt_1, r),shift);
270 #else
271          /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because
272             short blocks don't have the same energy as long */
273          r = 2.f*celt_exp2(-Ediff);
274          if (LM==3)
275             r *= 1.41421356f;
276          r = MIN16(thresh, r);
277          r = r*sqrt_1;
278 #endif
279          X = _X+c*size+(m->eBands[i]<<LM);
280          for (k=0;k<1<<LM;k++)
281          {
282             /* Detect collapse */
283             if (!(collapse_masks[i*C+c]&1<<k))
284             {
285                /* Fill with noise */
286                for (j=0;j<N0;j++)
287                {
288                   seed = lcg_rand(seed);
289                   X[(j<<LM)+k] = (seed&0x8000 ? r : -r);
290                }
291                renormalize = 1;
292             }
293          }
294          /* We just added some energy, so we need to renormalise */
295          if (renormalize)
296             renormalise_vector(X, N0<<LM, Q15ONE);
297       } while (++c<C);
298    }
299 }
300
301 static void intensity_stereo(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bank, int bandID, int N)
302 {
303    int i = bandID;
304    int j;
305    opus_val16 a1, a2;
306    opus_val16 left, right;
307    opus_val16 norm;
308 #ifdef FIXED_POINT
309    int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
310 #endif
311    left = VSHR32(bank[i],shift);
312    right = VSHR32(bank[i+m->nbEBands],shift);
313    norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
314    a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
315    a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);
316    for (j=0;j<N;j++)
317    {
318       celt_norm r, l;
319       l = X[j];
320       r = Y[j];
321       X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
322       /* Side is not encoded, no need to calculate */
323    }
324 }
325
326 static void stereo_split(celt_norm *X, celt_norm *Y, int N)
327 {
328    int j;
329    for (j=0;j<N;j++)
330    {
331       celt_norm r, l;
332       l = MULT16_16_Q15(QCONST16(.70710678f,15), X[j]);
333       r = MULT16_16_Q15(QCONST16(.70710678f,15), Y[j]);
334       X[j] = l+r;
335       Y[j] = r-l;
336    }
337 }
338
339 static void stereo_merge(celt_norm *X, celt_norm *Y, opus_val16 mid, int N)
340 {
341    int j;
342    opus_val32 xp=0, side=0;
343    opus_val32 El, Er;
344    opus_val16 mid2;
345 #ifdef FIXED_POINT
346    int kl, kr;
347 #endif
348    opus_val32 t, lgain, rgain;
349
350    /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
351    for (j=0;j<N;j++)
352    {
353       xp = MAC16_16(xp, X[j], Y[j]);
354       side = MAC16_16(side, Y[j], Y[j]);
355    }
356    /* Compensating for the mid normalization */
357    xp = MULT16_32_Q15(mid, xp);
358    /* mid and side are in Q15, not Q14 like X and Y */
359    mid2 = SHR32(mid, 1);
360    El = MULT16_16(mid2, mid2) + side - 2*xp;
361    Er = MULT16_16(mid2, mid2) + side + 2*xp;
362    if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28))
363    {
364       for (j=0;j<N;j++)
365          Y[j] = X[j];
366       return;
367    }
368
369 #ifdef FIXED_POINT
370    kl = celt_ilog2(El)>>1;
371    kr = celt_ilog2(Er)>>1;
372 #endif
373    t = VSHR32(El, (kl-7)<<1);
374    lgain = celt_rsqrt_norm(t);
375    t = VSHR32(Er, (kr-7)<<1);
376    rgain = celt_rsqrt_norm(t);
377
378 #ifdef FIXED_POINT
379    if (kl < 7)
380       kl = 7;
381    if (kr < 7)
382       kr = 7;
383 #endif
384
385    for (j=0;j<N;j++)
386    {
387       celt_norm r, l;
388       /* Apply mid scaling (side is already scaled) */
389       l = MULT16_16_Q15(mid, X[j]);
390       r = Y[j];
391       X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));
392       Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));
393    }
394 }
395
396 /* Decide whether we should spread the pulses in the current frame */
397 int spreading_decision(const CELTMode *m, celt_norm *X, int *average,
398       int last_decision, int *hf_average, int *tapset_decision, int update_hf,
399       int end, int _C, int M)
400 {
401    int i, c, N0;
402    int sum = 0, nbBands=0;
403    const int C = CHANNELS(_C);
404    const opus_int16 * restrict eBands = m->eBands;
405    int decision;
406    int hf_sum=0;
407
408    N0 = M*m->shortMdctSize;
409
410    if (M*(eBands[end]-eBands[end-1]) <= 8)
411       return SPREAD_NONE;
412    c=0; do {
413       for (i=0;i<end;i++)
414       {
415          int j, N, tmp=0;
416          int tcount[3] = {0,0,0};
417          celt_norm * restrict x = X+M*eBands[i]+c*N0;
418          N = M*(eBands[i+1]-eBands[i]);
419          if (N<=8)
420             continue;
421          /* Compute rough CDF of |x[j]| */
422          for (j=0;j<N;j++)
423          {
424             opus_val32 x2N; /* Q13 */
425
426             x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
427             if (x2N < QCONST16(0.25f,13))
428                tcount[0]++;
429             if (x2N < QCONST16(0.0625f,13))
430                tcount[1]++;
431             if (x2N < QCONST16(0.015625f,13))
432                tcount[2]++;
433          }
434
435          /* Only include four last bands (8 kHz and up) */
436          if (i>m->nbEBands-4)
437             hf_sum += 32*(tcount[1]+tcount[0])/N;
438          tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
439          sum += tmp*256;
440          nbBands++;
441       }
442    } while (++c<C);
443
444    if (update_hf)
445    {
446       if (hf_sum)
447          hf_sum /= C*(4-m->nbEBands+end);
448       *hf_average = (*hf_average+hf_sum)>>1;
449       hf_sum = *hf_average;
450       if (*tapset_decision==2)
451          hf_sum += 4;
452       else if (*tapset_decision==0)
453          hf_sum -= 4;
454       if (hf_sum > 22)
455          *tapset_decision=2;
456       else if (hf_sum > 18)
457          *tapset_decision=1;
458       else
459          *tapset_decision=0;
460    }
461    /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/
462    sum /= nbBands;
463    /* Recursive averaging */
464    sum = (sum+*average)>>1;
465    *average = sum;
466    /* Hysteresis */
467    sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2;
468    if (sum < 80)
469    {
470       decision = SPREAD_AGGRESSIVE;
471    } else if (sum < 256)
472    {
473       decision = SPREAD_NORMAL;
474    } else if (sum < 384)
475    {
476       decision = SPREAD_LIGHT;
477    } else {
478       decision = SPREAD_NONE;
479    }
480 #ifdef FUZZING
481    decision = rand()&0x3;
482    *tapset_decision=rand()%3;
483 #endif
484    return decision;
485 }
486
487 #ifdef MEASURE_NORM_MSE
488
489 float MSE[30] = {0};
490 int nbMSEBands = 0;
491 int MSECount[30] = {0};
492
493 void dump_norm_mse(void)
494 {
495    int i;
496    for (i=0;i<nbMSEBands;i++)
497    {
498       printf ("%g ", MSE[i]/MSECount[i]);
499    }
500    printf ("\n");
501 }
502
503 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
504 {
505    static int init = 0;
506    int i;
507    if (!init)
508    {
509       atexit(dump_norm_mse);
510       init = 1;
511    }
512    for (i=0;i<m->nbEBands;i++)
513    {
514       int j;
515       int c;
516       float g;
517       if (bandE0[i]<10 || (C==2 && bandE0[i+m->nbEBands]<1))
518          continue;
519       c=0; do {
520          g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
521          for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
522             MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
523       } while (++c<C);
524       MSECount[i]+=C;
525    }
526    nbMSEBands = m->nbEBands;
527 }
528
529 #endif
530
531 /* Indexing table for converting from natural Hadamard to ordery Hadamard
532    This is essentially a bit-reversed Gray, on top of which we've added
533    an inversion of the order because we want the DC at the end rather than
534    the beginning. The lines are for N=2, 4, 8, 16 */
535 static const int ordery_table[] = {
536        1,  0,
537        3,  0,  2,  1,
538        7,  0,  4,  3,  6,  1,  5,  2,
539       15,  0,  8,  7, 12,  3, 11,  4, 14,  1,  9,  6, 13,  2, 10,  5,
540 };
541
542 static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
543 {
544    int i,j;
545    VARDECL(celt_norm, tmp);
546    int N;
547    SAVE_STACK;
548    N = N0*stride;
549    ALLOC(tmp, N, celt_norm);
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
568 static 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
592 void 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
607 static 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. */
636 static 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, int resynth, 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
657    longBlocks = B0==1;
658
659    N_B /= B;
660    N_B0 = N_B;
661
662    split = stereo = Y != NULL;
663
664    /* Special case for one sample */
665    if (N==1)
666    {
667       int c;
668       celt_norm *x = X;
669       c=0; do {
670          int sign=0;
671          if (*remaining_bits>=1<<BITRES)
672          {
673             if (encode)
674             {
675                sign = x[0]<0;
676                ec_enc_bits(ec, sign, 1);
677             } else {
678                sign = ec_dec_bits(ec, 1);
679             }
680             *remaining_bits -= 1<<BITRES;
681             b-=1<<BITRES;
682          }
683          if (resynth)
684             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
685          x = Y;
686       } while (++c<1+stereo);
687       if (lowband_out)
688          lowband_out[0] = SHR16(X[0],4);
689       return 1;
690    }
691
692    if (!stereo && level == 0)
693    {
694       int k;
695       if (tf_change>0)
696          recombine = tf_change;
697       /* Band recombining to increase frequency resolution */
698
699       if (lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
700       {
701          int j;
702          for (j=0;j<N;j++)
703             lowband_scratch[j] = lowband[j];
704          lowband = lowband_scratch;
705       }
706
707       for (k=0;k<recombine;k++)
708       {
709          static const unsigned char bit_interleave_table[16]={
710            0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
711          };
712          if (encode)
713             haar1(X, N>>k, 1<<k);
714          if (lowband)
715             haar1(lowband, N>>k, 1<<k);
716          fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
717       }
718       B>>=recombine;
719       N_B<<=recombine;
720
721       /* Increasing the time resolution */
722       while ((N_B&1) == 0 && tf_change<0)
723       {
724          if (encode)
725             haar1(X, N_B, B);
726          if (lowband)
727             haar1(lowband, N_B, B);
728          fill |= fill<<B;
729          B <<= 1;
730          N_B >>= 1;
731          time_divide++;
732          tf_change++;
733       }
734       B0=B;
735       N_B0 = N_B;
736
737       /* Reorganize the samples in time order instead of frequency order */
738       if (B0>1)
739       {
740          if (encode)
741             deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
742          if (lowband)
743             deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
744       }
745    }
746
747    /* If we need 1.5 more bit than we can produce, split the band in two. */
748    cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
749    if (!stereo && LM != -1 && b > cache[cache[0]]+12 && N>2)
750    {
751       if (LM>0 || (N&1)==0)
752       {
753          N >>= 1;
754          Y = X+N;
755          split = 1;
756          LM -= 1;
757          if (B==1)
758             fill = fill&1|fill<<1;
759          B = (B+1)>>1;
760       }
761    }
762
763    if (split)
764    {
765       int qn;
766       int itheta=0;
767       int mbits, sbits, delta;
768       int qalloc;
769       int pulse_cap;
770       int offset;
771       int orig_fill;
772       opus_int32 tell;
773
774       /* Decide on the resolution to give to the split parameter theta */
775       pulse_cap = m->logN[i]+(LM<<BITRES);
776       offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
777       qn = compute_qn(N, b, offset, pulse_cap, stereo);
778       if (stereo && i>=intensity)
779          qn = 1;
780       if (encode)
781       {
782          /* theta is the atan() of the ratio between the (normalized)
783             side and mid. With just that parameter, we can re-scale both
784             mid and side because we know that 1) they have unit norm and
785             2) they are orthogonal. */
786          itheta = stereo_itheta(X, Y, stereo, N);
787       }
788       tell = ec_tell_frac(ec);
789       if (qn!=1)
790       {
791          if (encode)
792             itheta = (itheta*qn+8192)>>14;
793
794          /* Entropy coding of the angle. We use a uniform pdf for the
795             time split, a step for stereo, and a triangular one for the rest. */
796          if (stereo && N>2)
797          {
798             int p0 = 3;
799             int x = itheta;
800             int x0 = qn/2;
801             int ft = p0*(x0+1) + x0;
802             /* Use a probability of p0 up to itheta=8192 and then use 1 after */
803             if (encode)
804             {
805                ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
806             } else {
807                int fs;
808                fs=ec_decode(ec,ft);
809                if (fs<(x0+1)*p0)
810                   x=fs/p0;
811                else
812                   x=x0+1+(fs-(x0+1)*p0);
813                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);
814                itheta = x;
815             }
816          } else if (B0>1 || stereo) {
817             /* Uniform pdf */
818             if (encode)
819                ec_enc_uint(ec, itheta, qn+1);
820             else
821                itheta = ec_dec_uint(ec, qn+1);
822          } else {
823             int fs=1, ft;
824             ft = ((qn>>1)+1)*((qn>>1)+1);
825             if (encode)
826             {
827                int fl;
828
829                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
830                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
831                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
832
833                ec_encode(ec, fl, fl+fs, ft);
834             } else {
835                /* Triangular pdf */
836                int fl=0;
837                int fm;
838                fm = ec_decode(ec, ft);
839
840                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
841                {
842                   itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;
843                   fs = itheta + 1;
844                   fl = itheta*(itheta + 1)>>1;
845                }
846                else
847                {
848                   itheta = (2*(qn + 1)
849                    - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;
850                   fs = qn + 1 - itheta;
851                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
852                }
853
854                ec_dec_update(ec, fl, fl+fs, ft);
855             }
856          }
857          itheta = (opus_int32)itheta*16384/qn;
858          if (encode && stereo)
859          {
860             if (itheta==0)
861                intensity_stereo(m, X, Y, bandE, i, N);
862             else
863                stereo_split(X, Y, N);
864          }
865          /* TODO: Renormalising X and Y *may* help fixed-point a bit at very high rate.
866                   Let's do that at higher complexity */
867       } else if (stereo) {
868          if (encode)
869          {
870             inv = itheta > 8192;
871             if (inv)
872             {
873                int j;
874                for (j=0;j<N;j++)
875                   Y[j] = -Y[j];
876             }
877             intensity_stereo(m, X, Y, bandE, i, N);
878          }
879          if (b>2<<BITRES && *remaining_bits > 2<<BITRES)
880          {
881             if (encode)
882                ec_enc_bit_logp(ec, inv, 2);
883             else
884                inv = ec_dec_bit_logp(ec, 2);
885          } else
886             inv = 0;
887          itheta = 0;
888       }
889       qalloc = ec_tell_frac(ec) - tell;
890       b -= qalloc;
891
892       orig_fill = fill;
893       if (itheta == 0)
894       {
895          imid = 32767;
896          iside = 0;
897          fill &= (1<<B)-1;
898          delta = -16384;
899       } else if (itheta == 16384)
900       {
901          imid = 0;
902          iside = 32767;
903          fill &= (1<<B)-1<<B;
904          delta = 16384;
905       } else {
906          imid = bitexact_cos(itheta);
907          iside = bitexact_cos(16384-itheta);
908          /* This is the mid vs side allocation that minimizes squared error
909             in that band. */
910          delta = FRAC_MUL16(N-1<<7,bitexact_log2tan(iside,imid));
911       }
912
913 #ifdef FIXED_POINT
914       mid = imid;
915       side = iside;
916 #else
917       mid = (1.f/32768)*imid;
918       side = (1.f/32768)*iside;
919 #endif
920
921       /* This is a special case for N=2 that only works for stereo and takes
922          advantage of the fact that mid and side are orthogonal to encode
923          the side with just one bit. */
924       if (N==2 && stereo)
925       {
926          int c;
927          int sign=0;
928          celt_norm *x2, *y2;
929          mbits = b;
930          sbits = 0;
931          /* Only need one bit for the side */
932          if (itheta != 0 && itheta != 16384)
933             sbits = 1<<BITRES;
934          mbits -= sbits;
935          c = itheta > 8192;
936          *remaining_bits -= qalloc+sbits;
937
938          x2 = c ? Y : X;
939          y2 = c ? X : Y;
940          if (sbits)
941          {
942             if (encode)
943             {
944                /* Here we only need to encode a sign for the side */
945                sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
946                ec_enc_bits(ec, sign, 1);
947             } else {
948                sign = ec_dec_bits(ec, 1);
949             }
950          }
951          sign = 1-2*sign;
952          /* We use orig_fill here because we want to fold the side, but if
953              itheta==16384, we'll have cleared the low bits of fill. */
954          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);
955          /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
956              and there's no need to worry about mixing with the other channel. */
957          y2[0] = -sign*x2[1];
958          y2[1] = sign*x2[0];
959          if (resynth)
960          {
961             celt_norm tmp;
962             X[0] = MULT16_16_Q15(mid, X[0]);
963             X[1] = MULT16_16_Q15(mid, X[1]);
964             Y[0] = MULT16_16_Q15(side, Y[0]);
965             Y[1] = MULT16_16_Q15(side, Y[1]);
966             tmp = X[0];
967             X[0] = SUB16(tmp,Y[0]);
968             Y[0] = ADD16(tmp,Y[0]);
969             tmp = X[1];
970             X[1] = SUB16(tmp,Y[1]);
971             Y[1] = ADD16(tmp,Y[1]);
972          }
973       } else {
974          /* "Normal" split code */
975          celt_norm *next_lowband2=NULL;
976          celt_norm *next_lowband_out1=NULL;
977          int next_level=0;
978          opus_int32 rebalance;
979
980          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
981          if (B0>1 && !stereo && (itheta&0x3fff))
982          {
983             if (itheta > 8192)
984                /* Rough approximation for pre-echo masking */
985                delta -= delta>>(4-LM);
986             else
987                /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
988                delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
989          }
990          mbits = IMAX(0, IMIN(b, (b-delta)/2));
991          sbits = b-mbits;
992          *remaining_bits -= qalloc;
993
994          if (lowband && !stereo)
995             next_lowband2 = lowband+N; /* >32-bit split case */
996
997          /* Only stereo needs to pass on lowband_out. Otherwise, it's
998             handled at the end */
999          if (stereo)
1000             next_lowband_out1 = lowband_out;
1001          else
1002             next_level = level+1;
1003
1004          rebalance = *remaining_bits;
1005          if (mbits >= sbits)
1006          {
1007             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1008                mid for folding later */
1009             cm = quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1010                   lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
1011                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1012             rebalance = mbits - (rebalance-*remaining_bits);
1013             if (rebalance > 3<<BITRES && itheta!=0)
1014                sbits += rebalance - (3<<BITRES);
1015
1016             /* For a stereo split, the high bits of fill are always zero, so no
1017                folding will be done to the side. */
1018             cm |= quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1019                   next_lowband2, resynth, ec, remaining_bits, LM, NULL,
1020                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<(B0>>1&stereo-1);
1021          } else {
1022             /* For a stereo split, the high bits of fill are always zero, so no
1023                folding will be done to the side. */
1024             cm = quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1025                   next_lowband2, resynth, ec, remaining_bits, LM, NULL,
1026                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<(B0>>1&stereo-1);
1027             rebalance = sbits - (rebalance-*remaining_bits);
1028             if (rebalance > 3<<BITRES && itheta!=16384)
1029                mbits += rebalance - (3<<BITRES);
1030             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1031                mid for folding later */
1032             cm |= quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1033                   lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
1034                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1035          }
1036       }
1037
1038    } else {
1039       /* This is the basic no-split case */
1040       q = bits2pulses(m, i, LM, b);
1041       curr_bits = pulses2bits(m, i, LM, q);
1042       *remaining_bits -= curr_bits;
1043
1044       /* Ensures we can never bust the budget */
1045       while (*remaining_bits < 0 && q > 0)
1046       {
1047          *remaining_bits += curr_bits;
1048          q--;
1049          curr_bits = pulses2bits(m, i, LM, q);
1050          *remaining_bits -= curr_bits;
1051       }
1052
1053       if (q!=0)
1054       {
1055          int K = get_pulses(q);
1056
1057          /* Finally do the actual quantization */
1058          if (encode)
1059             cm = alg_quant(X, N, K, spread, B, resynth, ec, gain);
1060          else
1061             cm = alg_unquant(X, N, K, spread, B, ec, gain);
1062       } else {
1063          /* If there's no pulse, fill the band anyway */
1064          int j;
1065          if (resynth)
1066          {
1067             unsigned cm_mask;
1068             /*B can be as large as 16, so this shift might overflow an int on a
1069                16-bit platform; use a long to get defined behavior.*/
1070             cm_mask = (unsigned)(1UL<<B)-1;
1071             fill &= cm_mask;
1072             if (!fill)
1073             {
1074                for (j=0;j<N;j++)
1075                   X[j] = 0;
1076             } else {
1077                if (lowband == NULL)
1078                {
1079                   /* Noise */
1080                   for (j=0;j<N;j++)
1081                   {
1082                      *seed = lcg_rand(*seed);
1083                      X[j] = (opus_int32)(*seed)>>20;
1084                   }
1085                   cm = cm_mask;
1086                } else {
1087                   /* Folded spectrum */
1088                   for (j=0;j<N;j++)
1089                   {
1090                      opus_val16 tmp;
1091                      *seed = lcg_rand(*seed);
1092                      /* About 48 dB below the "normal" folding level */
1093                      tmp = QCONST16(1.0f/256, 10);
1094                      tmp = (*seed)&0x8000 ? tmp : -tmp;
1095                      X[j] = lowband[j]+tmp;
1096                   }
1097                   cm = fill;
1098                }
1099                renormalise_vector(X, N, gain);
1100             }
1101          }
1102       }
1103    }
1104
1105    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1106    if (resynth)
1107    {
1108       if (stereo)
1109       {
1110          if (N!=2)
1111             stereo_merge(X, Y, mid, N);
1112          if (inv)
1113          {
1114             int j;
1115             for (j=0;j<N;j++)
1116                Y[j] = -Y[j];
1117          }
1118       } else if (level == 0)
1119       {
1120          int k;
1121
1122          /* Undo the sample reorganization going from time order to frequency order */
1123          if (B0>1)
1124             interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1125
1126          /* Undo time-freq changes that we did earlier */
1127          N_B = N_B0;
1128          B = B0;
1129          for (k=0;k<time_divide;k++)
1130          {
1131             B >>= 1;
1132             N_B <<= 1;
1133             cm |= cm>>B;
1134             haar1(X, N_B, B);
1135          }
1136
1137          for (k=0;k<recombine;k++)
1138          {
1139             static const unsigned char bit_deinterleave_table[16]={
1140               0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1141               0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1142             };
1143             cm = bit_deinterleave_table[cm];
1144             haar1(X, N0>>k, 1<<k);
1145          }
1146          B<<=recombine;
1147
1148          /* Scale output for later folding */
1149          if (lowband_out)
1150          {
1151             int j;
1152             opus_val16 n;
1153             n = celt_sqrt(SHL32(EXTEND32(N0),22));
1154             for (j=0;j<N0;j++)
1155                lowband_out[j] = MULT16_16_Q15(n,X[j]);
1156          }
1157          cm &= (1<<B)-1;
1158       }
1159    }
1160    return cm;
1161 }
1162
1163 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1164       celt_norm *_X, celt_norm *_Y, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1165       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res, int resynth,
1166       opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed)
1167 {
1168    int i;
1169    opus_int32 remaining_bits;
1170    const opus_int16 * restrict eBands = m->eBands;
1171    celt_norm * restrict norm, * restrict norm2;
1172    VARDECL(celt_norm, _norm);
1173    VARDECL(celt_norm, lowband_scratch);
1174    int B;
1175    int M;
1176    int lowband_offset;
1177    int update_lowband = 1;
1178    int C = _Y != NULL ? 2 : 1;
1179    SAVE_STACK;
1180
1181    M = 1<<LM;
1182    B = shortBlocks ? M : 1;
1183    ALLOC(_norm, C*M*eBands[m->nbEBands], celt_norm);
1184    ALLOC(lowband_scratch, M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]), celt_norm);
1185    norm = _norm;
1186    norm2 = norm + M*eBands[m->nbEBands];
1187
1188    lowband_offset = 0;
1189    for (i=start;i<end;i++)
1190    {
1191       opus_int32 tell;
1192       int b;
1193       int N;
1194       opus_int32 curr_balance;
1195       int effective_lowband=-1;
1196       celt_norm * restrict X, * restrict Y;
1197       int tf_change=0;
1198       unsigned x_cm;
1199       unsigned y_cm;
1200
1201       X = _X+M*eBands[i];
1202       if (_Y!=NULL)
1203          Y = _Y+M*eBands[i];
1204       else
1205          Y = NULL;
1206       N = M*eBands[i+1]-M*eBands[i];
1207       tell = ec_tell_frac(ec);
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