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