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