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