rate.c 20.5 KB
Newer Older
Jean-Marc Valin's avatar
Jean-Marc Valin committed
1
2
3
/* Copyright (c) 2007-2008 CSIRO
   Copyright (c) 2007-2009 Xiph.Org Foundation
   Written by Jean-Marc Valin */
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
32
/*
   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.
*/

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

37
38
39
40
#include <math.h>
#include "modes.h"
#include "cwrs.h"
#include "arch.h"
Jean-Marc Valin's avatar
Jean-Marc Valin committed
41
#include "os_support.h"
42
43

#include "entcode.h"
44
#include "rate.h"
45

Jean-Marc Valin's avatar
Jean-Marc Valin committed
46

47
48
49
50
51
52
53
54
static const unsigned char LOG2_FRAC_TABLE[24]={
   0,
   8,13,
  16,19,21,23,
  24,26,27,28,29,30,31,32,
  32,33,34,34,35,36,36,37,37
};

55
#ifdef CUSTOM_MODES
56

Jean-Marc Valin's avatar
Jean-Marc Valin committed
57
58
59
/*Determines if V(N,K) fits in a 32-bit unsigned integer.
  N and K are themselves limited to 15 bits.*/
static int fits_in32(int _n, int _k)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
60
{
Jean-Marc Valin's avatar
Jean-Marc Valin committed
61
62
63
64
65
66
67
   static const celt_int16 maxN[15] = {
      32767, 32767, 32767, 1476, 283, 109,  60,  40,
       29,  24,  20,  18,  16,  14,  13};
   static const celt_int16 maxK[15] = {
      32767, 32767, 32767, 32767, 1172, 238,  95,  53,
       36,  27,  22,  18,  16,  15,  13};
   if (_n>=14)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
68
   {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
69
70
      if (_k>=14)
         return 0;
71
      else
Jean-Marc Valin's avatar
Jean-Marc Valin committed
72
73
74
         return _n <= maxN[_k];
   } else {
      return _k <= maxK[_n];
Jean-Marc Valin's avatar
Jean-Marc Valin committed
75
   }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
76
77
78
79
}

void compute_pulse_cache(CELTMode *m, int LM)
{
80
   int C;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
81
   int i;
82
   int j;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
83
84
85
86
87
   int curr=0;
   int nbEntries=0;
   int entryN[100], entryK[100], entryI[100];
   const celt_int16 *eBands = m->eBands;
   PulseCache *cache = &m->cache;
88
89
   celt_int16 *cindex;
   unsigned char *bits;
90
   unsigned char *cap;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
91

92
93
   cindex = celt_alloc(sizeof(cache->index[0])*m->nbEBands*(LM+2));
   cache->index = cindex;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
94

95
   /* Scan for all unique band sizes */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
96
   for (i=0;i<=LM+1;i++)
97
   {
98
      for (j=0;j<m->nbEBands;j++)
99
      {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
100
101
         int k;
         int N = (eBands[j+1]-eBands[j])<<i>>1;
102
         cindex[i*m->nbEBands+j] = -1;
103
         /* Find other bands that have the same size */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
104
         for (k=0;k<=i;k++)
105
         {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
106
            int n;
107
            for (n=0;n<m->nbEBands && (k!=i || n<j);n++)
108
            {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
109
110
               if (N == (eBands[n+1]-eBands[n])<<k>>1)
               {
111
                  cindex[i*m->nbEBands+j] = cindex[k*m->nbEBands+n];
Jean-Marc Valin's avatar
Jean-Marc Valin committed
112
113
                  break;
               }
114
115
            }
         }
116
         if (cache->index[i*m->nbEBands+j] == -1 && N!=0)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
117
118
119
120
         {
            int K;
            entryN[nbEntries] = N;
            K = 0;
121
            while (fits_in32(N,get_pulses(K+1)) && K<MAX_PSEUDO)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
122
123
               K++;
            entryK[nbEntries] = K;
124
            cindex[i*m->nbEBands+j] = curr;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
125
126
127
128
129
130
131
            entryI[nbEntries] = curr;

            curr += K+1;
            nbEntries++;
         }
      }
   }
132
133
134
   bits = celt_alloc(sizeof(unsigned char)*curr);
   cache->bits = bits;
   cache->size = curr;
135
   /* Compute the cache for all unique sizes */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
136
137
   for (i=0;i<nbEntries;i++)
   {
138
      unsigned char *ptr = bits+entryI[i];
139
      celt_int16 tmp[MAX_PULSES+1];
Jean-Marc Valin's avatar
Jean-Marc Valin committed
140
141
142
143
      get_required_bits(tmp, entryN[i], get_pulses(entryK[i]), BITRES);
      for (j=1;j<=entryK[i];j++)
         ptr[j] = tmp[get_pulses(j)]-1;
      ptr[0] = entryK[i];
144
   }
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159

   /* Compute the maximum rate for each band at which we'll reliably use as
       many bits as we ask for. */
   cache->caps = cap = celt_alloc(sizeof(cache->caps[0])*(LM+1)*2*m->nbEBands);
   for (i=0;i<=LM;i++)
   {
      for (C=1;C<=2;C++)
      {
         for (j=0;j<m->nbEBands;j++)
         {
            int N0;
            int max_bits;
            N0 = m->eBands[j+1]-m->eBands[j];
            /* N=1 bands only have a sign bit and fine bits. */
            if (N0<<i == 1)
160
               max_bits = C*(1+MAX_FINE_BITS)<<BITRES;
161
162
163
164
165
166
167
168
169
170
171
172
            else
            {
               const unsigned char *pcache;
               celt_int32           num;
               celt_int32           den;
               int                  LM0;
               int                  N;
               int                  offset;
               int                  ndof;
               int                  qb;
               int                  k;
               LM0 = 0;
173
174
175
               /* Even-sized bands bigger than N=2 can be split one more
                   time. */
               if (N0 > 2 && !(N0&1))
176
177
178
179
               {
                  N0>>=1;
                  LM0--;
               }
180
181
               /* N0=1 bands can't be split down to N<2. */
               else if (N0 <= 1)
182
               {
183
                  LM0=IMIN(i,1);
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
                  N0<<=LM0;
               }
               /* Compute the cost for the lowest-level PVQ of a fully split
                   band. */
               pcache = bits + cindex[(LM0+1)*m->nbEBands+j];
               max_bits = pcache[pcache[0]]+1;
               /* Add in the cost of coding regular splits. */
               N = N0;
               for(k=0;k<i-LM0;k++){
                  max_bits <<= 1;
                  /* Offset the number of qtheta bits by log2(N)/2
                      + QTHETA_OFFSET compared to their "fair share" of
                      total/N */
                  offset = (m->logN[j]+(LM0+k<<BITRES)>>1)-QTHETA_OFFSET;
                  /* The number of qtheta bits we'll allocate if the remainder
199
200
201
202
203
204
                      is to be max_bits.
                     The average measured cost for theta is 0.89701 times qb,
                      approximated here as 459/512. */
                  num=459*(celt_int32)((2*N-1)*offset+max_bits);
                  den=((celt_int32)(2*N-1)<<9)-459;
                  qb = IMIN((num+(den>>1))/den, 57);
205
                  celt_assert(qb >= 0);
206
                  max_bits += qb;
207
208
209
210
211
212
                  N <<= 1;
               }
               /* Add in the cost of a stereo split, if necessary. */
               if (C==2)
               {
                  max_bits <<= 1;
213
                  offset = (m->logN[j]+(i<<BITRES)>>1)-(N==2?QTHETA_OFFSET_TWOPHASE:QTHETA_OFFSET);
214
                  ndof = 2*N-1-(N==2);
215
216
217
218
219
                  /* The average measured cost for theta with the step PDF is
                      0.95164 times qb, approximated here as 487/512. */
                  num = (N==2?512:487)*(celt_int32)(max_bits+ndof*offset);
                  den = ((celt_int32)ndof<<9)-(N==2?512:487);
                  qb = IMIN((num+(den>>1))/den, (N==2?64:61));
220
                  celt_assert(qb >= 0);
221
                  max_bits += qb;
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
               }
               /* Add the fine bits we'll use. */
               /* Compensate for the extra DoF in stereo */
               ndof = C*N + ((C==2 && N>2) ? 1 : 0);
               /* Offset the number of fine bits by log2(N)/2 + FINE_OFFSET
                   compared to their "fair share" of total/N */
               offset = (m->logN[j] + (i<<BITRES)>>1)-FINE_OFFSET;
               /* N=2 is the only point that doesn't match the curve */
               if (N==2)
                  offset += 1<<BITRES>>2;
               /* The number of fine bits we'll allocate if the remainder is
                   to be max_bits. */
               num = max_bits+ndof*offset;
               den = ndof-1<<BITRES;
               qb = IMIN((num+(den>>1))/den, MAX_FINE_BITS);
               celt_assert(qb >= 0);
               max_bits += C*qb<<BITRES;
            }
240
241
242
243
            max_bits = (4*max_bits/(C*(m->eBands[j+1]-m->eBands[j]<<i)))-64;
            celt_assert(max_bits >= 0);
            celt_assert(max_bits < 256);
            *cap++ = (unsigned char)max_bits;
244
245
246
         }
      }
   }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
247
248
}

249
#endif /* CUSTOM_MODES */
250

251

252
#define ALLOC_STEPS 6
253

254
static inline int interp_bits2pulses(const CELTMode *m, int start, int end, int skip_start,
Timothy B. Terriberry's avatar
Timothy B. Terriberry committed
255
      const int *bits1, const int *bits2, const int *thresh, const int *cap, celt_int32 total, celt_int32 *_balance,
256
      int skip_rsv, int *intensity, int intensity_rsv, int *dual_stereo, int dual_stereo_rsv, int *bits,
257
      int *ebits, int *fine_priority, int _C, int LM, ec_ctx *ec, int encode, int prev)
258
{
Timothy B. Terriberry's avatar
Timothy B. Terriberry committed
259
   celt_int32 psum;
260
   int lo, hi;
261
   int i, j;
262
   int logM;
263
   const int C = CHANNELS(_C);
264
   int stereo;
265
   int codedBands=-1;
266
   int alloc_floor;
Timothy B. Terriberry's avatar
Timothy B. Terriberry committed
267
   celt_int32 left, percoeff;
268
   int done;
269
   int balance;
270
   SAVE_STACK;
271

272
   alloc_floor = C<<BITRES;
273
   stereo = C>1;
274

275
   logM = LM<<BITRES;
276
   lo = 0;
277
278
   hi = 1<<ALLOC_STEPS;
   for (i=0;i<ALLOC_STEPS;i++)
279
280
   {
      int mid = (lo+hi)>>1;
281
      psum = 0;
282
      done = 0;
283
      for (j=end;j-->start;)
284
      {
Timothy B. Terriberry's avatar
Timothy B. Terriberry committed
285
         int tmp = bits1[j] + (mid*(celt_int32)bits2[j]>>ALLOC_STEPS);
286
         if (tmp >= thresh[j] || done)
287
288
         {
            done = 1;
289
            /* Don't allocate more than we can actually use */
290
            psum += IMIN(tmp, cap[j]);
291
         } else {
292
293
294
            if (tmp >= alloc_floor)
               psum += alloc_floor;
         }
295
      }
296
      if (psum > total)
297
298
299
         hi = mid;
      else
         lo = mid;
300
   }
301
   psum = 0;
302
   /*printf ("interp bisection gave %d\n", lo);*/
303
   done = 0;
304
   for (j=end;j-->start;)
305
   {
306
      int tmp = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
307
      if (tmp < thresh[j] && !done)
308
      {
309
310
         if (tmp >= alloc_floor)
            tmp = alloc_floor;
311
312
         else
            tmp = 0;
313
314
      } else
         done = 1;
315
      /* Don't allocate more than we can actually use */
316
      tmp = IMIN(tmp, cap[j]);
317
318
      bits[j] = tmp;
      psum += tmp;
319
   }
320

321
322
   /* Decide which bands to skip, working backwards from the end. */
   for (codedBands=end;;codedBands--)
323
   {
324
325
326
      int band_width;
      int band_bits;
      int rem;
327
      j = codedBands-1;
328
329
330
331
332
333
      /* Never skip the first band, nor a band that has been boosted by
          dynalloc.
         In the first case, we'd be coding a bit to signal we're going to waste
          all the other bits.
         In the second case, we'd be coding a bit to redistribute all the bits
          we just signaled should be cocentrated in this band. */
334
      if (j<=skip_start)
335
      {
336
337
         /* Give the bit we reserved to end skipping back. */
         total += skip_rsv;
338
         break;
339
      }
340
341
342
343
344
      /*Figure out how many left-over bits we would be adding to this band.
        This can include bits we've stolen back from higher, skipped bands.*/
      left = total-psum;
      percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
      left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
345
      rem = IMAX(left-(m->eBands[j]-m->eBands[start]),0);
346
      band_width = m->eBands[codedBands]-m->eBands[j];
Timothy B. Terriberry's avatar
Timothy B. Terriberry committed
347
      band_bits = (int)(bits[j] + percoeff*band_width + rem);
348
349
      /*Only code a skip decision if we're above the threshold for this band.
        Otherwise it is force-skipped.
350
351
        This ensures that we have enough bits to code the skip flag.*/
      if (band_bits >= IMAX(thresh[j], alloc_floor+(1<<BITRES)))
352
      {
353
         if (encode)
354
         {
355
356
357
            /*This if() block is the only part of the allocation function that
               is not a mandatory part of the bitstream: any bands we choose to
               skip here must be explicitly signaled.*/
358
359
360
            /*Choose a threshold with some hysteresis to keep bands from
               fluctuating in and out.*/
            if (band_bits > ((j<prev?7:9)*band_width<<LM<<BITRES)>>4)
361
            {
362
               ec_enc_bit_logp(ec, 1, 1);
363
               break;
364
            }
365
366
            ec_enc_bit_logp(ec, 0, 1);
         } else if (ec_dec_bit_logp(ec, 1)) {
367
            break;
368
369
         }
         /*We used a bit to skip this band.*/
370
371
         psum += 1<<BITRES;
         band_bits -= 1<<BITRES;
372
      }
373
      /*Reclaim the bits originally allocated to this band.*/
374
375
376
377
      psum -= bits[j]+intensity_rsv;
      if (intensity_rsv > 0)
         intensity_rsv = LOG2_FRAC_TABLE[j-start];
      psum += intensity_rsv;
378
      if (band_bits >= alloc_floor)
379
      {
380
         /*If we have enough for a fine energy bit per channel, use it.*/
381
         psum += alloc_floor;
382
         bits[j] = alloc_floor;
383
      } else {
384
         /*Otherwise this band gets nothing at all.*/
385
         bits[j] = 0;
386
      }
387
   }
388

389
390
391
392
393
   celt_assert(codedBands > start);
   /* Code the intensity and dual stereo parameters. */
   if (intensity_rsv > 0)
   {
      if (encode)
394
      {
395
         *intensity = IMIN(*intensity, codedBands);
396
         ec_enc_uint(ec, *intensity-start, codedBands+1-start);
397
      }
398
      else
399
         *intensity = start+ec_dec_uint(ec, codedBands+1-start);
400
401
402
403
404
405
406
407
408
409
410
   }
   else
      *intensity = 0;
   if (*intensity <= start)
   {
      total += dual_stereo_rsv;
      dual_stereo_rsv = 0;
   }
   if (dual_stereo_rsv > 0)
   {
      if (encode)
411
         ec_enc_bit_logp(ec, *dual_stereo, 1);
412
      else
413
         *dual_stereo = ec_dec_bit_logp(ec, 1);
414
415
416
417
418
419
420
421
422
   }
   else
      *dual_stereo = 0;

   /* Allocate the remaining bits */
   left = total-psum;
   percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
   left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
   for (j=start;j<codedBands;j++)
423
      bits[j] += ((int)percoeff*(m->eBands[j+1]-m->eBands[j]));
424
425
   for (j=start;j<codedBands;j++)
   {
Timothy B. Terriberry's avatar
Timothy B. Terriberry committed
426
      int tmp = (int)IMIN(left, m->eBands[j+1]-m->eBands[j]);
427
428
      bits[j] += tmp;
      left -= tmp;
429
   }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
430
   /*for (j=0;j<end;j++)printf("%d ", bits[j]);printf("\n");*/
431
432

   balance = 0;
433
   for (j=start;j<codedBands;j++)
434
   {
435
      int N0, N, den;
436
      int offset;
437
      int NClogN;
438
      int excess;
439

440
      celt_assert(bits[j] >= 0);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
441
      N0 = m->eBands[j+1]-m->eBands[j];
442
      N=N0<<LM;
443
      bits[j] += balance;
444

445
446
      if (N>1)
      {
447
448
449
         excess = IMAX(bits[j]-cap[j],0);
         bits[j] -= excess;

450
         /* Compensate for the extra DoF in stereo */
451
         den=(C*N+ ((C==2 && N>2 && !*dual_stereo && j<*intensity) ? 1 : 0));
452

453
454
         NClogN = den*(m->logN[j] + logM);

455
456
         /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
            compared to their "fair share" of total/N */
457
         offset = (NClogN>>1)-den*FINE_OFFSET;
458

459
460
         /* N=2 is the only point that doesn't match the curve */
         if (N==2)
461
            offset += den<<BITRES>>2;
462

463
464
465
466
467
468
         /* Changing the offset for allocating the second and third
             fine energy bit */
         if (bits[j] + offset < den*2<<BITRES)
            offset += NClogN>>2;
         else if (bits[j] + offset < den*3<<BITRES)
            offset += NClogN>>3;
469

470
471
         /* Divide with rounding */
         ebits[j] = IMAX(0, (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES));
472

473
474
         /* Make sure not to bust */
         if (C*ebits[j] > (bits[j]>>BITRES))
475
            ebits[j] = bits[j] >> stereo >> BITRES;
476

477
478
         /* More than that is useless because that's about as far as PVQ can go */
         ebits[j] = IMIN(ebits[j], MAX_FINE_BITS);
479

480
481
482
483
         /* If we rounded down or capped this band, make it a candidate for the
             final fine energy pass */
         fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;

484
485
486
         /* Remove the allocated fine bits; the rest are assigned to PVQ */
         bits[j] -= C*ebits[j]<<BITRES;

487
488
      } else {
         /* For N=1, all bits go to fine energy except for a single sign bit */
489
490
491
492
         excess = IMAX(0,bits[j]-(C<<BITRES));
         bits[j] -= excess;
         ebits[j] = 0;
         fine_priority[j] = 1;
493
      }
494

495
496
497
498
      /* Fine energy can't take advantage of the re-balancing in
          quant_all_bands().
         Instead, do the re-balancing here.*/
      if(excess > 0)
499
      {
500
501
502
503
504
505
506
         int extra_fine;
         int extra_bits;
         extra_fine = IMIN(excess >> stereo+BITRES, MAX_FINE_BITS-ebits[j]);
         ebits[j] += extra_fine;
         extra_bits = extra_fine*C<<BITRES;
         fine_priority[j] = extra_bits >= excess-balance;
         excess -= extra_bits;
507
      }
508
      balance = excess;
509

Jean-Marc Valin's avatar
Jean-Marc Valin committed
510
      celt_assert(bits[j] >= 0);
511
      celt_assert(ebits[j] >= 0);
512
   }
513
514
515
   /* Save any remaining bits over the cap for the rebalancing in
       quant_all_bands(). */
   *_balance = balance;
516

517
518
519
   /* The skipped bands use all their bits for fine energy. */
   for (;j<end;j++)
   {
520
      ebits[j] = bits[j] >> stereo >> BITRES;
521
522
      celt_assert(C*ebits[j]<<BITRES == bits[j]);
      bits[j] = 0;
523
      fine_priority[j] = ebits[j]<1;
524
   }
525
   RESTORE_STACK;
526
   return codedBands;
527
528
}

529
int compute_allocation(const CELTMode *m, int start, int end, const int *offsets, const int *cap, int alloc_trim, int *intensity, int *dual_stereo,
530
      int total, celt_int32 *balance, int *pulses, int *ebits, int *fine_priority, int _C, int LM, ec_ctx *ec, int encode, int prev)
531
{
532
   int lo, hi, len, j;
533
   const int C = CHANNELS(_C);
534
   int codedBands;
535
   int skip_start;
536
   int skip_rsv;
537
538
   int intensity_rsv;
   int dual_stereo_rsv;
539
540
   VARDECL(int, bits1);
   VARDECL(int, bits2);
541
542
   VARDECL(int, thresh);
   VARDECL(int, trim_offset);
543
   SAVE_STACK;
544
   
545
   total = IMAX(total, 0);
546
   len = m->nbEBands;
547
   skip_start = start;
548
549
550
   /* Reserve a bit to signal the end of manually skipped bands. */
   skip_rsv = total >= 1<<BITRES ? 1<<BITRES : 0;
   total -= skip_rsv;
551
552
553
554
555
556
557
558
559
560
561
562
563
564
   /* Reserve bits for the intensity and dual stereo parameters. */
   intensity_rsv = dual_stereo_rsv = 0;
   if (C==2)
   {
      intensity_rsv = LOG2_FRAC_TABLE[end-start];
      if (intensity_rsv>total)
         intensity_rsv = 0;
      else
      {
         total -= intensity_rsv;
         dual_stereo_rsv = total>=1<<BITRES ? 1<<BITRES : 0;
         total -= dual_stereo_rsv;
      }
   }
565
566
   ALLOC(bits1, len, int);
   ALLOC(bits2, len, int);
567
568
569
570
   ALLOC(thresh, len, int);
   ALLOC(trim_offset, len, int);

   for (j=start;j<end;j++)
571
572
   {
      /* Below this threshold, we're sure not to allocate any PVQ bits */
573
      thresh[j] = IMAX((C)<<BITRES, (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4);
574
      /* Tilt of the allocation curve */
575
      trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(m->nbEBands-j-1)
576
            <<(LM+BITRES)>>6;
577
578
579
580
581
      /* Giving less resolution to single-coefficient bands because they get
         more benefit from having one coarse value per coefficient*/
      if ((m->eBands[j+1]-m->eBands[j])<<LM==1)
         trim_offset[j] -= C<<BITRES;
   }
582
   lo = 1;
583
   hi = m->nbAllocVectors - 1;
584
   do
585
   {
586
      int done = 0;
587
      int psum = 0;
588
      int mid = (lo+hi) >> 1;
589
      for (j=end;j-->start;)
590
      {
591
         int N = m->eBands[j+1]-m->eBands[j];
592
         bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
593
         if (bits1[j] > 0)
594
            bits1[j] = IMAX(0, bits1[j] + trim_offset[j]);
595
         bits1[j] += offsets[j];
596
597
598
599
         if (bits1[j] >= thresh[j] || done)
         {
            done = 1;
            /* Don't allocate more than we can actually use */
600
            psum += IMIN(bits1[j], cap[j]);
601
602
603
604
         } else {
            if (bits1[j] >= C<<BITRES)
               psum += C<<BITRES;
         }
605
      }
606
      if (psum > total)
607
         hi = mid - 1;
608
      else
609
         lo = mid + 1;
610
      /*printf ("lo = %d, hi = %d\n", lo, hi);*/
611
   }
612
613
   while (lo <= hi);
   hi = lo--;
614
   /*printf ("interp between %d and %d\n", lo, hi);*/
615
   for (j=start;j<end;j++)
616
   {
617
      int N = m->eBands[j+1]-m->eBands[j];
618
      bits1[j] = C*N*m->allocVectors[lo*len+j]<<LM>>2;
619
      bits2[j] = hi>=m->nbAllocVectors ?
620
            cap[j] : C*N*m->allocVectors[hi*len+j]<<LM>>2;
621
      if (bits1[j] > 0)
622
623
624
625
626
627
         bits1[j] = IMAX(0, bits1[j] + trim_offset[j]);
      if (bits2[j] > 0)
         bits2[j] = IMAX(0, bits2[j] + trim_offset[j]);
      if (lo > 0)
         bits1[j] += offsets[j];
      bits2[j] += offsets[j];
628
629
      if (offsets[j]>0)
         skip_start = j;
630
      bits2[j] = IMAX(0,bits2[j]-bits1[j]);
631
   }
632
   codedBands = interp_bits2pulses(m, start, end, skip_start, bits1, bits2, thresh, cap,
633
         total, balance, skip_rsv, intensity, intensity_rsv, dual_stereo, dual_stereo_rsv,
634
         pulses, ebits, fine_priority, C, LM, ec, encode, prev);
635
   RESTORE_STACK;
636
   return codedBands;
637
638
}