First step in quant_band() cleanup: N=1 case.
[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 static unsigned quant_band_n1(int encode, celt_norm *X, celt_norm *Y, int b,
652       opus_int32 *remaining_bits, ec_ctx *ec, celt_norm *lowband_out)
653 {
654 #ifdef RESYNTH
655    int resynth = 1;
656 #else
657    int resynth = !encode;
658 #endif
659    int c;
660    int stereo;
661    celt_norm *x = X;
662    stereo = Y != NULL;
663    c=0; do {
664       int sign=0;
665       if (*remaining_bits>=1<<BITRES)
666       {
667          if (encode)
668          {
669             sign = x[0]<0;
670             ec_enc_bits(ec, sign, 1);
671          } else {
672             sign = ec_dec_bits(ec, 1);
673          }
674          *remaining_bits -= 1<<BITRES;
675          b-=1<<BITRES;
676       }
677       if (resynth)
678          x[0] = sign ? -NORM_SCALING : NORM_SCALING;
679       x = Y;
680    } while (++c<1+stereo);
681    if (lowband_out)
682       lowband_out[0] = SHR16(X[0],4);
683    return 1;
684 }
685
686 /* This function is responsible for encoding and decoding a band for both
687    the mono and stereo case. Even in the mono case, it can split the band
688    in two and transmit the energy difference with the two half-bands. It
689    can be called recursively so bands can end up being split in 8 parts. */
690 static unsigned quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
691       int N, int b, int spread, int B, int intensity, int tf_change, celt_norm *lowband, ec_ctx *ec,
692       opus_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level,
693       opus_uint32 *seed, opus_val16 gain, celt_norm *lowband_scratch, int fill)
694 {
695    const unsigned char *cache;
696    int q;
697    int curr_bits;
698    int stereo, split;
699    int imid=0, iside=0;
700    int N0=N;
701    int N_B=N;
702    int N_B0;
703    int B0=B;
704    int time_divide=0;
705    int recombine=0;
706    int inv = 0;
707    opus_val16 mid=0, side=0;
708    int longBlocks;
709    unsigned cm=0;
710 #ifdef RESYNTH
711    int resynth = 1;
712 #else
713    int resynth = !encode;
714 #endif
715
716    longBlocks = B0==1;
717
718    N_B /= B;
719    N_B0 = N_B;
720
721    split = stereo = Y != NULL;
722
723    /* Special case for one sample */
724    if (N==1)
725    {
726       return quant_band_n1(encode, X, Y, b, remaining_bits, ec, lowband_out);
727    }
728
729    if (!stereo && level == 0)
730    {
731       int k;
732       if (tf_change>0)
733          recombine = tf_change;
734       /* Band recombining to increase frequency resolution */
735
736       if (lowband_scratch && lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
737       {
738          int j;
739          for (j=0;j<N;j++)
740             lowband_scratch[j] = lowband[j];
741          lowband = lowband_scratch;
742       }
743
744       for (k=0;k<recombine;k++)
745       {
746          static const unsigned char bit_interleave_table[16]={
747            0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
748          };
749          if (encode)
750             haar1(X, N>>k, 1<<k);
751          if (lowband)
752             haar1(lowband, N>>k, 1<<k);
753          fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
754       }
755       B>>=recombine;
756       N_B<<=recombine;
757
758       /* Increasing the time resolution */
759       while ((N_B&1) == 0 && tf_change<0)
760       {
761          if (encode)
762             haar1(X, N_B, B);
763          if (lowband)
764             haar1(lowband, N_B, B);
765          fill |= fill<<B;
766          B <<= 1;
767          N_B >>= 1;
768          time_divide++;
769          tf_change++;
770       }
771       B0=B;
772       N_B0 = N_B;
773
774       /* Reorganize the samples in time order instead of frequency order */
775       if (B0>1)
776       {
777          if (encode)
778             deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
779          if (lowband)
780             deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
781       }
782    }
783
784    /* If we need 1.5 more bit than we can produce, split the band in two. */
785    cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
786    if (!stereo && LM != -1 && b > cache[cache[0]]+12 && N>2)
787    {
788       N >>= 1;
789       Y = X+N;
790       split = 1;
791       LM -= 1;
792       if (B==1)
793          fill = (fill&1)|(fill<<1);
794       B = (B+1)>>1;
795    }
796
797    if (split)
798    {
799       int qn;
800       int itheta=0;
801       int mbits, sbits, delta;
802       int qalloc;
803       int pulse_cap;
804       int offset;
805       int orig_fill;
806       opus_int32 tell;
807
808       /* Decide on the resolution to give to the split parameter theta */
809       pulse_cap = m->logN[i]+LM*(1<<BITRES);
810       offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
811       qn = compute_qn(N, b, offset, pulse_cap, stereo);
812       if (stereo && i>=intensity)
813          qn = 1;
814       if (encode)
815       {
816          /* theta is the atan() of the ratio between the (normalized)
817             side and mid. With just that parameter, we can re-scale both
818             mid and side because we know that 1) they have unit norm and
819             2) they are orthogonal. */
820          itheta = stereo_itheta(X, Y, stereo, N);
821       }
822       tell = ec_tell_frac(ec);
823       if (qn!=1)
824       {
825          if (encode)
826             itheta = (itheta*qn+8192)>>14;
827
828          /* Entropy coding of the angle. We use a uniform pdf for the
829             time split, a step for stereo, and a triangular one for the rest. */
830          if (stereo && N>2)
831          {
832             int p0 = 3;
833             int x = itheta;
834             int x0 = qn/2;
835             int ft = p0*(x0+1) + x0;
836             /* Use a probability of p0 up to itheta=8192 and then use 1 after */
837             if (encode)
838             {
839                ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
840             } else {
841                int fs;
842                fs=ec_decode(ec,ft);
843                if (fs<(x0+1)*p0)
844                   x=fs/p0;
845                else
846                   x=x0+1+(fs-(x0+1)*p0);
847                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);
848                itheta = x;
849             }
850          } else if (B0>1 || stereo) {
851             /* Uniform pdf */
852             if (encode)
853                ec_enc_uint(ec, itheta, qn+1);
854             else
855                itheta = ec_dec_uint(ec, qn+1);
856          } else {
857             int fs=1, ft;
858             ft = ((qn>>1)+1)*((qn>>1)+1);
859             if (encode)
860             {
861                int fl;
862
863                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
864                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
865                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
866
867                ec_encode(ec, fl, fl+fs, ft);
868             } else {
869                /* Triangular pdf */
870                int fl=0;
871                int fm;
872                fm = ec_decode(ec, ft);
873
874                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
875                {
876                   itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;
877                   fs = itheta + 1;
878                   fl = itheta*(itheta + 1)>>1;
879                }
880                else
881                {
882                   itheta = (2*(qn + 1)
883                    - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;
884                   fs = qn + 1 - itheta;
885                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
886                }
887
888                ec_dec_update(ec, fl, fl+fs, ft);
889             }
890          }
891          itheta = (opus_int32)itheta*16384/qn;
892          if (encode && stereo)
893          {
894             if (itheta==0)
895                intensity_stereo(m, X, Y, bandE, i, N);
896             else
897                stereo_split(X, Y, N);
898          }
899          /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.
900                   Let's do that at higher complexity */
901       } else if (stereo) {
902          if (encode)
903          {
904             inv = itheta > 8192;
905             if (inv)
906             {
907                int j;
908                for (j=0;j<N;j++)
909                   Y[j] = -Y[j];
910             }
911             intensity_stereo(m, X, Y, bandE, i, N);
912          }
913          if (b>2<<BITRES && *remaining_bits > 2<<BITRES)
914          {
915             if (encode)
916                ec_enc_bit_logp(ec, inv, 2);
917             else
918                inv = ec_dec_bit_logp(ec, 2);
919          } else
920             inv = 0;
921          itheta = 0;
922       }
923       qalloc = ec_tell_frac(ec) - tell;
924       b -= qalloc;
925
926       orig_fill = fill;
927       if (itheta == 0)
928       {
929          imid = 32767;
930          iside = 0;
931          fill &= (1<<B)-1;
932          delta = -16384;
933       } else if (itheta == 16384)
934       {
935          imid = 0;
936          iside = 32767;
937          fill &= ((1<<B)-1)<<B;
938          delta = 16384;
939       } else {
940          imid = bitexact_cos((opus_int16)itheta);
941          iside = bitexact_cos((opus_int16)(16384-itheta));
942          /* This is the mid vs side allocation that minimizes squared error
943             in that band. */
944          delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid));
945       }
946
947 #ifdef FIXED_POINT
948       mid = imid;
949       side = iside;
950 #else
951       mid = (1.f/32768)*imid;
952       side = (1.f/32768)*iside;
953 #endif
954
955       /* This is a special case for N=2 that only works for stereo and takes
956          advantage of the fact that mid and side are orthogonal to encode
957          the side with just one bit. */
958       if (N==2 && stereo)
959       {
960          int c;
961          int sign=0;
962          celt_norm *x2, *y2;
963          mbits = b;
964          sbits = 0;
965          /* Only need one bit for the side */
966          if (itheta != 0 && itheta != 16384)
967             sbits = 1<<BITRES;
968          mbits -= sbits;
969          c = itheta > 8192;
970          *remaining_bits -= qalloc+sbits;
971
972          x2 = c ? Y : X;
973          y2 = c ? X : Y;
974          if (sbits)
975          {
976             if (encode)
977             {
978                /* Here we only need to encode a sign for the side */
979                sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
980                ec_enc_bits(ec, sign, 1);
981             } else {
982                sign = ec_dec_bits(ec, 1);
983             }
984          }
985          sign = 1-2*sign;
986          /* We use orig_fill here because we want to fold the side, but if
987              itheta==16384, we'll have cleared the low bits of fill. */
988          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);
989          /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
990              and there's no need to worry about mixing with the other channel. */
991          y2[0] = -sign*x2[1];
992          y2[1] = sign*x2[0];
993          if (resynth)
994          {
995             celt_norm tmp;
996             X[0] = MULT16_16_Q15(mid, X[0]);
997             X[1] = MULT16_16_Q15(mid, X[1]);
998             Y[0] = MULT16_16_Q15(side, Y[0]);
999             Y[1] = MULT16_16_Q15(side, Y[1]);
1000             tmp = X[0];
1001             X[0] = SUB16(tmp,Y[0]);
1002             Y[0] = ADD16(tmp,Y[0]);
1003             tmp = X[1];
1004             X[1] = SUB16(tmp,Y[1]);
1005             Y[1] = ADD16(tmp,Y[1]);
1006          }
1007       } else {
1008          /* "Normal" split code */
1009          celt_norm *next_lowband2=NULL;
1010          celt_norm *next_lowband_out1=NULL;
1011          int next_level=0;
1012          opus_int32 rebalance;
1013
1014          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
1015          if (B0>1 && !stereo && (itheta&0x3fff))
1016          {
1017             if (itheta > 8192)
1018                /* Rough approximation for pre-echo masking */
1019                delta -= delta>>(4-LM);
1020             else
1021                /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
1022                delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
1023          }
1024          mbits = IMAX(0, IMIN(b, (b-delta)/2));
1025          sbits = b-mbits;
1026          *remaining_bits -= qalloc;
1027
1028          if (lowband && !stereo)
1029             next_lowband2 = lowband+N; /* >32-bit split case */
1030
1031          /* Only stereo needs to pass on lowband_out. Otherwise, it's
1032             handled at the end */
1033          if (stereo)
1034             next_lowband_out1 = lowband_out;
1035          else
1036             next_level = level+1;
1037
1038          rebalance = *remaining_bits;
1039          if (mbits >= sbits)
1040          {
1041             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1042                mid for folding later */
1043             cm = quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1044                   lowband, ec, remaining_bits, LM, next_lowband_out1,
1045                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1046             rebalance = mbits - (rebalance-*remaining_bits);
1047             if (rebalance > 3<<BITRES && itheta!=0)
1048                sbits += rebalance - (3<<BITRES);
1049
1050             /* For a stereo split, the high bits of fill are always zero, so no
1051                folding will be done to the side. */
1052             cm |= quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1053                   next_lowband2, ec, remaining_bits, LM, NULL,
1054                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<((B0>>1)&(stereo-1));
1055          } else {
1056             /* For a stereo split, the high bits of fill are always zero, so no
1057                folding will be done to the side. */
1058             cm = quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1059                   next_lowband2, ec, remaining_bits, LM, NULL,
1060                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<((B0>>1)&(stereo-1));
1061             rebalance = sbits - (rebalance-*remaining_bits);
1062             if (rebalance > 3<<BITRES && itheta!=16384)
1063                mbits += rebalance - (3<<BITRES);
1064             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1065                mid for folding later */
1066             cm |= quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1067                   lowband, ec, remaining_bits, LM, next_lowband_out1,
1068                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1069          }
1070       }
1071
1072    } else {
1073       /* This is the basic no-split case */
1074       q = bits2pulses(m, i, LM, b);
1075       curr_bits = pulses2bits(m, i, LM, q);
1076       *remaining_bits -= curr_bits;
1077
1078       /* Ensures we can never bust the budget */
1079       while (*remaining_bits < 0 && q > 0)
1080       {
1081          *remaining_bits += curr_bits;
1082          q--;
1083          curr_bits = pulses2bits(m, i, LM, q);
1084          *remaining_bits -= curr_bits;
1085       }
1086
1087       if (q!=0)
1088       {
1089          int K = get_pulses(q);
1090
1091          /* Finally do the actual quantization */
1092          if (encode)
1093          {
1094             cm = alg_quant(X, N, K, spread, B, ec
1095 #ifdef RESYNTH
1096                  , gain
1097 #endif
1098                  );
1099          } else {
1100             cm = alg_unquant(X, N, K, spread, B, ec, gain);
1101          }
1102       } else {
1103          /* If there's no pulse, fill the band anyway */
1104          int j;
1105          if (resynth)
1106          {
1107             unsigned cm_mask;
1108             /*B can be as large as 16, so this shift might overflow an int on a
1109                16-bit platform; use a long to get defined behavior.*/
1110             cm_mask = (unsigned)(1UL<<B)-1;
1111             fill &= cm_mask;
1112             if (!fill)
1113             {
1114                for (j=0;j<N;j++)
1115                   X[j] = 0;
1116             } else {
1117                if (lowband == NULL)
1118                {
1119                   /* Noise */
1120                   for (j=0;j<N;j++)
1121                   {
1122                      *seed = celt_lcg_rand(*seed);
1123                      X[j] = (celt_norm)((opus_int32)*seed>>20);
1124                   }
1125                   cm = cm_mask;
1126                } else {
1127                   /* Folded spectrum */
1128                   for (j=0;j<N;j++)
1129                   {
1130                      opus_val16 tmp;
1131                      *seed = celt_lcg_rand(*seed);
1132                      /* About 48 dB below the "normal" folding level */
1133                      tmp = QCONST16(1.0f/256, 10);
1134                      tmp = (*seed)&0x8000 ? tmp : -tmp;
1135                      X[j] = lowband[j]+tmp;
1136                   }
1137                   cm = fill;
1138                }
1139                renormalise_vector(X, N, gain);
1140             }
1141          }
1142       }
1143    }
1144
1145    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1146    if (resynth)
1147    {
1148       if (stereo)
1149       {
1150          if (N!=2)
1151             stereo_merge(X, Y, mid, N);
1152          if (inv)
1153          {
1154             int j;
1155             for (j=0;j<N;j++)
1156                Y[j] = -Y[j];
1157          }
1158       } else if (level == 0)
1159       {
1160          int k;
1161
1162          /* Undo the sample reorganization going from time order to frequency order */
1163          if (B0>1)
1164             interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1165
1166          /* Undo time-freq changes that we did earlier */
1167          N_B = N_B0;
1168          B = B0;
1169          for (k=0;k<time_divide;k++)
1170          {
1171             B >>= 1;
1172             N_B <<= 1;
1173             cm |= cm>>B;
1174             haar1(X, N_B, B);
1175          }
1176
1177          for (k=0;k<recombine;k++)
1178          {
1179             static const unsigned char bit_deinterleave_table[16]={
1180               0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1181               0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1182             };
1183             cm = bit_deinterleave_table[cm];
1184             haar1(X, N0>>k, 1<<k);
1185          }
1186          B<<=recombine;
1187
1188          /* Scale output for later folding */
1189          if (lowband_out)
1190          {
1191             int j;
1192             opus_val16 n;
1193             n = celt_sqrt(SHL32(EXTEND32(N0),22));
1194             for (j=0;j<N0;j++)
1195                lowband_out[j] = MULT16_16_Q15(n,X[j]);
1196          }
1197          cm &= (1<<B)-1;
1198       }
1199    }
1200    return cm;
1201 }
1202
1203 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1204       celt_norm *X_, celt_norm *Y_, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1205       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res,
1206       opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed)
1207 {
1208    int i;
1209    opus_int32 remaining_bits;
1210    const opus_int16 * OPUS_RESTRICT eBands = m->eBands;
1211    celt_norm * OPUS_RESTRICT norm, * OPUS_RESTRICT norm2;
1212    VARDECL(celt_norm, _norm);
1213    celt_norm *lowband_scratch;
1214    int B;
1215    int M;
1216    int lowband_offset;
1217    int update_lowband = 1;
1218    int C = Y_ != NULL ? 2 : 1;
1219    int norm_offset;
1220 #ifdef RESYNTH
1221    int resynth = 1;
1222 #else
1223    int resynth = !encode;
1224 #endif
1225    SAVE_STACK;
1226
1227    M = 1<<LM;
1228    B = shortBlocks ? M : 1;
1229    norm_offset = M*eBands[start];
1230    /* No need to allocate norm for the last band because we don't need an
1231       output in that band */
1232    ALLOC(_norm, C*(M*eBands[m->nbEBands-1]-norm_offset), celt_norm);
1233    norm = _norm;
1234    norm2 = norm + M*eBands[m->nbEBands-1]-norm_offset;
1235    /* We can use the last band as scratch space because we don't need that
1236       scratch space for the last band */
1237    lowband_scratch = X_+M*eBands[m->nbEBands-1];
1238
1239    lowband_offset = 0;
1240    for (i=start;i<end;i++)
1241    {
1242       opus_int32 tell;
1243       int b;
1244       int N;
1245       opus_int32 curr_balance;
1246       int effective_lowband=-1;
1247       celt_norm * OPUS_RESTRICT X, * OPUS_RESTRICT Y;
1248       int tf_change=0;
1249       unsigned x_cm;
1250       unsigned y_cm;
1251       int last;
1252
1253       last = (i==end-1);
1254
1255       X = X_+M*eBands[i];
1256       if (Y_!=NULL)
1257          Y = Y_+M*eBands[i];
1258       else
1259          Y = NULL;
1260       N = M*eBands[i+1]-M*eBands[i];
1261       tell = ec_tell_frac(ec);
1262
1263       /* Compute how many bits we want to allocate to this band */
1264       if (i != start)
1265          balance -= tell;
1266       remaining_bits = total_bits-tell-1;
1267       if (i <= codedBands-1)
1268       {
1269          curr_balance = balance / IMIN(3, codedBands-i);
1270          b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1271       } else {
1272          b = 0;
1273       }
1274
1275       if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1276             lowband_offset = i;
1277
1278       tf_change = tf_res[i];
1279       if (i>=m->effEBands)
1280       {
1281          X=norm;
1282          if (Y_!=NULL)
1283             Y = norm;
1284          lowband_scratch = NULL;
1285       }
1286       if (i==end-1)
1287          lowband_scratch = NULL;
1288
1289       /* Get a conservative estimate of the collapse_mask's for the bands we're
1290           going to be folding from. */
1291       if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1292       {
1293          int fold_start;
1294          int fold_end;
1295          int fold_i;
1296          /* This ensures we never repeat spectral content within one band */
1297          effective_lowband = IMAX(0, M*eBands[lowband_offset]-norm_offset-N);
1298          fold_start = lowband_offset;
1299          while(M*eBands[--fold_start] > effective_lowband+norm_offset);
1300          fold_end = lowband_offset-1;
1301          while(M*eBands[++fold_end] < effective_lowband+norm_offset+N);
1302          x_cm = y_cm = 0;
1303          fold_i = fold_start; do {
1304            x_cm |= collapse_masks[fold_i*C+0];
1305            y_cm |= collapse_masks[fold_i*C+C-1];
1306          } while (++fold_i<fold_end);
1307       }
1308       /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1309           always) be non-zero.*/
1310       else
1311          x_cm = y_cm = (1<<B)-1;
1312
1313       if (dual_stereo && i==intensity)
1314       {
1315          int j;
1316
1317          /* Switch off dual stereo to do intensity */
1318          dual_stereo = 0;
1319          if (resynth)
1320             for (j=0;j<M*eBands[i]-norm_offset;j++)
1321                norm[j] = HALF32(norm[j]+norm2[j]);
1322       }
1323       if (dual_stereo)
1324       {
1325          x_cm = quant_band(encode, m, i, X, NULL, N, b/2, spread, B, intensity, tf_change,
1326                effective_lowband != -1 ? norm+effective_lowband : NULL, ec, &remaining_bits, LM,
1327                last?NULL:norm+M*eBands[i]-norm_offset, bandE, 0, seed, Q15ONE, lowband_scratch, x_cm);
1328          y_cm = quant_band(encode, m, i, Y, NULL, N, b/2, spread, B, intensity, tf_change,
1329                effective_lowband != -1 ? norm2+effective_lowband : NULL, ec, &remaining_bits, LM,
1330                last?NULL:norm2+M*eBands[i]-norm_offset, bandE, 0, seed, Q15ONE, lowband_scratch, y_cm);
1331       } else {
1332          x_cm = quant_band(encode, m, i, X, Y, N, b, spread, B, intensity, tf_change,
1333                effective_lowband != -1 ? norm+effective_lowband : NULL, ec, &remaining_bits, LM,
1334                last?NULL:norm+M*eBands[i]-norm_offset, bandE, 0, seed, Q15ONE, lowband_scratch, x_cm|y_cm);
1335          y_cm = x_cm;
1336       }
1337       collapse_masks[i*C+0] = (unsigned char)x_cm;
1338       collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1339       balance += pulses[i] + tell;
1340
1341       /* Update the folding position only as long as we have 1 bit/sample depth */
1342       update_lowband = b>(N<<BITRES);
1343    }
1344    RESTORE_STACK;
1345 }
1346