Re-enables test cases, including the rotation test
[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    return decision;
481 }
482
483 #ifdef MEASURE_NORM_MSE
484
485 float MSE[30] = {0};
486 int nbMSEBands = 0;
487 int MSECount[30] = {0};
488
489 void dump_norm_mse(void)
490 {
491    int i;
492    for (i=0;i<nbMSEBands;i++)
493    {
494       printf ("%g ", MSE[i]/MSECount[i]);
495    }
496    printf ("\n");
497 }
498
499 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
500 {
501    static int init = 0;
502    int i;
503    if (!init)
504    {
505       atexit(dump_norm_mse);
506       init = 1;
507    }
508    for (i=0;i<m->nbEBands;i++)
509    {
510       int j;
511       int c;
512       float g;
513       if (bandE0[i]<10 || (C==2 && bandE0[i+m->nbEBands]<1))
514          continue;
515       c=0; do {
516          g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
517          for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
518             MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
519       } while (++c<C);
520       MSECount[i]+=C;
521    }
522    nbMSEBands = m->nbEBands;
523 }
524
525 #endif
526
527 /* Indexing table for converting from natural Hadamard to ordery Hadamard
528    This is essentially a bit-reversed Gray, on top of which we've added
529    an inversion of the order because we want the DC at the end rather than
530    the beginning. The lines are for N=2, 4, 8, 16 */
531 static const int ordery_table[] = {
532        1,  0,
533        3,  0,  2,  1,
534        7,  0,  4,  3,  6,  1,  5,  2,
535       15,  0,  8,  7, 12,  3, 11,  4, 14,  1,  9,  6, 13,  2, 10,  5,
536 };
537
538 static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
539 {
540    int i,j;
541    VARDECL(celt_norm, tmp);
542    int N;
543    SAVE_STACK;
544    N = N0*stride;
545    ALLOC(tmp, N, celt_norm);
546    if (hadamard)
547    {
548       const int *ordery = ordery_table+stride-2;
549       for (i=0;i<stride;i++)
550       {
551          for (j=0;j<N0;j++)
552             tmp[ordery[i]*N0+j] = X[j*stride+i];
553       }
554    } else {
555       for (i=0;i<stride;i++)
556          for (j=0;j<N0;j++)
557             tmp[i*N0+j] = X[j*stride+i];
558    }
559    for (j=0;j<N;j++)
560       X[j] = tmp[j];
561    RESTORE_STACK;
562 }
563
564 static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
565 {
566    int i,j;
567    VARDECL(celt_norm, tmp);
568    int N;
569    SAVE_STACK;
570    N = N0*stride;
571    ALLOC(tmp, N, celt_norm);
572    if (hadamard)
573    {
574       const int *ordery = ordery_table+stride-2;
575       for (i=0;i<stride;i++)
576          for (j=0;j<N0;j++)
577             tmp[j*stride+i] = X[ordery[i]*N0+j];
578    } else {
579       for (i=0;i<stride;i++)
580          for (j=0;j<N0;j++)
581             tmp[j*stride+i] = X[i*N0+j];
582    }
583    for (j=0;j<N;j++)
584       X[j] = tmp[j];
585    RESTORE_STACK;
586 }
587
588 void haar1(celt_norm *X, int N0, int stride)
589 {
590    int i, j;
591    N0 >>= 1;
592    for (i=0;i<stride;i++)
593       for (j=0;j<N0;j++)
594       {
595          celt_norm tmp1, tmp2;
596          tmp1 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*2*j+i]);
597          tmp2 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*(2*j+1)+i]);
598          X[stride*2*j+i] = tmp1 + tmp2;
599          X[stride*(2*j+1)+i] = tmp1 - tmp2;
600       }
601 }
602
603 static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo)
604 {
605    static const opus_int16 exp2_table8[8] =
606       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
607    int qn, qb;
608    int N2 = 2*N-1;
609    if (stereo && N==2)
610       N2--;
611    /* The upper limit ensures that in a stereo split with itheta==16384, we'll
612        always have enough bits left over to code at least one pulse in the
613        side; otherwise it would collapse, since it doesn't get folded. */
614    qb = IMIN(b-pulse_cap-(4<<BITRES), (b+N2*offset)/N2);
615
616    qb = IMIN(8<<BITRES, qb);
617
618    if (qb<(1<<BITRES>>1)) {
619       qn = 1;
620    } else {
621       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
622       qn = (qn+1)>>1<<1;
623    }
624    celt_assert(qn <= 256);
625    return qn;
626 }
627
628 /* This function is responsible for encoding and decoding a band for both
629    the mono and stereo case. Even in the mono case, it can split the band
630    in two and transmit the energy difference with the two half-bands. It
631    can be called recursively so bands can end up being split in 8 parts. */
632 static unsigned quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
633       int N, int b, int spread, int B, int intensity, int tf_change, celt_norm *lowband, int resynth, ec_ctx *ec,
634       opus_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level,
635       opus_uint32 *seed, opus_val16 gain, celt_norm *lowband_scratch, int fill)
636 {
637    const unsigned char *cache;
638    int q;
639    int curr_bits;
640    int stereo, split;
641    int imid=0, iside=0;
642    int N0=N;
643    int N_B=N;
644    int N_B0;
645    int B0=B;
646    int time_divide=0;
647    int recombine=0;
648    int inv = 0;
649    opus_val16 mid=0, side=0;
650    int longBlocks;
651    unsigned cm=0;
652
653    longBlocks = B0==1;
654
655    N_B /= B;
656    N_B0 = N_B;
657
658    split = stereo = Y != NULL;
659
660    /* Special case for one sample */
661    if (N==1)
662    {
663       int c;
664       celt_norm *x = X;
665       c=0; do {
666          int sign=0;
667          if (*remaining_bits>=1<<BITRES)
668          {
669             if (encode)
670             {
671                sign = x[0]<0;
672                ec_enc_bits(ec, sign, 1);
673             } else {
674                sign = ec_dec_bits(ec, 1);
675             }
676             *remaining_bits -= 1<<BITRES;
677             b-=1<<BITRES;
678          }
679          if (resynth)
680             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
681          x = Y;
682       } while (++c<1+stereo);
683       if (lowband_out)
684          lowband_out[0] = SHR16(X[0],4);
685       return 1;
686    }
687
688    if (!stereo && level == 0)
689    {
690       int k;
691       if (tf_change>0)
692          recombine = tf_change;
693       /* Band recombining to increase frequency resolution */
694
695       if (lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
696       {
697          int j;
698          for (j=0;j<N;j++)
699             lowband_scratch[j] = lowband[j];
700          lowband = lowband_scratch;
701       }
702
703       for (k=0;k<recombine;k++)
704       {
705          static const unsigned char bit_interleave_table[16]={
706            0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
707          };
708          if (encode)
709             haar1(X, N>>k, 1<<k);
710          if (lowband)
711             haar1(lowband, N>>k, 1<<k);
712          fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
713       }
714       B>>=recombine;
715       N_B<<=recombine;
716
717       /* Increasing the time resolution */
718       while ((N_B&1) == 0 && tf_change<0)
719       {
720          if (encode)
721             haar1(X, N_B, B);
722          if (lowband)
723             haar1(lowband, N_B, B);
724          fill |= fill<<B;
725          B <<= 1;
726          N_B >>= 1;
727          time_divide++;
728          tf_change++;
729       }
730       B0=B;
731       N_B0 = N_B;
732
733       /* Reorganize the samples in time order instead of frequency order */
734       if (B0>1)
735       {
736          if (encode)
737             deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
738          if (lowband)
739             deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
740       }
741    }
742
743    /* If we need 1.5 more bit than we can produce, split the band in two. */
744    cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
745    if (!stereo && LM != -1 && b > cache[cache[0]]+12 && N>2)
746    {
747       if (LM>0 || (N&1)==0)
748       {
749          N >>= 1;
750          Y = X+N;
751          split = 1;
752          LM -= 1;
753          if (B==1)
754             fill = fill&1|fill<<1;
755          B = (B+1)>>1;
756       }
757    }
758
759    if (split)
760    {
761       int qn;
762       int itheta=0;
763       int mbits, sbits, delta;
764       int qalloc;
765       int pulse_cap;
766       int offset;
767       int orig_fill;
768       opus_int32 tell;
769
770       /* Decide on the resolution to give to the split parameter theta */
771       pulse_cap = m->logN[i]+(LM<<BITRES);
772       offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
773       qn = compute_qn(N, b, offset, pulse_cap, stereo);
774       if (stereo && i>=intensity)
775          qn = 1;
776       if (encode)
777       {
778          /* theta is the atan() of the ratio between the (normalized)
779             side and mid. With just that parameter, we can re-scale both
780             mid and side because we know that 1) they have unit norm and
781             2) they are orthogonal. */
782          itheta = stereo_itheta(X, Y, stereo, N);
783       }
784       tell = ec_tell_frac(ec);
785       if (qn!=1)
786       {
787          if (encode)
788             itheta = (itheta*qn+8192)>>14;
789
790          /* Entropy coding of the angle. We use a uniform pdf for the
791             time split, a step for stereo, and a triangular one for the rest. */
792          if (stereo && N>2)
793          {
794             int p0 = 3;
795             int x = itheta;
796             int x0 = qn/2;
797             int ft = p0*(x0+1) + x0;
798             /* Use a probability of p0 up to itheta=8192 and then use 1 after */
799             if (encode)
800             {
801                ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
802             } else {
803                int fs;
804                fs=ec_decode(ec,ft);
805                if (fs<(x0+1)*p0)
806                   x=fs/p0;
807                else
808                   x=x0+1+(fs-(x0+1)*p0);
809                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);
810                itheta = x;
811             }
812          } else if (B0>1 || stereo) {
813             /* Uniform pdf */
814             if (encode)
815                ec_enc_uint(ec, itheta, qn+1);
816             else
817                itheta = ec_dec_uint(ec, qn+1);
818          } else {
819             int fs=1, ft;
820             ft = ((qn>>1)+1)*((qn>>1)+1);
821             if (encode)
822             {
823                int fl;
824
825                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
826                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
827                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
828
829                ec_encode(ec, fl, fl+fs, ft);
830             } else {
831                /* Triangular pdf */
832                int fl=0;
833                int fm;
834                fm = ec_decode(ec, ft);
835
836                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
837                {
838                   itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;
839                   fs = itheta + 1;
840                   fl = itheta*(itheta + 1)>>1;
841                }
842                else
843                {
844                   itheta = (2*(qn + 1)
845                    - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;
846                   fs = qn + 1 - itheta;
847                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
848                }
849
850                ec_dec_update(ec, fl, fl+fs, ft);
851             }
852          }
853          itheta = (opus_int32)itheta*16384/qn;
854          if (encode && stereo)
855          {
856             if (itheta==0)
857                intensity_stereo(m, X, Y, bandE, i, N);
858             else
859                stereo_split(X, Y, N);
860          }
861          /* TODO: Renormalising X and Y *may* help fixed-point a bit at very high rate.
862                   Let's do that at higher complexity */
863       } else if (stereo) {
864          if (encode)
865          {
866             inv = itheta > 8192;
867             if (inv)
868             {
869                int j;
870                for (j=0;j<N;j++)
871                   Y[j] = -Y[j];
872             }
873             intensity_stereo(m, X, Y, bandE, i, N);
874          }
875          if (b>2<<BITRES && *remaining_bits > 2<<BITRES)
876          {
877             if (encode)
878                ec_enc_bit_logp(ec, inv, 2);
879             else
880                inv = ec_dec_bit_logp(ec, 2);
881          } else
882             inv = 0;
883          itheta = 0;
884       }
885       qalloc = ec_tell_frac(ec) - tell;
886       b -= qalloc;
887
888       orig_fill = fill;
889       if (itheta == 0)
890       {
891          imid = 32767;
892          iside = 0;
893          fill &= (1<<B)-1;
894          delta = -16384;
895       } else if (itheta == 16384)
896       {
897          imid = 0;
898          iside = 32767;
899          fill &= (1<<B)-1<<B;
900          delta = 16384;
901       } else {
902          imid = bitexact_cos(itheta);
903          iside = bitexact_cos(16384-itheta);
904          /* This is the mid vs side allocation that minimizes squared error
905             in that band. */
906          delta = FRAC_MUL16(N-1<<7,bitexact_log2tan(iside,imid));
907       }
908
909 #ifdef FIXED_POINT
910       mid = imid;
911       side = iside;
912 #else
913       mid = (1.f/32768)*imid;
914       side = (1.f/32768)*iside;
915 #endif
916
917       /* This is a special case for N=2 that only works for stereo and takes
918          advantage of the fact that mid and side are orthogonal to encode
919          the side with just one bit. */
920       if (N==2 && stereo)
921       {
922          int c;
923          int sign=0;
924          celt_norm *x2, *y2;
925          mbits = b;
926          sbits = 0;
927          /* Only need one bit for the side */
928          if (itheta != 0 && itheta != 16384)
929             sbits = 1<<BITRES;
930          mbits -= sbits;
931          c = itheta > 8192;
932          *remaining_bits -= qalloc+sbits;
933
934          x2 = c ? Y : X;
935          y2 = c ? X : Y;
936          if (sbits)
937          {
938             if (encode)
939             {
940                /* Here we only need to encode a sign for the side */
941                sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
942                ec_enc_bits(ec, sign, 1);
943             } else {
944                sign = ec_dec_bits(ec, 1);
945             }
946          }
947          sign = 1-2*sign;
948          /* We use orig_fill here because we want to fold the side, but if
949              itheta==16384, we'll have cleared the low bits of fill. */
950          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);
951          /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
952              and there's no need to worry about mixing with the other channel. */
953          y2[0] = -sign*x2[1];
954          y2[1] = sign*x2[0];
955          if (resynth)
956          {
957             celt_norm tmp;
958             X[0] = MULT16_16_Q15(mid, X[0]);
959             X[1] = MULT16_16_Q15(mid, X[1]);
960             Y[0] = MULT16_16_Q15(side, Y[0]);
961             Y[1] = MULT16_16_Q15(side, Y[1]);
962             tmp = X[0];
963             X[0] = SUB16(tmp,Y[0]);
964             Y[0] = ADD16(tmp,Y[0]);
965             tmp = X[1];
966             X[1] = SUB16(tmp,Y[1]);
967             Y[1] = ADD16(tmp,Y[1]);
968          }
969       } else {
970          /* "Normal" split code */
971          celt_norm *next_lowband2=NULL;
972          celt_norm *next_lowband_out1=NULL;
973          int next_level=0;
974          opus_int32 rebalance;
975
976          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
977          if (B0>1 && !stereo && (itheta&0x3fff))
978          {
979             if (itheta > 8192)
980                /* Rough approximation for pre-echo masking */
981                delta -= delta>>(4-LM);
982             else
983                /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
984                delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
985          }
986          mbits = IMAX(0, IMIN(b, (b-delta)/2));
987          sbits = b-mbits;
988          *remaining_bits -= qalloc;
989
990          if (lowband && !stereo)
991             next_lowband2 = lowband+N; /* >32-bit split case */
992
993          /* Only stereo needs to pass on lowband_out. Otherwise, it's
994             handled at the end */
995          if (stereo)
996             next_lowband_out1 = lowband_out;
997          else
998             next_level = level+1;
999
1000          rebalance = *remaining_bits;
1001          if (mbits >= sbits)
1002          {
1003             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1004                mid for folding later */
1005             cm = quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1006                   lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
1007                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1008             rebalance = mbits - (rebalance-*remaining_bits);
1009             if (rebalance > 3<<BITRES && itheta!=0)
1010                sbits += rebalance - (3<<BITRES);
1011
1012             /* For a stereo split, the high bits of fill are always zero, so no
1013                folding will be done to the side. */
1014             cm |= quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1015                   next_lowband2, resynth, ec, remaining_bits, LM, NULL,
1016                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<(B0>>1&stereo-1);
1017          } else {
1018             /* For a stereo split, the high bits of fill are always zero, so no
1019                folding will be done to the side. */
1020             cm = quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1021                   next_lowband2, resynth, ec, remaining_bits, LM, NULL,
1022                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<(B0>>1&stereo-1);
1023             rebalance = sbits - (rebalance-*remaining_bits);
1024             if (rebalance > 3<<BITRES && itheta!=16384)
1025                mbits += rebalance - (3<<BITRES);
1026             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1027                mid for folding later */
1028             cm |= quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1029                   lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
1030                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1031          }
1032       }
1033
1034    } else {
1035       /* This is the basic no-split case */
1036       q = bits2pulses(m, i, LM, b);
1037       curr_bits = pulses2bits(m, i, LM, q);
1038       *remaining_bits -= curr_bits;
1039
1040       /* Ensures we can never bust the budget */
1041       while (*remaining_bits < 0 && q > 0)
1042       {
1043          *remaining_bits += curr_bits;
1044          q--;
1045          curr_bits = pulses2bits(m, i, LM, q);
1046          *remaining_bits -= curr_bits;
1047       }
1048
1049       if (q!=0)
1050       {
1051          int K = get_pulses(q);
1052
1053          /* Finally do the actual quantization */
1054          if (encode)
1055             cm = alg_quant(X, N, K, spread, B, resynth, ec, gain);
1056          else
1057             cm = alg_unquant(X, N, K, spread, B, ec, gain);
1058       } else {
1059          /* If there's no pulse, fill the band anyway */
1060          int j;
1061          if (resynth)
1062          {
1063             unsigned cm_mask;
1064             /*B can be as large as 16, so this shift might overflow an int on a
1065                16-bit platform; use a long to get defined behavior.*/
1066             cm_mask = (unsigned)(1UL<<B)-1;
1067             fill &= cm_mask;
1068             if (!fill)
1069             {
1070                for (j=0;j<N;j++)
1071                   X[j] = 0;
1072             } else {
1073                if (lowband == NULL)
1074                {
1075                   /* Noise */
1076                   for (j=0;j<N;j++)
1077                   {
1078                      *seed = lcg_rand(*seed);
1079                      X[j] = (opus_int32)(*seed)>>20;
1080                   }
1081                   cm = cm_mask;
1082                } else {
1083                   /* Folded spectrum */
1084                   for (j=0;j<N;j++)
1085                   {
1086                      opus_val16 tmp;
1087                      *seed = lcg_rand(*seed);
1088                      /* About 48 dB below the "normal" folding level */
1089                      tmp = QCONST16(1.0f/256, 10);
1090                      tmp = (*seed)&0x8000 ? tmp : -tmp;
1091                      X[j] = lowband[j]+tmp;
1092                   }
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
1144          /* Scale output for later folding */
1145          if (lowband_out)
1146          {
1147             int j;
1148             opus_val16 n;
1149             n = celt_sqrt(SHL32(EXTEND32(N0),22));
1150             for (j=0;j<N0;j++)
1151                lowband_out[j] = MULT16_16_Q15(n,X[j]);
1152          }
1153          cm &= (1<<B)-1;
1154       }
1155    }
1156    return cm;
1157 }
1158
1159 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1160       celt_norm *_X, celt_norm *_Y, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1161       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res, int resynth,
1162       opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed)
1163 {
1164    int i;
1165    opus_int32 remaining_bits;
1166    const opus_int16 * restrict eBands = m->eBands;
1167    celt_norm * restrict norm, * restrict norm2;
1168    VARDECL(celt_norm, _norm);
1169    VARDECL(celt_norm, lowband_scratch);
1170    int B;
1171    int M;
1172    int lowband_offset;
1173    int update_lowband = 1;
1174    int C = _Y != NULL ? 2 : 1;
1175    SAVE_STACK;
1176
1177    M = 1<<LM;
1178    B = shortBlocks ? M : 1;
1179    ALLOC(_norm, C*M*eBands[m->nbEBands], celt_norm);
1180    ALLOC(lowband_scratch, M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]), celt_norm);
1181    norm = _norm;
1182    norm2 = norm + M*eBands[m->nbEBands];
1183
1184    lowband_offset = 0;
1185    for (i=start;i<end;i++)
1186    {
1187       opus_int32 tell;
1188       int b;
1189       int N;
1190       opus_int32 curr_balance;
1191       int effective_lowband=-1;
1192       celt_norm * restrict X, * restrict Y;
1193       int tf_change=0;
1194       unsigned x_cm;
1195       unsigned y_cm;
1196
1197       X = _X+M*eBands[i];
1198       if (_Y!=NULL)
1199          Y = _Y+M*eBands[i];
1200       else
1201          Y = NULL;
1202       N = M*eBands[i+1]-M*eBands[i];
1203       tell = ec_tell_frac(ec);
1204
1205       /* Compute how many bits we want to allocate to this band */
1206       if (i != start)
1207          balance -= tell;
1208       remaining_bits = total_bits-tell-1;
1209       if (i <= codedBands-1)
1210       {
1211          curr_balance = balance / IMIN(3, codedBands-i);
1212          b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1213       } else {
1214          b = 0;
1215       }
1216
1217       if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1218             lowband_offset = i;
1219
1220       tf_change = tf_res[i];
1221       if (i>=m->effEBands)
1222       {
1223          X=norm;
1224          if (_Y!=NULL)
1225             Y = norm;
1226       }
1227
1228       /* Get a conservative estimate of the collapse_mask's for the bands we're
1229           going to be folding from. */
1230       if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1231       {
1232          int fold_start;
1233          int fold_end;
1234          int fold_i;
1235          /* This ensures we never repeat spectral content within one band */
1236          effective_lowband = IMAX(M*eBands[start], M*eBands[lowband_offset]-N);
1237          fold_start = lowband_offset;
1238          while(M*eBands[--fold_start] > effective_lowband);
1239          fold_end = lowband_offset-1;
1240          while(M*eBands[++fold_end] < effective_lowband+N);
1241          x_cm = y_cm = 0;
1242          fold_i = fold_start; do {
1243            x_cm |= collapse_masks[fold_i*C+0];
1244            y_cm |= collapse_masks[fold_i*C+C-1];
1245          } while (++fold_i<fold_end);
1246       }
1247       /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1248           always) be non-zero.*/
1249       else
1250          x_cm = y_cm = (1<<B)-1;
1251
1252       if (dual_stereo && i==intensity)
1253       {
1254          int j;
1255
1256          /* Switch off dual stereo to do intensity */
1257          dual_stereo = 0;
1258          for (j=M*eBands[start];j<M*eBands[i];j++)
1259             norm[j] = HALF32(norm[j]+norm2[j]);
1260       }
1261       if (dual_stereo)
1262       {
1263          x_cm = quant_band(encode, m, i, X, NULL, N, b/2, spread, B, intensity, tf_change,
1264                effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1265                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm);
1266          y_cm = quant_band(encode, m, i, Y, NULL, N, b/2, spread, B, intensity, tf_change,
1267                effective_lowband != -1 ? norm2+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1268                norm2+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, y_cm);
1269       } else {
1270          x_cm = quant_band(encode, m, i, X, Y, N, b, spread, B, intensity, tf_change,
1271                effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1272                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm|y_cm);
1273          y_cm = x_cm;
1274       }
1275       collapse_masks[i*C+0] = (unsigned char)x_cm;
1276       collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1277       balance += pulses[i] + tell;
1278
1279       /* Update the folding position only as long as we have 1 bit/sample depth */
1280       update_lowband = b>(N<<BITRES);
1281    }
1282    RESTORE_STACK;
1283 }
1284