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