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