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