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