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