bands.c 27.4 KB
Newer Older
1
/* (C) 2007-2008 Jean-Marc Valin, CSIRO
2
   (C) 2008-2009 Gregory Maxwell */
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
/*
   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:
   
   - Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
   
   - Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
   
   - Neither the name of the Xiph.org Foundation nor the names of its
   contributors may be used to endorse or promote products derived from
   this software without specific prior written permission.
   
   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

32
33
34
35
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

36
37
#include <math.h>
#include "bands.h"
38
#include "modes.h"
39
#include "vq.h"
Jean-Marc Valin's avatar
Jean-Marc Valin committed
40
#include "cwrs.h"
41
#include "stack_alloc.h"
42
#include "os_support.h"
43
#include "mathops.h"
44
#include "rate.h"
45

46
47
const celt_word16_t sqrtC_1[2] = {QCONST16(1.f, 14), QCONST16(1.414214f, 14)};

48
49


50
#ifdef FIXED_POINT
51
52
/* Compute the amplitude (sqrt energy) in each of the bands */
void compute_band_energies(const CELTMode *m, const celt_sig_t *X, celt_ener_t *bank)
53
{
54
   int i, c, N;
55
   const celt_int16_t *eBands = m->eBands;
56
   const int C = CHANNELS(m);
57
   N = FRAMESIZE(m);
58
   for (c=0;c<C;c++)
59
   {
60
61
62
      for (i=0;i<m->nbEBands;i++)
      {
         int j;
63
64
         celt_word32_t maxval=0;
         celt_word32_t sum = 0;
65
66
         
         j=eBands[i]; do {
67
68
            maxval = MAX32(maxval, X[j+c*N]);
            maxval = MAX32(maxval, -X[j+c*N]);
69
70
         } while (++j<eBands[i+1]);
         
71
72
         if (maxval > 0)
         {
73
            int shift = celt_ilog2(maxval)-10;
74
            j=eBands[i]; do {
75
76
               sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
                                   EXTRACT16(VSHR32(X[j+c*N],shift)));
77
            } while (++j<eBands[i+1]);
78
79
            /* We're adding one here to make damn sure we never end up with a pitch vector that's
               larger than unity norm */
80
            bank[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
81
         } else {
82
            bank[i+c*m->nbEBands] = EPSILON;
83
         }
84
         /*printf ("%f ", bank[i+c*m->nbEBands]);*/
85
86
87
88
89
90
      }
   }
   /*printf ("\n");*/
}

/* Normalise each band such that the energy is one. */
91
void normalise_bands(const CELTMode *m, const celt_sig_t * restrict freq, celt_norm_t * restrict X, const celt_ener_t *bank)
92
{
93
   int i, c, N;
94
   const celt_int16_t *eBands = m->eBands;
95
   const int C = CHANNELS(m);
96
   N = FRAMESIZE(m);
97
98
   for (c=0;c<C;c++)
   {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
99
      i=0; do {
100
101
102
         celt_word16_t g;
         int j,shift;
         celt_word16_t E;
103
104
         shift = celt_zlog2(bank[i+c*m->nbEBands])-13;
         E = VSHR32(bank[i+c*m->nbEBands], shift);
105
         g = EXTRACT16(celt_rcp(SHL32(E,3)));
106
         j=eBands[i]; do {
107
            X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
108
         } while (++j<eBands[i+1]);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
109
      } while (++i<m->nbEBands);
110
111
112
   }
}

113
#else /* FIXED_POINT */
114
115
116
/* Compute the amplitude (sqrt energy) in each of the bands */
void compute_band_energies(const CELTMode *m, const celt_sig_t *X, celt_ener_t *bank)
{
117
   int i, c, N;
118
   const celt_int16_t *eBands = m->eBands;
119
   const int C = CHANNELS(m);
120
   N = FRAMESIZE(m);
121
122
123
124
125
   for (c=0;c<C;c++)
   {
      for (i=0;i<m->nbEBands;i++)
      {
         int j;
126
         celt_word32_t sum = 1e-10;
127
         for (j=eBands[i];j<eBands[i+1];j++)
128
            sum += X[j+c*N]*X[j+c*N];
129
130
         bank[i+c*m->nbEBands] = sqrt(sum);
         /*printf ("%f ", bank[i+c*m->nbEBands]);*/
131
132
133
134
135
      }
   }
   /*printf ("\n");*/
}

136
137
138
#ifdef EXP_PSY
void compute_noise_energies(const CELTMode *m, const celt_sig_t *X, const celt_word16_t *tonality, celt_ener_t *bank)
{
139
   int i, c, N;
140
141
   const celt_int16_t *eBands = m->eBands;
   const int C = CHANNELS(m);
142
   N = FRAMESIZE(m);
143
144
145
146
147
148
149
   for (c=0;c<C;c++)
   {
      for (i=0;i<m->nbEBands;i++)
      {
         int j;
         celt_word32_t sum = 1e-10;
         for (j=eBands[i];j<eBands[i+1];j++)
150
            sum += X[j*C+c]*X[j+c*N]*tonality[j];
151
152
         bank[i+c*m->nbEBands] = sqrt(sum);
         /*printf ("%f ", bank[i+c*m->nbEBands]);*/
153
154
155
156
157
158
      }
   }
   /*printf ("\n");*/
}
#endif

159
/* Normalise each band such that the energy is one. */
160
void normalise_bands(const CELTMode *m, const celt_sig_t * restrict freq, celt_norm_t * restrict X, const celt_ener_t *bank)
161
{
162
   int i, c, N;
163
   const celt_int16_t *eBands = m->eBands;
164
   const int C = CHANNELS(m);
165
   N = FRAMESIZE(m);
166
167
168
169
170
   for (c=0;c<C;c++)
   {
      for (i=0;i<m->nbEBands;i++)
      {
         int j;
171
         celt_word16_t g = 1.f/(1e-10+bank[i+c*m->nbEBands]);
172
         for (j=eBands[i];j<eBands[i+1];j++)
173
            X[j+c*N] = freq[j+c*N]*g;
174
175
176
177
      }
   }
}

178
179
#endif /* FIXED_POINT */

180
void renormalise_bands(const CELTMode *m, celt_norm_t * restrict X)
181
{
182
183
184
185
186
187
   int i, c;
   const celt_int16_t *eBands = m->eBands;
   const int C = CHANNELS(m);
   for (c=0;c<C;c++)
   {
      i=0; do {
188
         renormalise_vector(X+eBands[i]+c*m->eBands[m->nbEBands+1], Q15ONE, eBands[i+1]-eBands[i], 1);
189
190
      } while (++i<m->nbEBands);
   }
191
192
}

193
/* De-normalise the energy to produce the synthesis from the unit-energy bands */
194
void denormalise_bands(const CELTMode *m, const celt_norm_t * restrict X, celt_sig_t * restrict freq, const celt_ener_t *bank)
195
{
196
   int i, c, N;
197
   const celt_int16_t *eBands = m->eBands;
198
   const int C = CHANNELS(m);
199
   N = FRAMESIZE(m);
200
201
   if (C>2)
      celt_fatal("denormalise_bands() not implemented for >2 channels");
202
   for (c=0;c<C;c++)
203
   {
204
205
206
      for (i=0;i<m->nbEBands;i++)
      {
         int j;
207
         celt_word32_t g = SHR32(bank[i+c*m->nbEBands],1);
208
         j=eBands[i]; do {
209
            freq[j+c*N] = SHL32(MULT16_32_Q15(X[j+c*N], g),2);
210
         } while (++j<eBands[i+1]);
211
      }
212
213
      for (i=eBands[m->nbEBands];i<eBands[m->nbEBands+1];i++)
         freq[i+c*N] = 0;
214
215
216
   }
}

217
int compute_pitch_gain(const CELTMode *m, const celt_sig_t *X, const celt_sig_t *P, int norm_rate, int *gain_id)
218
{
219
   int j, c;
220
   celt_word16_t g;
221
   celt_word16_t delta;
222
   const int C = CHANNELS(m);
223
   celt_word32_t Sxy=0, Sxx=0, Syy=0;
224
   int len = m->pitchEnd;
225
   const int N = FRAMESIZE(m);
226
227
228
#ifdef FIXED_POINT
   int shift = 0;
   celt_word32_t maxabs=0;
229
230

   for (c=0;c<C;c++)
231
   {
232
233
234
235
236
      for (j=0;j<len;j++)
      {
         maxabs = MAX32(maxabs, ABS32(X[j+c*N]));
         maxabs = MAX32(maxabs, ABS32(P[j+c*N]));
      }
237
238
239
240
   }
   shift = celt_ilog2(maxabs)-12;
   if (shift<0)
      shift = 0;
241
#endif
242
243
   delta = PDIV32_16(Q15ONE, len);
   for (c=0;c<C;c++)
244
   {
245
246
247
248
      celt_word16_t gg = Q15ONE;
      for (j=0;j<len;j++)
      {
         celt_word16_t Xj, Pj;
249
250
         Xj = EXTRACT16(SHR32(X[j+c*N], shift));
         Pj = MULT16_16_P15(gg,EXTRACT16(SHR32(P[j+c*N], shift)));
251
252
253
254
255
         Sxy = MAC16_16(Sxy, Xj, Pj);
         Sxx = MAC16_16(Sxx, Pj, Pj);
         Syy = MAC16_16(Syy, Xj, Xj);
         gg = SUB16(gg, delta);
      }
256
   }
257
#ifdef FIXED_POINT
258
259
   {
      celt_word32_t num, den;
260
261
262
263
      celt_word16_t fact;
      fact = MULT16_16(QCONST16(.04, 14), norm_rate);
      if (fact < QCONST16(1., 14))
         fact = QCONST16(1., 14);
264
265
266
267
268
269
      num = Sxy;
      den = EPSILON+Sxx+MULT16_32_Q15(QCONST16(.03,15),Syy);
      shift = celt_ilog2(Sxy)-16;
      if (shift < 0)
         shift = 0;
      g = DIV32(SHL32(SHR32(num,shift),14),SHR32(den,shift));
270
      if (Sxy < MULT16_32_Q15(fact, MULT16_16(celt_sqrt(EPSILON+Sxx),celt_sqrt(EPSILON+Syy))))
271
         g = 0;
272
      /* This MUST round down so that we don't over-estimate the gain */
273
274
      *gain_id = EXTRACT16(SHR32(MULT16_16(20,(g-QCONST16(.5,14))),14));
   }
275
#else
276
277
278
279
280
281
282
   {
      float fact = .04*norm_rate;
      if (fact < 1)
         fact = 1;
      g = Sxy/(.1+Sxx+.03*Syy);
      if (Sxy < .5*fact*celt_sqrt(1+Sxx*Syy))
         g = 0;
283
      /* This MUST round down so that we don't over-estimate the gain */
284
285
      *gain_id = floor(20*(g-.5));
   }
286
#endif
287
288
289
290
291
292
293
294
295
296
297
   if (*gain_id < 0)
   {
      *gain_id = 0;
      return 0;
   } else {
      if (*gain_id > 15)
         *gain_id = 15;
      return 1;
   }
}

298
void apply_pitch(const CELTMode *m, celt_sig_t *X, const celt_sig_t *P, int gain_id, int pred)
299
{
300
   int j, c, N;
301
   celt_word16_t gain;
302
   celt_word16_t delta;
303
   const int C = CHANNELS(m);
304
   int len = m->pitchEnd;
305
306

   N = FRAMESIZE(m);
307
   gain = ADD16(QCONST16(.5,14), MULT16_16_16(QCONST16(.05,14),gain_id));
308
   delta = PDIV32_16(gain, len);
309
310
   if (pred)
      gain = -gain;
311
312
313
   else
      delta = -delta;
   for (c=0;c<C;c++)
314
   {
315
316
317
      celt_word16_t gg = gain;
      for (j=0;j<len;j++)
      {
318
         X[j+c*N] += SHL32(MULT16_32_Q15(gg,P[j+c*N]),1);
319
320
         gg = ADD16(gg, delta);
      }
321
322
   }
}
323

324
325
#ifndef DISABLE_STEREO

326
static void stereo_band_mix(const CELTMode *m, celt_norm_t *X, celt_norm_t *Y, const celt_ener_t *bank, int stereo_mode, int bandID, int dir)
327
328
329
{
   int i = bandID;
   const celt_int16_t *eBands = m->eBands;
330
331
332
   int j;
   celt_word16_t a1, a2;
   if (stereo_mode==0)
333
   {
334
335
336
337
338
339
      /* Do mid-side when not doing intensity stereo */
      a1 = QCONST16(.70711f,14);
      a2 = dir*QCONST16(.70711f,14);
   } else {
      celt_word16_t left, right;
      celt_word16_t norm;
340
#ifdef FIXED_POINT
341
      int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
342
#endif
343
344
345
346
347
      left = VSHR32(bank[i],shift);
      right = VSHR32(bank[i+m->nbEBands],shift);
      norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
      a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
      a2 = dir*DIV32_16(SHL32(EXTEND32(right),14),norm);
348
   }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
349
   for (j=0;j<eBands[i+1]-eBands[i];j++)
350
351
   {
      celt_norm_t r, l;
352
353
354
355
      l = X[j];
      r = Y[j];
      X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
      Y[j] = MULT16_16_Q14(a1,r) - MULT16_16_Q14(a2,l);
356
357
358
   }
}

359

360
361
#endif /* DISABLE_STEREO */

362
363
int folding_decision(const CELTMode *m, celt_norm_t *X, celt_word16_t *average, int *last_decision)
{
364
   int i, c, N0;
365
366
   int NR=0;
   celt_word32_t ratio = EPSILON;
367
   const int C = CHANNELS(m);
368
   const celt_int16_t * restrict eBands = m->eBands;
369
370
371
   
   N0 = FRAMESIZE(m);

372
373
   for (c=0;c<C;c++)
   {
374
375
376
377
378
379
   for (i=0;i<m->nbEBands;i++)
   {
      int j, N;
      int max_i=0;
      celt_word16_t max_val=EPSILON;
      celt_word32_t floor_ener=EPSILON;
380
      celt_norm_t * restrict x = X+eBands[i]+c*N0;
381
382
383
      N = eBands[i+1]-eBands[i];
      for (j=0;j<N;j++)
      {
384
         if (ABS16(x[j])>max_val)
385
         {
386
            max_val = ABS16(x[j]);
387
388
389
390
391
392
393
394
395
396
397
398
            max_i = j;
         }
      }
#if 0
      for (j=0;j<N;j++)
      {
         if (abs(j-max_i)>2)
            floor_ener += x[j]*x[j];
      }
#else
      floor_ener = QCONST32(1.,28)-MULT16_16(max_val,max_val);
      if (max_i < N-1)
399
         floor_ener -= MULT16_16(x[(max_i+1)], x[(max_i+1)]);
400
      if (max_i < N-2)
401
         floor_ener -= MULT16_16(x[(max_i+2)], x[(max_i+2)]);
402
      if (max_i > 0)
403
         floor_ener -= MULT16_16(x[(max_i-1)], x[(max_i-1)]);
404
      if (max_i > 1)
405
         floor_ener -= MULT16_16(x[(max_i-2)], x[(max_i-2)]);
406
407
      floor_ener = MAX32(floor_ener, EPSILON);
#endif
408
      if (N>7)
409
410
411
412
413
414
415
416
417
      {
         celt_word16_t r;
         celt_word16_t den = celt_sqrt(floor_ener);
         den = MAX32(QCONST16(.02, 15), den);
         r = DIV32_16(SHL32(EXTEND32(max_val),8),den);
         ratio = ADD32(ratio, EXTEND32(r));
         NR++;
      }
   }
418
   }
419
420
421
422
423
424
425
426
427
428
429
430
431
   if (NR>0)
      ratio = DIV32_16(ratio, NR);
   ratio = ADD32(HALF32(ratio), HALF32(*average));
   if (!*last_decision)
   {
      *last_decision = (ratio < QCONST16(1.8,8));
   } else {
      *last_decision = (ratio < QCONST16(3.,8));
   }
   *average = EXTRACT16(ratio);
   return *last_decision;
}

432
/* Quantisation of the residual */
433
void quant_bands(const CELTMode *m, celt_norm_t * restrict X, const celt_ener_t *bandE, int *pulses, int shortBlocks, int fold, int total_bits, ec_enc *enc)
434
{
435
   int i, j, remaining_bits, balance;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
436
   const celt_int16_t * restrict eBands = m->eBands;
437
   celt_norm_t * restrict norm;
438
   VARDECL(celt_norm_t, _norm);
439
   int B;
440
441
   SAVE_STACK;

442
   B = shortBlocks ? m->nbShortMdcts : 1;
443
   ALLOC(_norm, eBands[m->nbEBands+1], celt_norm_t);
444
   norm = _norm;
445

446
   balance = 0;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
447
   for (i=0;i<m->nbEBands;i++)
448
   {
449
      int tell;
450
      int N;
451
      int q;
452
      celt_word16_t n;
453
      const celt_int16_t * const *BPbits;
454
455
456
      
      int curr_balance, curr_bits;
      
457
      N = eBands[i+1]-eBands[i];
458
      BPbits = m->bits;
459

460
      tell = ec_enc_tell(enc, BITRES);
461
462
      if (i != 0)
         balance -= tell;
463
      remaining_bits = (total_bits<<BITRES)-tell-1;
464
465
466
467
      curr_balance = (m->nbEBands-i);
      if (curr_balance > 3)
         curr_balance = 3;
      curr_balance = balance / curr_balance;
468
469
      q = bits2pulses(m, BPbits[i], N, pulses[i]+curr_balance);
      curr_bits = pulses2bits(BPbits[i], N, q);
470
      remaining_bits -= curr_bits;
471
      while (remaining_bits < 0 && q > 0)
472
473
      {
         remaining_bits += curr_bits;
474
         q--;
475
         curr_bits = pulses2bits(BPbits[i], N, q);
476
         remaining_bits -= curr_bits;
477
478
479
      }
      balance += pulses[i] + tell;
      
480
      n = SHL16(celt_sqrt(eBands[i+1]-eBands[i]),11);
481

482
      if (q > 0)
483
      {
484
         int spread = fold ? B : 0;
485
         alg_quant(X+eBands[i], eBands[i+1]-eBands[i], q, spread, enc);
486
      } else {
487
         intra_fold(m, eBands[i+1]-eBands[i], norm, X+eBands[i], eBands[i], B);
488
      }
489
      for (j=eBands[i];j<eBands[i+1];j++)
490
         norm[j] = MULT16_16_Q15(n,X[j]);
491
   }
492
   RESTORE_STACK;
493
494
}

495
496
#ifndef DISABLE_STEREO

497
void quant_bands_stereo(const CELTMode *m, celt_norm_t *_X, const celt_ener_t *bandE, int *pulses, int shortBlocks, int fold, int total_bits, ec_enc *enc)
498
499
500
{
   int i, j, remaining_bits, balance;
   const celt_int16_t * restrict eBands = m->eBands;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
501
   celt_norm_t * restrict norm;
502
503
504
505
506
507
   VARDECL(celt_norm_t, _norm);
   int B;
   celt_word16_t mid, side;
   SAVE_STACK;

   B = shortBlocks ? m->nbShortMdcts : 1;
508
   ALLOC(_norm, eBands[m->nbEBands+1], celt_norm_t);
509
510
511
512
513
   norm = _norm;

   balance = 0;
   for (i=0;i<m->nbEBands;i++)
   {
514
      int c;
515
516
517
518
519
520
521
522
523
524
      int tell;
      int q1, q2;
      celt_word16_t n;
      const celt_int16_t * const *BPbits;
      int b, qb;
      int N;
      int curr_balance, curr_bits;
      int imid, iside, itheta;
      int mbits, sbits, delta;
      int qalloc;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
525
      celt_norm_t * restrict X, * restrict Y;
526
      
Jean-Marc Valin's avatar
Jean-Marc Valin committed
527
528
      X = _X+eBands[i];
      Y = X+eBands[m->nbEBands+1];
529
530
531
      BPbits = m->bits;

      N = eBands[i+1]-eBands[i];
532
      tell = ec_enc_tell(enc, BITRES);
533
534
535
536
537
538
539
      if (i != 0)
         balance -= tell;
      remaining_bits = (total_bits<<BITRES)-tell-1;
      curr_balance = (m->nbEBands-i);
      if (curr_balance > 3)
         curr_balance = 3;
      curr_balance = balance / curr_balance;
540
      b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
541
542
543
      if (b<0)
         b = 0;

544
      qb = (b-2*(N-1)*(QTHETA_OFFSET-log2_frac(N,BITRES)))/(32*(N-1));
545
546
547
548
      if (qb > (b>>BITRES)-1)
         qb = (b>>BITRES)-1;
      if (qb<0)
         qb = 0;
549
550
      if (qb>14)
         qb = 14;
551

552
      stereo_band_mix(m, X, Y, bandE, qb==0, i, 1);
553

Jean-Marc Valin's avatar
Jean-Marc Valin committed
554
555
      mid = renormalise_vector(X, Q15ONE, N, 1);
      side = renormalise_vector(Y, Q15ONE, N, 1);
556
557
558
559
560
#ifdef FIXED_POINT
      itheta = MULT16_16_Q15(QCONST16(0.63662,15),celt_atan2p(side, mid));
#else
      itheta = floor(.5+16384*0.63662*atan2(side,mid));
#endif
561
      qalloc = log2_frac((1<<qb)+1,BITRES);
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
      if (qb==0)
      {
         itheta=0;
      } else {
         int shift;
         shift = 14-qb;
         itheta = (itheta+(1<<shift>>1))>>shift;
         ec_enc_uint(enc, itheta, (1<<qb)+1);
         itheta <<= shift;
      }
      if (itheta == 0)
      {
         imid = 32767;
         iside = 0;
         delta = -10000;
      } else if (itheta == 16384)
      {
         imid = 0;
         iside = 32767;
         delta = 10000;
      } else {
         imid = bitexact_cos(itheta);
         iside = bitexact_cos(16384-itheta);
585
         delta = (N-1)*(log2_frac(iside,BITRES+2)-log2_frac(imid,BITRES+2))>>2;
586
      }
587
588
589
590
591
592
593
      n = SHL16(celt_sqrt((eBands[i+1]-eBands[i])),11);

      if (N==2)
      {
         int c2;
         int sign=1;
         celt_norm_t v[2], w[2];
594
         celt_norm_t *x2, *y2;
595
596
597
598
599
600
601
602
         mbits = b-qalloc;
         sbits = 0;
         if (itheta != 0 && itheta != 16384)
            sbits = 1<<BITRES;
         mbits -= sbits;
         c = itheta > 8192 ? 1 : 0;
         c2 = 1-c;

Jean-Marc Valin's avatar
Jean-Marc Valin committed
603
604
         x2 = X;
         y2 = Y;
605
606
607
608
609
610
611
612
613
614
615
616
         if (c==0)
         {
            v[0] = x2[0];
            v[1] = x2[1];
            w[0] = y2[0];
            w[1] = y2[1];
         } else {
            v[0] = y2[0];
            v[1] = y2[1];
            w[0] = x2[0];
            w[1] = x2[1];
         }
617
618
619
620
621
622
623
624
625
626
627
628
         q1 = bits2pulses(m, BPbits[i], N, mbits);
         curr_bits = pulses2bits(BPbits[i], N, q1)+qalloc+sbits;
         remaining_bits -= curr_bits;
         while (remaining_bits < 0 && q1 > 0)
         {
            remaining_bits += curr_bits;
            q1--;
            curr_bits = pulses2bits(BPbits[i], N, q1)+qalloc;
            remaining_bits -= curr_bits;
         }

         if (q1 > 0)
629
         {
630
            int spread = fold ? B : 0;
631
            alg_quant(v, N, q1, spread, enc);
632
         } else {
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
            v[0] = QCONST16(1.f, 14);
            v[1] = 0;
         }
         if (sbits)
         {
            if (v[0]*w[1] - v[1]*w[0] > 0)
               sign = 1;
            else
               sign = -1;
            ec_enc_bits(enc, sign==1, 1);
         } else {
            sign = 1;
         }
         w[0] = -sign*v[1];
         w[1] = sign*v[0];
         if (c==0)
         {
            x2[0] = v[0];
            x2[1] = v[1];
652
653
            y2[0] = w[0];
            y2[1] = w[1];
654
655
656
         } else {
            x2[0] = w[0];
            x2[1] = w[1];
657
658
            y2[0] = v[0];
            y2[1] = v[1];
659
660
661
         }
      } else {
         
662
663
664
665
666
667
668
669
670
671
672
         mbits = (b-qalloc/2-delta)/2;
         if (mbits > b-qalloc)
            mbits = b-qalloc;
         if (mbits<0)
            mbits=0;
         sbits = b-qalloc-mbits;
         q1 = bits2pulses(m, BPbits[i], N, mbits);
         q2 = bits2pulses(m, BPbits[i], N, sbits);
         curr_bits = pulses2bits(BPbits[i], N, q1)+pulses2bits(BPbits[i], N, q2)+qalloc;
         remaining_bits -= curr_bits;
         while (remaining_bits < 0 && (q1 > 0 || q2 > 0))
673
         {
674
675
676
677
678
679
680
681
682
683
            remaining_bits += curr_bits;
            if (q1>q2)
            {
               q1--;
               curr_bits = pulses2bits(BPbits[i], N, q1)+pulses2bits(BPbits[i], N, q2)+qalloc;
            } else {
               q2--;
               curr_bits = pulses2bits(BPbits[i], N, q1)+pulses2bits(BPbits[i], N, q2)+qalloc;
            }
            remaining_bits -= curr_bits;
684
685
         }

686
687
         if (q1 > 0) {
            int spread = fold ? B : 0;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
688
            alg_quant(X, N, q1, spread, enc);
689
         } else {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
690
            intra_fold(m, eBands[i+1]-eBands[i], norm, X, eBands[i], B);
691
692
693
         }
         if (q2 > 0) {
            int spread = fold ? B : 0;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
694
            alg_quant(Y, N, q2, spread, enc);
695
         } else
Jean-Marc Valin's avatar
Jean-Marc Valin committed
696
            for (j=0;j<N;j++)
697
               Y[j] = 0;
698
699
700
      }
      
      balance += pulses[i] + tell;
701
702
703
704
705
706
707
708

#ifdef FIXED_POINT
      mid = imid;
      side = iside;
#else
      mid = (1./32768)*imid;
      side = (1./32768)*iside;
#endif
709
      for (j=0;j<N;j++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
710
         norm[eBands[i]+j] = MULT16_16_Q15(n,X[j]);
711

712
      for (j=0;j<N;j++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
713
         X[j] = MULT16_16_Q15(X[j], mid);
714
      for (j=0;j<N;j++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
715
         Y[j] = MULT16_16_Q15(Y[j], side);
716

717
      stereo_band_mix(m, X, Y, bandE, 0, i, -1);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
718
719
      renormalise_vector(X, Q15ONE, N, 1);
      renormalise_vector(Y, Q15ONE, N, 1);
720
721
722
   }
   RESTORE_STACK;
}
723
#endif /* DISABLE_STEREO */
724

725
/* Decoding of the residual */
726
void unquant_bands(const CELTMode *m, celt_norm_t * restrict X, const celt_ener_t *bandE, int *pulses, int shortBlocks, int fold, int total_bits, ec_dec *dec)
727
{
728
   int i, j, remaining_bits, balance;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
729
   const celt_int16_t * restrict eBands = m->eBands;
730
731
   celt_norm_t * restrict norm;
   VARDECL(celt_norm_t, _norm);
732
   int B;
733
734
   SAVE_STACK;

735
   B = shortBlocks ? m->nbShortMdcts : 1;
736
   ALLOC(_norm, eBands[m->nbEBands+1], celt_norm_t);
737
   norm = _norm;
738

739
   balance = 0;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
740
   for (i=0;i<m->nbEBands;i++)
741
   {
742
      int tell;
743
      int N;
744
      int q;
745
      celt_word16_t n;
746
      const celt_int16_t * const *BPbits;
747
748
      
      int curr_balance, curr_bits;
749
750

      N = eBands[i+1]-eBands[i];
751
      BPbits = m->bits;
752

753
      tell = ec_dec_tell(dec, BITRES);
754
755
      if (i != 0)
         balance -= tell;
756
      remaining_bits = (total_bits<<BITRES)-tell-1;
757
758
759
760
      curr_balance = (m->nbEBands-i);
      if (curr_balance > 3)
         curr_balance = 3;
      curr_balance = balance / curr_balance;
761
762
      q = bits2pulses(m, BPbits[i], N, pulses[i]+curr_balance);
      curr_bits = pulses2bits(BPbits[i], N, q);
763
      remaining_bits -= curr_bits;
764
      while (remaining_bits < 0 && q > 0)
765
766
      {
         remaining_bits += curr_bits;
767
         q--;
768
         curr_bits = pulses2bits(BPbits[i], N, q);
769
         remaining_bits -= curr_bits;
770
771
772
      }
      balance += pulses[i] + tell;

773
      n = SHL16(celt_sqrt(eBands[i+1]-eBands[i]),11);
774

775
      if (q > 0)
776
      {
777
         int spread = fold ? B : 0;
778
         alg_unquant(X+eBands[i], eBands[i+1]-eBands[i], q, spread, dec);
779
      } else {
780
         intra_fold(m, eBands[i+1]-eBands[i], norm, X+eBands[i], eBands[i], B);
781
      }
782
      for (j=eBands[i];j<eBands[i+1];j++)
783
         norm[j] = MULT16_16_Q15(n,X[j]);
784
   }
785
   RESTORE_STACK;
786
}
787

788
789
#ifndef DISABLE_STEREO

790
void unquant_bands_stereo(const CELTMode *m, celt_norm_t *_X, const celt_ener_t *bandE, int *pulses, int shortBlocks, int fold, int total_bits, ec_dec *dec)
791
792
793
{
   int i, j, remaining_bits, balance;
   const celt_int16_t * restrict eBands = m->eBands;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
794
   celt_norm_t * restrict norm;
795
796
797
798
799
800
   VARDECL(celt_norm_t, _norm);
   int B;
   celt_word16_t mid, side;
   SAVE_STACK;

   B = shortBlocks ? m->nbShortMdcts : 1;
801
   ALLOC(_norm, eBands[m->nbEBands+1], celt_norm_t);
802
803
804
805
806
   norm = _norm;

   balance = 0;
   for (i=0;i<m->nbEBands;i++)
   {
807
      int c;
808
809
810
811
812
813
814
815
816
817
      int tell;
      int q1, q2;
      celt_word16_t n;
      const celt_int16_t * const *BPbits;
      int b, qb;
      int N;
      int curr_balance, curr_bits;
      int imid, iside, itheta;
      int mbits, sbits, delta;
      int qalloc;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
818
      celt_norm_t * restrict X, * restrict Y;
819
      
Jean-Marc Valin's avatar
Jean-Marc Valin committed
820
821
      X = _X+eBands[i];
      Y = X+eBands[m->nbEBands+1];
822
823
824
      BPbits = m->bits;

      N = eBands[i+1]-eBands[i];
825
      tell = ec_dec_tell(dec, BITRES);
826
827
828
829
830
831
832
      if (i != 0)
         balance -= tell;
      remaining_bits = (total_bits<<BITRES)-tell-1;
      curr_balance = (m->nbEBands-i);
      if (curr_balance > 3)
         curr_balance = 3;
      curr_balance = balance / curr_balance;
833
      b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
834
835
836
      if (b<0)
         b = 0;

837
      qb = (b-2*(N-1)*(QTHETA_OFFSET-log2_frac(N,BITRES)))/(32*(N-1));
838
839
      if (qb > (b>>BITRES)-1)
         qb = (b>>BITRES)-1;
840
841
      if (qb>14)
         qb = 14;
842
843
      if (qb<0)
         qb = 0;
844
      qalloc = log2_frac((1<<qb)+1,BITRES);
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
      if (qb==0)
      {
         itheta=0;
      } else {
         int shift;
         shift = 14-qb;
         itheta = ec_dec_uint(dec, (1<<qb)+1);
         itheta <<= shift;
      }
      if (itheta == 0)
      {
         imid = 32767;
         iside = 0;
         delta = -10000;
      } else if (itheta == 16384)
      {
         imid = 0;
         iside = 32767;
         delta = 10000;
      } else {
         imid = bitexact_cos(itheta);
         iside = bitexact_cos(16384-itheta);
867
         delta = (N-1)*(log2_frac(iside,BITRES+2)-log2_frac(imid,BITRES+2))>>2;
868
      }
869
870
871
872
873
874
875
      n = SHL16(celt_sqrt((eBands[i+1]-eBands[i])),11);

      if (N==2)
      {
         int c2;
         int sign=1;
         celt_norm_t v[2], w[2];
876
         celt_norm_t *x2, *y2;
877
878
879
880
881
882
883
884
         mbits = b-qalloc;
         sbits = 0;
         if (itheta != 0 && itheta != 16384)
            sbits = 1<<BITRES;
         mbits -= sbits;
         c = itheta > 8192 ? 1 : 0;
         c2 = 1-c;

Jean-Marc Valin's avatar
Jean-Marc Valin committed
885
886
         x2 = X;
         y2 = Y;
887
         v[0] = x2[c];
888
         v[1] = y2[c];
889
         w[0] = x2[c2];
890
         w[1] = y2[c2];
891
892
893
894
895
896
897
898
899
900
901
902
         q1 = bits2pulses(m, BPbits[i], N, mbits);
         curr_bits