Re-enables test cases, including the rotation test
[opus.git] / libcelt / rate.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Written by Jean-Marc Valin */
4 /*
5    Redistribution and use in source and binary forms, with or without
6    modification, are permitted provided that the following conditions
7    are met:
8
9    - Redistributions of source code must retain the above copyright
10    notice, this list of conditions and the following disclaimer.
11
12    - Redistributions in binary form must reproduce the above copyright
13    notice, this list of conditions and the following disclaimer in the
14    documentation and/or other materials provided with the distribution.
15
16    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
17    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
18    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
19    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
20    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
21    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
22    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
23    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
24    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
25    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29 #ifdef HAVE_CONFIG_H
30 #include "config.h"
31 #endif
32
33 #include <math.h>
34 #include "modes.h"
35 #include "cwrs.h"
36 #include "arch.h"
37 #include "os_support.h"
38
39 #include "entcode.h"
40 #include "rate.h"
41
42 static const unsigned char LOG2_FRAC_TABLE[24]={
43    0,
44    8,13,
45   16,19,21,23,
46   24,26,27,28,29,30,31,32,
47   32,33,34,34,35,36,36,37,37
48 };
49
50 #ifdef CUSTOM_MODES
51
52 /*Determines if V(N,K) fits in a 32-bit unsigned integer.
53   N and K are themselves limited to 15 bits.*/
54 static int fits_in32(int _n, int _k)
55 {
56    static const opus_int16 maxN[15] = {
57       32767, 32767, 32767, 1476, 283, 109,  60,  40,
58        29,  24,  20,  18,  16,  14,  13};
59    static const opus_int16 maxK[15] = {
60       32767, 32767, 32767, 32767, 1172, 238,  95,  53,
61        36,  27,  22,  18,  16,  15,  13};
62    if (_n>=14)
63    {
64       if (_k>=14)
65          return 0;
66       else
67          return _n <= maxN[_k];
68    } else {
69       return _k <= maxK[_n];
70    }
71 }
72
73 void compute_pulse_cache(CELTMode *m, int LM)
74 {
75    int C;
76    int i;
77    int j;
78    int curr=0;
79    int nbEntries=0;
80    int entryN[100], entryK[100], entryI[100];
81    const opus_int16 *eBands = m->eBands;
82    PulseCache *cache = &m->cache;
83    opus_int16 *cindex;
84    unsigned char *bits;
85    unsigned char *cap;
86
87    cindex = celt_alloc(sizeof(cache->index[0])*m->nbEBands*(LM+2));
88    cache->index = cindex;
89
90    /* Scan for all unique band sizes */
91    for (i=0;i<=LM+1;i++)
92    {
93       for (j=0;j<m->nbEBands;j++)
94       {
95          int k;
96          int N = (eBands[j+1]-eBands[j])<<i>>1;
97          cindex[i*m->nbEBands+j] = -1;
98          /* Find other bands that have the same size */
99          for (k=0;k<=i;k++)
100          {
101             int n;
102             for (n=0;n<m->nbEBands && (k!=i || n<j);n++)
103             {
104                if (N == (eBands[n+1]-eBands[n])<<k>>1)
105                {
106                   cindex[i*m->nbEBands+j] = cindex[k*m->nbEBands+n];
107                   break;
108                }
109             }
110          }
111          if (cache->index[i*m->nbEBands+j] == -1 && N!=0)
112          {
113             int K;
114             entryN[nbEntries] = N;
115             K = 0;
116             while (fits_in32(N,get_pulses(K+1)) && K<MAX_PSEUDO)
117                K++;
118             entryK[nbEntries] = K;
119             cindex[i*m->nbEBands+j] = curr;
120             entryI[nbEntries] = curr;
121
122             curr += K+1;
123             nbEntries++;
124          }
125       }
126    }
127    bits = celt_alloc(sizeof(unsigned char)*curr);
128    cache->bits = bits;
129    cache->size = curr;
130    /* Compute the cache for all unique sizes */
131    for (i=0;i<nbEntries;i++)
132    {
133       unsigned char *ptr = bits+entryI[i];
134       opus_int16 tmp[MAX_PULSES+1];
135       get_required_bits(tmp, entryN[i], get_pulses(entryK[i]), BITRES);
136       for (j=1;j<=entryK[i];j++)
137          ptr[j] = tmp[get_pulses(j)]-1;
138       ptr[0] = entryK[i];
139    }
140
141    /* Compute the maximum rate for each band at which we'll reliably use as
142        many bits as we ask for. */
143    cache->caps = cap = celt_alloc(sizeof(cache->caps[0])*(LM+1)*2*m->nbEBands);
144    for (i=0;i<=LM;i++)
145    {
146       for (C=1;C<=2;C++)
147       {
148          for (j=0;j<m->nbEBands;j++)
149          {
150             int N0;
151             int max_bits;
152             N0 = m->eBands[j+1]-m->eBands[j];
153             /* N=1 bands only have a sign bit and fine bits. */
154             if (N0<<i == 1)
155                max_bits = C*(1+MAX_FINE_BITS)<<BITRES;
156             else
157             {
158                const unsigned char *pcache;
159                opus_int32           num;
160                opus_int32           den;
161                int                  LM0;
162                int                  N;
163                int                  offset;
164                int                  ndof;
165                int                  qb;
166                int                  k;
167                LM0 = 0;
168                /* Even-sized bands bigger than N=2 can be split one more
169                    time. */
170                if (N0 > 2 && !(N0&1))
171                {
172                   N0>>=1;
173                   LM0--;
174                }
175                /* N0=1 bands can't be split down to N<2. */
176                else if (N0 <= 1)
177                {
178                   LM0=IMIN(i,1);
179                   N0<<=LM0;
180                }
181                /* Compute the cost for the lowest-level PVQ of a fully split
182                    band. */
183                pcache = bits + cindex[(LM0+1)*m->nbEBands+j];
184                max_bits = pcache[pcache[0]]+1;
185                /* Add in the cost of coding regular splits. */
186                N = N0;
187                for(k=0;k<i-LM0;k++){
188                   max_bits <<= 1;
189                   /* Offset the number of qtheta bits by log2(N)/2
190                       + QTHETA_OFFSET compared to their "fair share" of
191                       total/N */
192                   offset = (m->logN[j]+(LM0+k<<BITRES)>>1)-QTHETA_OFFSET;
193                   /* The number of qtheta bits we'll allocate if the remainder
194                       is to be max_bits.
195                      The average measured cost for theta is 0.89701 times qb,
196                       approximated here as 459/512. */
197                   num=459*(opus_int32)((2*N-1)*offset+max_bits);
198                   den=((opus_int32)(2*N-1)<<9)-459;
199                   qb = IMIN((num+(den>>1))/den, 57);
200                   celt_assert(qb >= 0);
201                   max_bits += qb;
202                   N <<= 1;
203                }
204                /* Add in the cost of a stereo split, if necessary. */
205                if (C==2)
206                {
207                   max_bits <<= 1;
208                   offset = (m->logN[j]+(i<<BITRES)>>1)-(N==2?QTHETA_OFFSET_TWOPHASE:QTHETA_OFFSET);
209                   ndof = 2*N-1-(N==2);
210                   /* The average measured cost for theta with the step PDF is
211                       0.95164 times qb, approximated here as 487/512. */
212                   num = (N==2?512:487)*(opus_int32)(max_bits+ndof*offset);
213                   den = ((opus_int32)ndof<<9)-(N==2?512:487);
214                   qb = IMIN((num+(den>>1))/den, (N==2?64:61));
215                   celt_assert(qb >= 0);
216                   max_bits += qb;
217                }
218                /* Add the fine bits we'll use. */
219                /* Compensate for the extra DoF in stereo */
220                ndof = C*N + ((C==2 && N>2) ? 1 : 0);
221                /* Offset the number of fine bits by log2(N)/2 + FINE_OFFSET
222                    compared to their "fair share" of total/N */
223                offset = (m->logN[j] + (i<<BITRES)>>1)-FINE_OFFSET;
224                /* N=2 is the only point that doesn't match the curve */
225                if (N==2)
226                   offset += 1<<BITRES>>2;
227                /* The number of fine bits we'll allocate if the remainder is
228                    to be max_bits. */
229                num = max_bits+ndof*offset;
230                den = ndof-1<<BITRES;
231                qb = IMIN((num+(den>>1))/den, MAX_FINE_BITS);
232                celt_assert(qb >= 0);
233                max_bits += C*qb<<BITRES;
234             }
235             max_bits = (4*max_bits/(C*(m->eBands[j+1]-m->eBands[j]<<i)))-64;
236             celt_assert(max_bits >= 0);
237             celt_assert(max_bits < 256);
238             *cap++ = (unsigned char)max_bits;
239          }
240       }
241    }
242 }
243
244 #endif /* CUSTOM_MODES */
245
246 #define ALLOC_STEPS 6
247
248 static inline int interp_bits2pulses(const CELTMode *m, int start, int end, int skip_start,
249       const int *bits1, const int *bits2, const int *thresh, const int *cap, opus_int32 total, opus_int32 *_balance,
250       int skip_rsv, int *intensity, int intensity_rsv, int *dual_stereo, int dual_stereo_rsv, int *bits,
251       int *ebits, int *fine_priority, int _C, int LM, ec_ctx *ec, int encode, int prev)
252 {
253    opus_int32 psum;
254    int lo, hi;
255    int i, j;
256    int logM;
257    const int C = CHANNELS(_C);
258    int stereo;
259    int codedBands=-1;
260    int alloc_floor;
261    opus_int32 left, percoeff;
262    int done;
263    int balance;
264    SAVE_STACK;
265
266    alloc_floor = C<<BITRES;
267    stereo = C>1;
268
269    logM = LM<<BITRES;
270    lo = 0;
271    hi = 1<<ALLOC_STEPS;
272    for (i=0;i<ALLOC_STEPS;i++)
273    {
274       int mid = (lo+hi)>>1;
275       psum = 0;
276       done = 0;
277       for (j=end;j-->start;)
278       {
279          int tmp = bits1[j] + (mid*(opus_int32)bits2[j]>>ALLOC_STEPS);
280          if (tmp >= thresh[j] || done)
281          {
282             done = 1;
283             /* Don't allocate more than we can actually use */
284             psum += IMIN(tmp, cap[j]);
285          } else {
286             if (tmp >= alloc_floor)
287                psum += alloc_floor;
288          }
289       }
290       if (psum > total)
291          hi = mid;
292       else
293          lo = mid;
294    }
295    psum = 0;
296    /*printf ("interp bisection gave %d\n", lo);*/
297    done = 0;
298    for (j=end;j-->start;)
299    {
300       int tmp = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
301       if (tmp < thresh[j] && !done)
302       {
303          if (tmp >= alloc_floor)
304             tmp = alloc_floor;
305          else
306             tmp = 0;
307       } else
308          done = 1;
309       /* Don't allocate more than we can actually use */
310       tmp = IMIN(tmp, cap[j]);
311       bits[j] = tmp;
312       psum += tmp;
313    }
314
315    /* Decide which bands to skip, working backwards from the end. */
316    for (codedBands=end;;codedBands--)
317    {
318       int band_width;
319       int band_bits;
320       int rem;
321       j = codedBands-1;
322       /* Never skip the first band, nor a band that has been boosted by
323           dynalloc.
324          In the first case, we'd be coding a bit to signal we're going to waste
325           all the other bits.
326          In the second case, we'd be coding a bit to redistribute all the bits
327           we just signaled should be cocentrated in this band. */
328       if (j<=skip_start)
329       {
330          /* Give the bit we reserved to end skipping back. */
331          total += skip_rsv;
332          break;
333       }
334       /*Figure out how many left-over bits we would be adding to this band.
335         This can include bits we've stolen back from higher, skipped bands.*/
336       left = total-psum;
337       percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
338       left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
339       rem = IMAX(left-(m->eBands[j]-m->eBands[start]),0);
340       band_width = m->eBands[codedBands]-m->eBands[j];
341       band_bits = (int)(bits[j] + percoeff*band_width + rem);
342       /*Only code a skip decision if we're above the threshold for this band.
343         Otherwise it is force-skipped.
344         This ensures that we have enough bits to code the skip flag.*/
345       if (band_bits >= IMAX(thresh[j], alloc_floor+(1<<BITRES)))
346       {
347          if (encode)
348          {
349             /*This if() block is the only part of the allocation function that
350                is not a mandatory part of the bitstream: any bands we choose to
351                skip here must be explicitly signaled.*/
352             /*Choose a threshold with some hysteresis to keep bands from
353                fluctuating in and out.*/
354             if (band_bits > ((j<prev?7:9)*band_width<<LM<<BITRES)>>4)
355             {
356                ec_enc_bit_logp(ec, 1, 1);
357                break;
358             }
359             ec_enc_bit_logp(ec, 0, 1);
360          } else if (ec_dec_bit_logp(ec, 1)) {
361             break;
362          }
363          /*We used a bit to skip this band.*/
364          psum += 1<<BITRES;
365          band_bits -= 1<<BITRES;
366       }
367       /*Reclaim the bits originally allocated to this band.*/
368       psum -= bits[j]+intensity_rsv;
369       if (intensity_rsv > 0)
370          intensity_rsv = LOG2_FRAC_TABLE[j-start];
371       psum += intensity_rsv;
372       if (band_bits >= alloc_floor)
373       {
374          /*If we have enough for a fine energy bit per channel, use it.*/
375          psum += alloc_floor;
376          bits[j] = alloc_floor;
377       } else {
378          /*Otherwise this band gets nothing at all.*/
379          bits[j] = 0;
380       }
381    }
382
383    celt_assert(codedBands > start);
384    /* Code the intensity and dual stereo parameters. */
385    if (intensity_rsv > 0)
386    {
387       if (encode)
388       {
389          *intensity = IMIN(*intensity, codedBands);
390          ec_enc_uint(ec, *intensity-start, codedBands+1-start);
391       }
392       else
393          *intensity = start+ec_dec_uint(ec, codedBands+1-start);
394    }
395    else
396       *intensity = 0;
397    if (*intensity <= start)
398    {
399       total += dual_stereo_rsv;
400       dual_stereo_rsv = 0;
401    }
402    if (dual_stereo_rsv > 0)
403    {
404       if (encode)
405          ec_enc_bit_logp(ec, *dual_stereo, 1);
406       else
407          *dual_stereo = ec_dec_bit_logp(ec, 1);
408    }
409    else
410       *dual_stereo = 0;
411
412    /* Allocate the remaining bits */
413    left = total-psum;
414    percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
415    left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
416    for (j=start;j<codedBands;j++)
417       bits[j] += ((int)percoeff*(m->eBands[j+1]-m->eBands[j]));
418    for (j=start;j<codedBands;j++)
419    {
420       int tmp = (int)IMIN(left, m->eBands[j+1]-m->eBands[j]);
421       bits[j] += tmp;
422       left -= tmp;
423    }
424    /*for (j=0;j<end;j++)printf("%d ", bits[j]);printf("\n");*/
425
426    balance = 0;
427    for (j=start;j<codedBands;j++)
428    {
429       int N0, N, den;
430       int offset;
431       int NClogN;
432       int excess;
433
434       celt_assert(bits[j] >= 0);
435       N0 = m->eBands[j+1]-m->eBands[j];
436       N=N0<<LM;
437       bits[j] += balance;
438
439       if (N>1)
440       {
441          excess = IMAX(bits[j]-cap[j],0);
442          bits[j] -= excess;
443
444          /* Compensate for the extra DoF in stereo */
445          den=(C*N+ ((C==2 && N>2 && !*dual_stereo && j<*intensity) ? 1 : 0));
446
447          NClogN = den*(m->logN[j] + logM);
448
449          /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
450             compared to their "fair share" of total/N */
451          offset = (NClogN>>1)-den*FINE_OFFSET;
452
453          /* N=2 is the only point that doesn't match the curve */
454          if (N==2)
455             offset += den<<BITRES>>2;
456
457          /* Changing the offset for allocating the second and third
458              fine energy bit */
459          if (bits[j] + offset < den*2<<BITRES)
460             offset += NClogN>>2;
461          else if (bits[j] + offset < den*3<<BITRES)
462             offset += NClogN>>3;
463
464          /* Divide with rounding */
465          ebits[j] = IMAX(0, (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES));
466
467          /* Make sure not to bust */
468          if (C*ebits[j] > (bits[j]>>BITRES))
469             ebits[j] = bits[j] >> stereo >> BITRES;
470
471          /* More than that is useless because that's about as far as PVQ can go */
472          ebits[j] = IMIN(ebits[j], MAX_FINE_BITS);
473
474          /* If we rounded down or capped this band, make it a candidate for the
475              final fine energy pass */
476          fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
477
478          /* Remove the allocated fine bits; the rest are assigned to PVQ */
479          bits[j] -= C*ebits[j]<<BITRES;
480
481       } else {
482          /* For N=1, all bits go to fine energy except for a single sign bit */
483          excess = IMAX(0,bits[j]-(C<<BITRES));
484          bits[j] -= excess;
485          ebits[j] = 0;
486          fine_priority[j] = 1;
487       }
488
489       /* Fine energy can't take advantage of the re-balancing in
490           quant_all_bands().
491          Instead, do the re-balancing here.*/
492       if(excess > 0)
493       {
494          int extra_fine;
495          int extra_bits;
496          extra_fine = IMIN(excess >> stereo+BITRES, MAX_FINE_BITS-ebits[j]);
497          ebits[j] += extra_fine;
498          extra_bits = extra_fine*C<<BITRES;
499          fine_priority[j] = extra_bits >= excess-balance;
500          excess -= extra_bits;
501       }
502       balance = excess;
503
504       celt_assert(bits[j] >= 0);
505       celt_assert(ebits[j] >= 0);
506    }
507    /* Save any remaining bits over the cap for the rebalancing in
508        quant_all_bands(). */
509    *_balance = balance;
510
511    /* The skipped bands use all their bits for fine energy. */
512    for (;j<end;j++)
513    {
514       ebits[j] = bits[j] >> stereo >> BITRES;
515       celt_assert(C*ebits[j]<<BITRES == bits[j]);
516       bits[j] = 0;
517       fine_priority[j] = ebits[j]<1;
518    }
519    RESTORE_STACK;
520    return codedBands;
521 }
522
523 int compute_allocation(const CELTMode *m, int start, int end, const int *offsets, const int *cap, int alloc_trim, int *intensity, int *dual_stereo,
524       opus_int32 total, opus_int32 *balance, int *pulses, int *ebits, int *fine_priority, int _C, int LM, ec_ctx *ec, int encode, int prev)
525 {
526    int lo, hi, len, j;
527    const int C = CHANNELS(_C);
528    int codedBands;
529    int skip_start;
530    int skip_rsv;
531    int intensity_rsv;
532    int dual_stereo_rsv;
533    VARDECL(int, bits1);
534    VARDECL(int, bits2);
535    VARDECL(int, thresh);
536    VARDECL(int, trim_offset);
537    SAVE_STACK;
538
539    total = IMAX(total, 0);
540    len = m->nbEBands;
541    skip_start = start;
542    /* Reserve a bit to signal the end of manually skipped bands. */
543    skip_rsv = total >= 1<<BITRES ? 1<<BITRES : 0;
544    total -= skip_rsv;
545    /* Reserve bits for the intensity and dual stereo parameters. */
546    intensity_rsv = dual_stereo_rsv = 0;
547    if (C==2)
548    {
549       intensity_rsv = LOG2_FRAC_TABLE[end-start];
550       if (intensity_rsv>total)
551          intensity_rsv = 0;
552       else
553       {
554          total -= intensity_rsv;
555          dual_stereo_rsv = total>=1<<BITRES ? 1<<BITRES : 0;
556          total -= dual_stereo_rsv;
557       }
558    }
559    ALLOC(bits1, len, int);
560    ALLOC(bits2, len, int);
561    ALLOC(thresh, len, int);
562    ALLOC(trim_offset, len, int);
563
564    for (j=start;j<end;j++)
565    {
566       /* Below this threshold, we're sure not to allocate any PVQ bits */
567       thresh[j] = IMAX((C)<<BITRES, (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4);
568       /* Tilt of the allocation curve */
569       trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(end-j-1)
570             <<(LM+BITRES)>>6;
571       /* Giving less resolution to single-coefficient bands because they get
572          more benefit from having one coarse value per coefficient*/
573       if ((m->eBands[j+1]-m->eBands[j])<<LM==1)
574          trim_offset[j] -= C<<BITRES;
575    }
576    lo = 1;
577    hi = m->nbAllocVectors - 1;
578    do
579    {
580       int done = 0;
581       int psum = 0;
582       int mid = (lo+hi) >> 1;
583       for (j=end;j-->start;)
584       {
585          int bitsj;
586          int N = m->eBands[j+1]-m->eBands[j];
587          bitsj = C*N*m->allocVectors[mid*len+j]<<LM>>2;
588          if (bitsj > 0)
589             bitsj = IMAX(0, bitsj + trim_offset[j]);
590          bitsj += offsets[j];
591          if (bitsj >= thresh[j] || done)
592          {
593             done = 1;
594             /* Don't allocate more than we can actually use */
595             psum += IMIN(bitsj, cap[j]);
596          } else {
597             if (bitsj >= C<<BITRES)
598                psum += C<<BITRES;
599          }
600       }
601       if (psum > total)
602          hi = mid - 1;
603       else
604          lo = mid + 1;
605       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
606    }
607    while (lo <= hi);
608    hi = lo--;
609    /*printf ("interp between %d and %d\n", lo, hi);*/
610    for (j=start;j<end;j++)
611    {
612       int bits1j, bits2j;
613       int N = m->eBands[j+1]-m->eBands[j];
614       bits1j = C*N*m->allocVectors[lo*len+j]<<LM>>2;
615       bits2j = hi>=m->nbAllocVectors ?
616             cap[j] : C*N*m->allocVectors[hi*len+j]<<LM>>2;
617       if (bits1j > 0)
618          bits1j = IMAX(0, bits1j + trim_offset[j]);
619       if (bits2j > 0)
620          bits2j = IMAX(0, bits2j + trim_offset[j]);
621       if (lo > 0)
622          bits1j += offsets[j];
623       bits2j += offsets[j];
624       if (offsets[j]>0)
625          skip_start = j;
626       bits2j = IMAX(0,bits2j-bits1j);
627       bits1[j] = bits1j;
628       bits2[j] = bits2j;
629    }
630    codedBands = interp_bits2pulses(m, start, end, skip_start, bits1, bits2, thresh, cap,
631          total, balance, skip_rsv, intensity, intensity_rsv, dual_stereo, dual_stereo_rsv,
632          pulses, ebits, fine_priority, C, LM, ec, encode, prev);
633    RESTORE_STACK;
634    return codedBands;
635 }
636