Adding a tiny amount of noise when folding to prevent/reduce numerical issues
[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 celt_uint32 lcg_rand(celt_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 celt_int16 bitexact_cos(celt_int16 x)
52 {
53    celt_int32 tmp;
54    celt_int16 x2;
55    tmp = (4096+((celt_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 celt_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          celt_word32 maxval=0;
91          celt_word32 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 celt_int16 *eBands = m->eBands;
122    const int C = CHANNELS(_C);
123    N = M*m->shortMdctSize;
124    c=0; do {
125       i=0; do {
126          celt_word16 g;
127          int j,shift;
128          celt_word16 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 celt_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          celt_word32 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 celt_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          celt_word16 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 celt_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          celt_word32 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[m->nbEBands];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, celt_word16 *logE, celt_word16 *prev1logE,
213       celt_word16 *prev2logE, int *pulses, celt_uint32 seed)
214 {
215    int c, i, j, k;
216    for (i=start;i<end;i++)
217    {
218       int N0;
219       celt_word16 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          celt_word32 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          celt_word16 prev1;
247          celt_word16 prev2;
248          celt_word16 Ediff;
249          celt_word16 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
302
303 static void intensity_stereo(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bank, int bandID, int N)
304 {
305    int i = bandID;
306    int j;
307    celt_word16 a1, a2;
308    celt_word16 left, right;
309    celt_word16 norm;
310 #ifdef FIXED_POINT
311    int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
312 #endif
313    left = VSHR32(bank[i],shift);
314    right = VSHR32(bank[i+m->nbEBands],shift);
315    norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
316    a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
317    a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);
318    for (j=0;j<N;j++)
319    {
320       celt_norm r, l;
321       l = X[j];
322       r = Y[j];
323       X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
324       /* Side is not encoded, no need to calculate */
325    }
326 }
327
328 static void stereo_split(celt_norm *X, celt_norm *Y, int N)
329 {
330    int j;
331    for (j=0;j<N;j++)
332    {
333       celt_norm r, l;
334       l = MULT16_16_Q15(QCONST16(.70710678f,15), X[j]);
335       r = MULT16_16_Q15(QCONST16(.70710678f,15), Y[j]);
336       X[j] = l+r;
337       Y[j] = r-l;
338    }
339 }
340
341 static void stereo_merge(celt_norm *X, celt_norm *Y, celt_word16 mid, int N)
342 {
343    int j;
344    celt_word32 xp=0, side=0;
345    celt_word32 El, Er;
346    celt_word16 mid2;
347 #ifdef FIXED_POINT
348    int kl, kr;
349 #endif
350    celt_word32 t, lgain, rgain;
351
352    /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
353    for (j=0;j<N;j++)
354    {
355       xp = MAC16_16(xp, X[j], Y[j]);
356       side = MAC16_16(side, Y[j], Y[j]);
357    }
358    /* Compensating for the mid normalization */
359    xp = MULT16_32_Q15(mid, xp);
360    /* mid and side are in Q15, not Q14 like X and Y */
361    mid2 = SHR32(mid, 1);
362    El = MULT16_16(mid2, mid2) + side - 2*xp;
363    Er = MULT16_16(mid2, mid2) + side + 2*xp;
364    if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28))
365    {
366       for (j=0;j<N;j++)
367          Y[j] = X[j];
368       return;
369    }
370
371 #ifdef FIXED_POINT
372    kl = celt_ilog2(El)>>1;
373    kr = celt_ilog2(Er)>>1;
374 #endif
375    t = VSHR32(El, (kl-7)<<1);
376    lgain = celt_rsqrt_norm(t);
377    t = VSHR32(Er, (kr-7)<<1);
378    rgain = celt_rsqrt_norm(t);
379
380 #ifdef FIXED_POINT
381    if (kl < 7)
382       kl = 7;
383    if (kr < 7)
384       kr = 7;
385 #endif
386
387    for (j=0;j<N;j++)
388    {
389       celt_norm r, l;
390       /* Apply mid scaling (side is already scaled) */
391       l = MULT16_16_Q15(mid, X[j]);
392       r = Y[j];
393       X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));
394       Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));
395    }
396 }
397
398 /* Decide whether we should spread the pulses in the current frame */
399 int spreading_decision(const CELTMode *m, celt_norm *X, int *average,
400       int last_decision, int *hf_average, int *tapset_decision, int update_hf,
401       int end, int _C, int M)
402 {
403    int i, c, N0;
404    int sum = 0, nbBands=0;
405    const int C = CHANNELS(_C);
406    const celt_int16 * restrict eBands = m->eBands;
407    int decision;
408    int hf_sum=0;
409    
410    N0 = M*m->shortMdctSize;
411
412    if (M*(eBands[end]-eBands[end-1]) <= 8)
413       return SPREAD_NONE;
414    c=0; do {
415       for (i=0;i<end;i++)
416       {
417          int j, N, tmp=0;
418          int tcount[3] = {0};
419          celt_norm * restrict x = X+M*eBands[i]+c*N0;
420          N = M*(eBands[i+1]-eBands[i]);
421          if (N<=8)
422             continue;
423          /* Compute rough CDF of |x[j]| */
424          for (j=0;j<N;j++)
425          {
426             celt_word32 x2N; /* Q13 */
427
428             x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
429             if (x2N < QCONST16(0.25f,13))
430                tcount[0]++;
431             if (x2N < QCONST16(0.0625f,13))
432                tcount[1]++;
433             if (x2N < QCONST16(0.015625f,13))
434                tcount[2]++;
435          }
436
437          /* Only include four last bands (8 kHz and up) */
438          if (i>m->nbEBands-4)
439             hf_sum += 32*(tcount[1]+tcount[0])/N;
440          tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
441          sum += tmp*256;
442          nbBands++;
443       }
444    } while (++c<C);
445
446    if (update_hf)
447    {
448       if (hf_sum)
449          hf_sum /= C*(4-m->nbEBands+end);
450       *hf_average = (*hf_average+hf_sum)>>1;
451       hf_sum = *hf_average;
452       if (*tapset_decision==2)
453          hf_sum += 4;
454       else if (*tapset_decision==0)
455          hf_sum -= 4;
456       if (hf_sum > 22)
457          *tapset_decision=2;
458       else if (hf_sum > 18)
459          *tapset_decision=1;
460       else
461          *tapset_decision=0;
462    }
463    /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/
464    sum /= nbBands;
465    /* Recursive averaging */
466    sum = (sum+*average)>>1;
467    *average = sum;
468    /* Hysteresis */
469    sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2;
470    if (sum < 80)
471    {
472       decision = SPREAD_AGGRESSIVE;
473    } else if (sum < 256)
474    {
475       decision = SPREAD_NORMAL;
476    } else if (sum < 384)
477    {
478       decision = SPREAD_LIGHT;
479    } else {
480       decision = SPREAD_NONE;
481    }
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 celt_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, int resynth, ec_ctx *ec,
636       celt_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level,
637       celt_uint32 *seed, celt_word16 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    celt_word16 mid=0, side=0;
652    int longBlocks;
653    unsigned cm=0;
654
655    longBlocks = B0==1;
656
657    N_B /= B;
658    N_B0 = N_B;
659
660    split = stereo = Y != NULL;
661
662    /* Special case for one sample */
663    if (N==1)
664    {
665       int c;
666       celt_norm *x = X;
667       c=0; do {
668          int sign=0;
669          if (*remaining_bits>=1<<BITRES)
670          {
671             if (encode)
672             {
673                sign = x[0]<0;
674                ec_enc_bits(ec, sign, 1);
675             } else {
676                sign = ec_dec_bits(ec, 1);
677             }
678             *remaining_bits -= 1<<BITRES;
679             b-=1<<BITRES;
680          }
681          if (resynth)
682             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
683          x = Y;
684       } while (++c<1+stereo);
685       if (lowband_out)
686          lowband_out[0] = SHR16(X[0],4);
687       return 1;
688    }
689
690    if (!stereo && level == 0)
691    {
692       int k;
693       if (tf_change>0)
694          recombine = tf_change;
695       /* Band recombining to increase frequency resolution */
696
697       if (lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
698       {
699          int j;
700          for (j=0;j<N;j++)
701             lowband_scratch[j] = lowband[j];
702          lowband = lowband_scratch;
703       }
704
705       for (k=0;k<recombine;k++)
706       {
707          static const unsigned char bit_interleave_table[16]={
708            0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
709          };
710          if (encode)
711             haar1(X, N>>k, 1<<k);
712          if (lowband)
713             haar1(lowband, N>>k, 1<<k);
714          fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
715       }
716       B>>=recombine;
717       N_B<<=recombine;
718
719       /* Increasing the time resolution */
720       while ((N_B&1) == 0 && tf_change<0)
721       {
722          if (encode)
723             haar1(X, N_B, B);
724          if (lowband)
725             haar1(lowband, N_B, B);
726          fill |= fill<<B;
727          B <<= 1;
728          N_B >>= 1;
729          time_divide++;
730          tf_change++;
731       }
732       B0=B;
733       N_B0 = N_B;
734
735       /* Reorganize the samples in time order instead of frequency order */
736       if (B0>1)
737       {
738          if (encode)
739             deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
740          if (lowband)
741             deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
742       }
743    }
744
745    /* If we need 1.5 more bit than we can produce, split the band in two. */
746    cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
747    if (!stereo && LM != -1 && b > cache[cache[0]]+12 && N>2)
748    {
749       if (LM>0 || (N&1)==0)
750       {
751          N >>= 1;
752          Y = X+N;
753          split = 1;
754          LM -= 1;
755          if (B==1)
756             fill = fill&1|fill<<1;
757          B = (B+1)>>1;
758       }
759    }
760
761    if (split)
762    {
763       int qn;
764       int itheta=0;
765       int mbits, sbits, delta;
766       int qalloc;
767       int pulse_cap;
768       int offset;
769       int orig_fill;
770       celt_int32 tell;
771
772       /* Decide on the resolution to give to the split parameter theta */
773       pulse_cap = m->logN[i]+(LM<<BITRES);
774       offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
775       qn = compute_qn(N, b, offset, pulse_cap, stereo);
776       if (stereo && i>=intensity)
777          qn = 1;
778       if (encode)
779       {
780          /* theta is the atan() of the ratio between the (normalized)
781             side and mid. With just that parameter, we can re-scale both
782             mid and side because we know that 1) they have unit norm and
783             2) they are orthogonal. */
784          itheta = stereo_itheta(X, Y, stereo, N);
785       }
786       tell = ec_tell_frac(ec);
787       if (qn!=1)
788       {
789          if (encode)
790             itheta = (itheta*qn+8192)>>14;
791
792          /* Entropy coding of the angle. We use a uniform pdf for the
793             time split, a step for stereo, and a triangular one for the rest. */
794          if (stereo && N>2)
795          {
796             int p0 = 3;
797             int x = itheta;
798             int x0 = qn/2;
799             int ft = p0*(x0+1) + x0;
800             /* Use a probability of p0 up to itheta=8192 and then use 1 after */
801             if (encode)
802             {
803                ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
804             } else {
805                int fs;
806                fs=ec_decode(ec,ft);
807                if (fs<(x0+1)*p0)
808                   x=fs/p0;
809                else
810                   x=x0+1+(fs-(x0+1)*p0);
811                ec_dec_update(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
812                itheta = x;
813             }
814          } else if (B0>1 || stereo) {
815             /* Uniform pdf */
816             if (encode)
817                ec_enc_uint(ec, itheta, qn+1);
818             else
819                itheta = ec_dec_uint(ec, qn+1);
820          } else {
821             int fs=1, ft;
822             ft = ((qn>>1)+1)*((qn>>1)+1);
823             if (encode)
824             {
825                int fl;
826
827                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
828                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
829                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
830
831                ec_encode(ec, fl, fl+fs, ft);
832             } else {
833                /* Triangular pdf */
834                int fl=0;
835                int fm;
836                fm = ec_decode(ec, ft);
837
838                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
839                {
840                   itheta = (isqrt32(8*(celt_uint32)fm + 1) - 1)>>1;
841                   fs = itheta + 1;
842                   fl = itheta*(itheta + 1)>>1;
843                }
844                else
845                {
846                   itheta = (2*(qn + 1)
847                    - isqrt32(8*(celt_uint32)(ft - fm - 1) + 1))>>1;
848                   fs = qn + 1 - itheta;
849                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
850                }
851
852                ec_dec_update(ec, fl, fl+fs, ft);
853             }
854          }
855          itheta = (celt_int32)itheta*16384/qn;
856          if (encode && stereo)
857          {
858             if (itheta==0)
859                intensity_stereo(m, X, Y, bandE, i, N);
860             else
861                stereo_split(X, Y, N);
862          }
863          /* TODO: Renormalising X and Y *may* help fixed-point a bit at very high rate.
864                   Let's do that at higher complexity */
865       } else if (stereo) {
866          if (encode)
867          {
868             inv = itheta > 8192;
869             if (inv)
870             {
871                int j;
872                for (j=0;j<N;j++)
873                   Y[j] = -Y[j];
874             }
875             intensity_stereo(m, X, Y, bandE, i, N);
876          }
877          if (b>2<<BITRES && *remaining_bits > 2<<BITRES)
878          {
879             if (encode)
880                ec_enc_bit_logp(ec, inv, 2);
881             else
882                inv = ec_dec_bit_logp(ec, 2);
883          } else
884             inv = 0;
885          itheta = 0;
886       }
887       qalloc = ec_tell_frac(ec) - tell;
888       b -= qalloc;
889
890       orig_fill = fill;
891       if (itheta == 0)
892       {
893          imid = 32767;
894          iside = 0;
895          fill &= (1<<B)-1;
896          delta = -16384;
897       } else if (itheta == 16384)
898       {
899          imid = 0;
900          iside = 32767;
901          fill &= (1<<B)-1<<B;
902          delta = 16384;
903       } else {
904          imid = bitexact_cos(itheta);
905          iside = bitexact_cos(16384-itheta);
906          /* This is the mid vs side allocation that minimizes squared error
907             in that band. */
908          delta = FRAC_MUL16(N-1<<7,bitexact_log2tan(iside,imid));
909       }
910
911 #ifdef FIXED_POINT
912       mid = imid;
913       side = iside;
914 #else
915       mid = (1.f/32768)*imid;
916       side = (1.f/32768)*iside;
917 #endif
918
919       /* This is a special case for N=2 that only works for stereo and takes
920          advantage of the fact that mid and side are orthogonal to encode
921          the side with just one bit. */
922       if (N==2 && stereo)
923       {
924          int c;
925          int sign=0;
926          celt_norm *x2, *y2;
927          mbits = b;
928          sbits = 0;
929          /* Only need one bit for the side */
930          if (itheta != 0 && itheta != 16384)
931             sbits = 1<<BITRES;
932          mbits -= sbits;
933          c = itheta > 8192;
934          *remaining_bits -= qalloc+sbits;
935
936          x2 = c ? Y : X;
937          y2 = c ? X : Y;
938          if (sbits)
939          {
940             if (encode)
941             {
942                /* Here we only need to encode a sign for the side */
943                sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
944                ec_enc_bits(ec, sign, 1);
945             } else {
946                sign = ec_dec_bits(ec, 1);
947             }
948          }
949          sign = 1-2*sign;
950          /* We use orig_fill here because we want to fold the side, but if
951              itheta==16384, we'll have cleared the low bits of fill. */
952          cm = quant_band(encode, m, i, x2, NULL, N, mbits, spread, B, intensity, tf_change, lowband, resynth, ec, remaining_bits, LM, lowband_out, NULL, level, seed, gain, lowband_scratch, orig_fill);
953          /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
954              and there's no need to worry about mixing with the other channel. */
955          y2[0] = -sign*x2[1];
956          y2[1] = sign*x2[0];
957          if (resynth)
958          {
959             celt_norm tmp;
960             X[0] = MULT16_16_Q15(mid, X[0]);
961             X[1] = MULT16_16_Q15(mid, X[1]);
962             Y[0] = MULT16_16_Q15(side, Y[0]);
963             Y[1] = MULT16_16_Q15(side, Y[1]);
964             tmp = X[0];
965             X[0] = SUB16(tmp,Y[0]);
966             Y[0] = ADD16(tmp,Y[0]);
967             tmp = X[1];
968             X[1] = SUB16(tmp,Y[1]);
969             Y[1] = ADD16(tmp,Y[1]);
970          }
971       } else {
972          /* "Normal" split code */
973          celt_norm *next_lowband2=NULL;
974          celt_norm *next_lowband_out1=NULL;
975          int next_level=0;
976          celt_int32 rebalance;
977
978          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
979          if (B0>1 && !stereo && (itheta&0x3fff))
980          {
981             if (itheta > 8192)
982                /* Rough approximation for pre-echo masking */
983                delta -= delta>>(4-LM);
984             else
985                /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
986                delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
987          }
988          mbits = IMAX(0, IMIN(b, (b-delta)/2));
989          sbits = b-mbits;
990          *remaining_bits -= qalloc;
991
992          if (lowband && !stereo)
993             next_lowband2 = lowband+N; /* >32-bit split case */
994
995          /* Only stereo needs to pass on lowband_out. Otherwise, it's
996             handled at the end */
997          if (stereo)
998             next_lowband_out1 = lowband_out;
999          else
1000             next_level = level+1;
1001
1002          rebalance = *remaining_bits;
1003          if (mbits >= sbits)
1004          {
1005             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1006                mid for folding later */
1007             cm = quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1008                   lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
1009                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1010             rebalance = mbits - (rebalance-*remaining_bits);
1011             if (rebalance > 3<<BITRES && itheta!=0)
1012                sbits += rebalance - (3<<BITRES);
1013
1014             /* For a stereo split, the high bits of fill are always zero, so no
1015                folding will be done to the side. */
1016             cm |= quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1017                   next_lowband2, resynth, ec, remaining_bits, LM, NULL,
1018                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<(B0>>1&stereo-1);
1019          } else {
1020             /* For a stereo split, the high bits of fill are always zero, so no
1021                folding will be done to the side. */
1022             cm = quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1023                   next_lowband2, resynth, ec, remaining_bits, LM, NULL,
1024                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<(B0>>1&stereo-1);
1025             rebalance = sbits - (rebalance-*remaining_bits);
1026             if (rebalance > 3<<BITRES && itheta!=16384)
1027                mbits += rebalance - (3<<BITRES);
1028             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1029                mid for folding later */
1030             cm |= quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1031                   lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
1032                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1033          }
1034       }
1035
1036    } else {
1037       /* This is the basic no-split case */
1038       q = bits2pulses(m, i, LM, b);
1039       curr_bits = pulses2bits(m, i, LM, q);
1040       *remaining_bits -= curr_bits;
1041
1042       /* Ensures we can never bust the budget */
1043       while (*remaining_bits < 0 && q > 0)
1044       {
1045          *remaining_bits += curr_bits;
1046          q--;
1047          curr_bits = pulses2bits(m, i, LM, q);
1048          *remaining_bits -= curr_bits;
1049       }
1050
1051       if (q!=0)
1052       {
1053          int K = get_pulses(q);
1054
1055          /* Finally do the actual quantization */
1056          if (encode)
1057             cm = alg_quant(X, N, K, spread, B, resynth, ec, gain);
1058          else
1059             cm = alg_unquant(X, N, K, spread, B, ec, gain);
1060       } else {
1061          /* If there's no pulse, fill the band anyway */
1062          int j;
1063          if (resynth)
1064          {
1065             unsigned cm_mask;
1066             /*B can be as large as 16, so this shift might overflow an int on a
1067                16-bit platform; use a long to get defined behavior.*/
1068             cm_mask = (unsigned)(1UL<<B)-1;
1069             fill &= cm_mask;
1070             if (!fill)
1071             {
1072                for (j=0;j<N;j++)
1073                   X[j] = 0;
1074             } else {
1075                if (lowband == NULL)
1076                {
1077                   /* Noise */
1078                   for (j=0;j<N;j++)
1079                   {
1080                      *seed = lcg_rand(*seed);
1081                      X[j] = (celt_int32)(*seed)>>20;
1082                   }
1083                   cm = cm_mask;
1084                } else {
1085                   /* Folded spectrum */
1086                   for (j=0;j<N;j++)
1087                   {
1088                      celt_word16 tmp;
1089                      *seed = lcg_rand(*seed);
1090                      /* About 48 dB below the "normal" folding level */
1091                      tmp = QCONST16(1.0f/256, 10);
1092                      tmp = (*seed)&0x8000 ? tmp : -tmp;
1093                      X[j] = lowband[j]+tmp;
1094                   }
1095                   cm = fill;
1096                }
1097                renormalise_vector(X, N, gain);
1098             }
1099          }
1100       }
1101    }
1102
1103    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1104    if (resynth)
1105    {
1106       if (stereo)
1107       {
1108          if (N!=2)
1109             stereo_merge(X, Y, mid, N);
1110          if (inv)
1111          {
1112             int j;
1113             for (j=0;j<N;j++)
1114                Y[j] = -Y[j];
1115          }
1116       } else if (level == 0)
1117       {
1118          int k;
1119
1120          /* Undo the sample reorganization going from time order to frequency order */
1121          if (B0>1)
1122             interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1123
1124          /* Undo time-freq changes that we did earlier */
1125          N_B = N_B0;
1126          B = B0;
1127          for (k=0;k<time_divide;k++)
1128          {
1129             B >>= 1;
1130             N_B <<= 1;
1131             cm |= cm>>B;
1132             haar1(X, N_B, B);
1133          }
1134
1135          for (k=0;k<recombine;k++)
1136          {
1137             static const unsigned char bit_deinterleave_table[16]={
1138               0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1139               0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1140             };
1141             cm = bit_deinterleave_table[cm];
1142             haar1(X, N0>>k, 1<<k);
1143          }
1144          B<<=recombine;
1145
1146          /* Scale output for later folding */
1147          if (lowband_out)
1148          {
1149             int j;
1150             celt_word16 n;
1151             n = celt_sqrt(SHL32(EXTEND32(N0),22));
1152             for (j=0;j<N0;j++)
1153                lowband_out[j] = MULT16_16_Q15(n,X[j]);
1154          }
1155          cm &= (1<<B)-1;
1156       }
1157    }
1158    return cm;
1159 }
1160
1161 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1162       celt_norm *_X, celt_norm *_Y, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1163       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res, int resynth,
1164       celt_int32 total_bits, celt_int32 balance, ec_ctx *ec, int LM, int codedBands, ec_uint32 *seed)
1165 {
1166    int i;
1167    celt_int32 remaining_bits;
1168    const celt_int16 * restrict eBands = m->eBands;
1169    celt_norm * restrict norm, * restrict norm2;
1170    VARDECL(celt_norm, _norm);
1171    VARDECL(celt_norm, lowband_scratch);
1172    int B;
1173    int M;
1174    int lowband_offset;
1175    int update_lowband = 1;
1176    int C = _Y != NULL ? 2 : 1;
1177    SAVE_STACK;
1178
1179    M = 1<<LM;
1180    B = shortBlocks ? M : 1;
1181    ALLOC(_norm, C*M*eBands[m->nbEBands], celt_norm);
1182    ALLOC(lowband_scratch, M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]), celt_norm);
1183    norm = _norm;
1184    norm2 = norm + M*eBands[m->nbEBands];
1185
1186    lowband_offset = 0;
1187    for (i=start;i<end;i++)
1188    {
1189       celt_int32 tell;
1190       int b;
1191       int N;
1192       celt_int32 curr_balance;
1193       int effective_lowband=-1;
1194       celt_norm * restrict X, * restrict Y;
1195       int tf_change=0;
1196       unsigned x_cm;
1197       unsigned y_cm;
1198
1199       X = _X+M*eBands[i];
1200       if (_Y!=NULL)
1201          Y = _Y+M*eBands[i];
1202       else
1203          Y = NULL;
1204       N = M*eBands[i+1]-M*eBands[i];
1205       tell = ec_tell_frac(ec);
1206
1207       /* Compute how many bits we want to allocate to this band */
1208       if (i != start)
1209          balance -= tell;
1210       remaining_bits = total_bits-tell-1;
1211       if (i <= codedBands-1)
1212       {
1213          curr_balance = balance / IMIN(3, codedBands-i);
1214          b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1215       } else {
1216          b = 0;
1217       }
1218
1219       if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1220             lowband_offset = i;
1221
1222       tf_change = tf_res[i];
1223       if (i>=m->effEBands)
1224       {
1225          X=norm;
1226          if (_Y!=NULL)
1227             Y = norm;
1228       }
1229
1230       /* Get a conservative estimate of the collapse_mask's for the bands we're
1231           going to be folding from. */
1232       if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1233       {
1234          int fold_start;
1235          int fold_end;
1236          int fold_i;
1237          /* This ensures we never repeat spectral content within one band */
1238          effective_lowband = IMAX(M*eBands[start], M*eBands[lowband_offset]-N);
1239          fold_start = lowband_offset;
1240          while(M*eBands[--fold_start] > effective_lowband);
1241          fold_end = lowband_offset-1;
1242          while(M*eBands[++fold_end] < effective_lowband+N);
1243          x_cm = y_cm = 0;
1244          fold_i = fold_start; do {
1245            x_cm |= collapse_masks[fold_i*C+0];
1246            y_cm |= collapse_masks[fold_i*C+C-1];
1247          } while (++fold_i<fold_end);
1248       }
1249       /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1250           always) be non-zero.*/
1251       else
1252          x_cm = y_cm = (1<<B)-1;
1253
1254       if (dual_stereo && i==intensity)
1255       {
1256          int j;
1257
1258          /* Switch off dual stereo to do intensity */
1259          dual_stereo = 0;
1260          for (j=M*eBands[start];j<M*eBands[i];j++)
1261             norm[j] = HALF32(norm[j]+norm2[j]);
1262       }
1263       if (dual_stereo)
1264       {
1265          x_cm = quant_band(encode, m, i, X, NULL, N, b/2, spread, B, intensity, tf_change,
1266                effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1267                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm);
1268          y_cm = quant_band(encode, m, i, Y, NULL, N, b/2, spread, B, intensity, tf_change,
1269                effective_lowband != -1 ? norm2+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1270                norm2+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, y_cm);
1271       } else {
1272          x_cm = quant_band(encode, m, i, X, Y, N, b, spread, B, intensity, tf_change,
1273                effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1274                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm|y_cm);
1275          y_cm = x_cm;
1276       }
1277       collapse_masks[i*C+0] = (unsigned char)x_cm;
1278       collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1279       balance += pulses[i] + tell;
1280
1281       /* Update the folding position only as long as we have 1 bit/sample depth */
1282       update_lowband = b>(N<<BITRES);
1283    }
1284    RESTORE_STACK;
1285 }
1286