Add coarse energy entropy model tuning.
[opus.git] / libcelt / modes.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Copyright (c) 2008 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    - Neither the name of the Xiph.org Foundation nor the names of its
18    contributors may be used to endorse or promote products derived from
19    this software without specific prior written permission.
20    
21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
25    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 #ifdef HAVE_CONFIG_H
35 #include "config.h"
36 #endif
37
38 #include "celt.h"
39 #include "modes.h"
40 #include "rate.h"
41 #include "os_support.h"
42 #include "stack_alloc.h"
43 #include "quant_bands.h"
44
45 static const celt_int16 eband5ms[] = {
46   0,  1,  2,  3,  4,  5,  6,  7,  8, 10, 12, 14, 16, 20, 24, 28, 34, 40, 48, 60, 78, 100
47 };
48
49 #if 1
50
51 #define BITALLOC_SIZE 9
52 /* Bit allocation table in units of 1/32 bit/sample (0.1875 dB SNR) */
53 static const unsigned char band_allocation[] = {
54 /*0  200 400 600 800  1k 1.2 1.4 1.6  2k 2.4 2.8 3.2  4k 4.8 5.6 6.8  8k 9.6 12k 15.6 */
55   0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
56  95, 90, 80, 75, 65, 60, 55, 50, 44, 40, 35, 30, 15,  1,  0,  0,  0,  0,  0,  0,  0,
57 100, 95, 90, 88, 85, 82, 78, 74, 70, 65, 60, 54, 45, 35, 25, 15,  1,  0,  0,  0,  0,
58 120,110,110,110,100, 96, 90, 88, 84, 76, 70, 65, 60, 45, 35, 25, 20,  1,  1,  0,  0,
59 135,125,125,125,115,112,104,104,100, 96, 83, 78, 70, 55, 46, 36, 32, 28, 20,  8,  0,
60 170,165,157,155,149,145,143,138,138,138,129,124,108, 96, 88, 83, 72, 56, 44, 28,  2,
61 192,192,160,160,160,160,160,160,160,160,150,140,120,110,100, 90, 80, 70, 60, 50, 20,
62 224,224,192,192,192,192,192,192,192,160,160,160,160,160,160,160,160,120, 80, 64, 40,
63 255,255,224,224,224,224,224,224,224,192,192,192,192,192,192,192,192,192,192,192,120,
64 };
65
66 #else
67
68 #if 1
69 /* Alternate tuning based on exp_tuning_knobs from 2010/10/21:
70  87 5 25 14 300
71 120 5 25 16 400
72 125 5 25 17 500
73 130 5 24 18 600
74 130 5 24 19 920
75 140 7 24 20 1240
76 170 7 20 21 2000
77 120 5 20 21 4000
78  */
79
80 #define BITALLOC_SIZE 9
81 /* Bit allocation table in units of 1/32 bit/sample (0.1875 dB SNR) */
82 static const unsigned char band_allocation[] = {
83   0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
84  99, 99, 90, 84, 78, 71, 65, 58, 51, 45, 39, 32, 26, 20,  1,  0,  0,  0,  0,  0,  0,
85 124,124,109, 98, 88, 77, 71, 65, 59, 52, 46, 40, 34, 27, 21, 15,  1,  0,  0,  0,  0,
86 132,132,117,107, 97, 87, 80, 74, 68, 62, 55, 49, 43, 37, 30, 23, 17,  1,  1,  0,  0,
87 138,138,123,112,102, 91, 85, 79, 73, 67, 61, 55, 49, 43, 37, 31, 25, 19,  1,  0,  0,
88 152,152,139,130,120,110,104, 98, 92, 86, 80, 74, 68, 62, 56, 50, 44, 38, 32,  8,  0,
89 167,167,154,145,136,127,117,108,102, 96, 90, 84, 78, 72, 66, 60, 54, 48, 42, 36,  2,
90 215,215,194,180,165,151,137,121,116,111,106,101, 96, 91, 86, 81, 76, 71, 66, 61, 56,
91 240,240,239,233,225,218,211,203,198,193,188,183,178,173,168,163,158,153,148,143,138,
92 };
93
94 #else
95 /* Alternate tuning (partially derived from Vorbis) */
96 #define BITALLOC_SIZE 11
97 /* Bit allocation table in units of 1/32 bit/sample (0.1875 dB SNR) */
98 static const unsigned char band_allocation[] = {
99   0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
100  85, 85, 75, 70, 63, 55, 50, 45, 40, 35, 33, 28, 15,  1,  0,  0,  0,  0,  0,  0,  0,
101 108,105, 96, 88, 83, 77, 72, 66, 59, 53, 49, 44, 41, 37, 30, 27,  1,  0,  0,  0,  0,
102 120,117,108,100, 95, 89, 84, 78, 71, 65, 61, 56, 53, 49, 42, 34, 30, 25,  1,  0,  0,
103 136,131,123,118,109, 99, 93, 87, 81, 75, 69, 66, 61, 56, 50, 45, 40, 35, 32,  1,  1,
104 151,148,138,131,122,113,105,102, 96, 92, 85, 82, 76, 68, 63, 58, 51, 44, 38, 27,  8,
105 171,168,158,147,139,130,123,119,111,108,103, 99, 91, 82, 78, 75, 66, 55, 48, 36, 12,
106 187,184,174,163,155,146,139,135,127,124,119,115,107, 98, 94, 91, 82, 71, 64, 52, 28,
107 203,200,191,181,174,166,159,156,149,147,143,139,132,124,121,119,110,100, 94, 83, 60,
108 219,216,207,197,190,183,176,173,166,164,161,157,150,142,139,138,129,119,113,102, 80,
109 229,229,224,222,223,224,224,225,222,221,221,220,220,219,218,200,178,154,146,130,102,
110 };
111 #endif
112 #endif
113
114 #ifdef STATIC_MODES
115 #ifdef FIXED_POINT
116 #include "static_modes_fixed.c"
117 #else
118 #include "static_modes_float.c"
119 #endif
120 #endif
121
122 #ifndef M_PI
123 #define M_PI 3.141592653
124 #endif
125
126
127 int celt_mode_info(const CELTMode *mode, int request, celt_int32 *value)
128 {
129    switch (request)
130    {
131       case CELT_GET_LOOKAHEAD:
132          *value = mode->overlap;
133          break;
134       case CELT_GET_BITSTREAM_VERSION:
135          *value = CELT_BITSTREAM_VERSION;
136          break;
137       case CELT_GET_SAMPLE_RATE:
138          *value = mode->Fs;
139          break;
140       default:
141          return CELT_UNIMPLEMENTED;
142    }
143    return CELT_OK;
144 }
145
146 #ifndef STATIC_MODES
147
148 /* Defining 25 critical bands for the full 0-20 kHz audio bandwidth
149    Taken from http://ccrma.stanford.edu/~jos/bbt/Bark_Frequency_Scale.html */
150 #define BARK_BANDS 25
151 static const celt_int16 bark_freq[BARK_BANDS+1] = {
152       0,   100,   200,   300,   400,
153     510,   630,   770,   920,  1080,
154    1270,  1480,  1720,  2000,  2320,
155    2700,  3150,  3700,  4400,  5300,
156    6400,  7700,  9500, 12000, 15500,
157   20000};
158
159 static celt_int16 *compute_ebands(celt_int32 Fs, int frame_size, int res, int *nbEBands)
160 {
161    celt_int16 *eBands;
162    int i, lin, low, high, nBark, offset=0;
163
164    /* All modes that have 2.5 ms short blocks use the same definition */
165    if (Fs == 400*(celt_int32)frame_size)
166    {
167       *nbEBands = sizeof(eband5ms)/sizeof(eband5ms[0])-1;
168       eBands = celt_alloc(sizeof(celt_int16)*(*nbEBands+1));
169       for (i=0;i<*nbEBands+1;i++)
170          eBands[i] = eband5ms[i];
171       return eBands;
172    }
173    /* Find the number of critical bands supported by our sampling rate */
174    for (nBark=1;nBark<BARK_BANDS;nBark++)
175     if (bark_freq[nBark+1]*2 >= Fs)
176        break;
177
178    /* Find where the linear part ends (i.e. where the spacing is more than min_width */
179    for (lin=0;lin<nBark;lin++)
180       if (bark_freq[lin+1]-bark_freq[lin] >= res)
181          break;
182
183    low = (bark_freq[lin]+res/2)/res;
184    high = nBark-lin;
185    *nbEBands = low+high;
186    eBands = celt_alloc(sizeof(celt_int16)*(*nbEBands+2));
187    
188    if (eBands==NULL)
189       return NULL;
190    
191    /* Linear spacing (min_width) */
192    for (i=0;i<low;i++)
193       eBands[i] = i;
194    if (low>0)
195       offset = eBands[low-1]*res - bark_freq[lin-1];
196    /* Spacing follows critical bands */
197    for (i=0;i<high;i++)
198    {
199       int target = bark_freq[lin+i];
200       eBands[i+low] = (target+(offset+res)/2)/res;
201       offset = eBands[i+low]*res - target;
202    }
203    /* Enforce the minimum spacing at the boundary */
204    for (i=0;i<*nbEBands;i++)
205       if (eBands[i] < i)
206          eBands[i] = i;
207    eBands[*nbEBands] = (bark_freq[nBark]+res/2)/res;
208    if (eBands[*nbEBands] > frame_size)
209       eBands[*nbEBands] = frame_size;
210    for (i=1;i<*nbEBands-1;i++)
211    {
212       if (eBands[i+1]-eBands[i] < eBands[i]-eBands[i-1])
213       {
214          eBands[i] -= (2*eBands[i]-eBands[i-1]-eBands[i+1])/2;
215       }
216    }
217    /*for (i=0;i<=*nbEBands+1;i++)
218       printf ("%d ", eBands[i]);
219    printf ("\n");
220    exit(1);*/
221    /* FIXME: Remove last band if too small */
222    return eBands;
223 }
224
225 static void compute_allocation_table(CELTMode *mode, int res)
226 {
227    int i, j;
228    unsigned char *allocVectors;
229    int maxBands = sizeof(eband5ms)/sizeof(eband5ms[0])-1;
230
231    mode->nbAllocVectors = BITALLOC_SIZE;
232    allocVectors = celt_alloc(sizeof(unsigned char)*(BITALLOC_SIZE*mode->nbEBands));
233    if (allocVectors==NULL)
234       return;
235
236    /* Check for standard mode */
237    if (mode->Fs == 400*(celt_int32)mode->shortMdctSize && mode->Fs >= 40000)
238    {
239       for (i=0;i<BITALLOC_SIZE*mode->nbEBands;i++)
240          allocVectors[i] = band_allocation[i];
241       mode->allocVectors = allocVectors;
242       return;
243    }
244    /* If not the standard mode, interpolate */
245
246    /* Compute per-codec-band allocation from per-critical-band matrix */
247    for (i=0;i<BITALLOC_SIZE;i++)
248    {
249       celt_int32 current = 0;
250       int eband = 0;
251       /* We may be looping over too many bands, but eband will stop being
252          incremented once we reach the last band */
253       for (j=0;j<maxBands;j++)
254       {
255          int edge, low, high;
256          celt_int32 alloc;
257          alloc = band_allocation[i*maxBands + j]*(mode->eBands[eband+1]-mode->eBands[eband])<<4;
258          low = eband5ms[j]*200;
259          high = eband5ms[j+1]*200;
260          edge = mode->eBands[eband+1]*res;
261          while (edge <= high && eband < mode->nbEBands)
262          {
263             celt_int32 num;
264             int den, bits;
265             int N = (mode->eBands[eband+1]-mode->eBands[eband]);
266             num = alloc * (edge-low);
267             den = high-low;
268             /* Divide with rounding */
269             bits = (2*num+den)/(2*den);
270             allocVectors[i*mode->nbEBands+eband] = (2*(current+bits)+(N<<4))/(2*N<<4);
271             /* Remove the part of the band we just allocated */
272             low = edge;
273             alloc -= bits;
274
275             /* Move to next eband */
276             current = 0;
277             eband++;
278             if (eband < mode->nbEBands)
279                edge = mode->eBands[eband+1]*res;
280          }
281          current += alloc;
282       }
283       if (eband < mode->nbEBands)
284       {
285          int N = (mode->eBands[eband+1]-mode->eBands[eband]);
286          allocVectors[i*mode->nbEBands+eband] = (2*current+(N<<4))/(2*N<<4);
287       }
288    }
289    /*printf ("\n");
290    for (i=0;i<BITALLOC_SIZE;i++)
291    {
292       for (j=0;j<mode->nbEBands;j++)
293          printf ("%d ", allocVectors[i*mode->nbEBands+j]);
294       printf ("\n");
295    }
296    exit(0);*/
297
298    mode->allocVectors = allocVectors;
299 }
300
301 #endif /* STATIC_MODES */
302
303 CELTMode *celt_mode_create(celt_int32 Fs, int frame_size, int *error)
304 {
305    int i;
306 #ifdef STATIC_MODES
307    for (i=0;i<TOTAL_MODES;i++)
308    {
309       if (Fs == static_mode_list[i]->Fs &&
310           frame_size == static_mode_list[i]->shortMdctSize*static_mode_list[i]->nbShortMdcts)
311       {
312          return (CELTMode*)static_mode_list[i];
313       }
314    }
315    if (error)
316       *error = CELT_BAD_ARG;
317    return NULL;
318 #else
319    int res;
320    CELTMode *mode=NULL;
321    celt_word16 *window;
322    celt_int16 *logN;
323    int LM;
324 #ifdef STDIN_TUNING
325    scanf("%d ", &MIN_BINS);
326    scanf("%d ", &BITALLOC_SIZE);
327    band_allocation = celt_alloc(sizeof(int)*BARK_BANDS*BITALLOC_SIZE);
328    for (i=0;i<BARK_BANDS*BITALLOC_SIZE;i++)
329    {
330       scanf("%d ", band_allocation+i);
331    }
332 #endif
333    ALLOC_STACK;
334 #if !defined(VAR_ARRAYS) && !defined(USE_ALLOCA)
335    if (global_stack==NULL)
336       goto failure;
337 #endif 
338
339    /* The good thing here is that permutation of the arguments will automatically be invalid */
340    
341    if (Fs < 8000 || Fs > 96000)
342    {
343       if (error)
344          *error = CELT_BAD_ARG;
345       return NULL;
346    }
347    if (frame_size < 40 || frame_size > 1024 || frame_size%2!=0)
348    {
349       if (error)
350          *error = CELT_BAD_ARG;
351       return NULL;
352    }
353    
354    mode = celt_alloc(sizeof(CELTMode));
355    if (mode==NULL)
356       goto failure;
357    mode->Fs = Fs;
358
359    /* Pre/de-emphasis depends on sampling rate. The "standard" pre-emphasis
360       is defined as A(z) = 1 - 0.85*z^-1 at 48 kHz. Other rates should
361       approximate that. */
362    if(Fs < 12000) /* 8 kHz */
363    {
364       mode->preemph[0] =  QCONST16(.35f, 15);
365       mode->preemph[1] = -QCONST16(.18f, 15);
366       mode->preemph[2] =  QCONST16(.272f, SIG_SHIFT);
367       mode->preemph[3] =  QCONST16(3.6765f, 13);
368    } else if(Fs < 24000) /* 16 kHz */
369    {
370       mode->preemph[0] =  QCONST16(.6f, 15);
371       mode->preemph[1] = -QCONST16(.18f, 15);
372       mode->preemph[2] =  QCONST16(.4425f, SIG_SHIFT);
373       mode->preemph[3] =  QCONST16(2.259887f, 13);
374    } else if(Fs < 40000) /* 32 kHz */
375    {
376       mode->preemph[0] =  QCONST16(.78f, 15);
377       mode->preemph[1] = -QCONST16(.1f, 15);
378       mode->preemph[2] =  QCONST16(.75f, SIG_SHIFT);
379       mode->preemph[3] =  QCONST16(1.33333333f, 13);
380    } else /* 48 kHz */
381    {
382       mode->preemph[0] =  QCONST16(.85f, 15);
383       mode->preemph[1] =  QCONST16(.0f, 15);
384       mode->preemph[2] =  QCONST16(1.f, SIG_SHIFT);
385       mode->preemph[3] =  QCONST16(1.f, 13);
386    }
387
388    if ((celt_int32)frame_size*75 >= Fs && (frame_size%16)==0)
389    {
390      LM = 3;
391    } else if ((celt_int32)frame_size*150 >= Fs && (frame_size%8)==0)
392    {
393      LM = 2;
394    } else if ((celt_int32)frame_size*300 >= Fs && (frame_size%4)==0)
395    {
396      LM = 1;
397    } else
398    {
399      LM = 0;
400    }
401
402    mode->maxLM = LM;
403    mode->nbShortMdcts = 1<<LM;
404    mode->shortMdctSize = frame_size/mode->nbShortMdcts;
405    res = (mode->Fs+mode->shortMdctSize)/(2*mode->shortMdctSize);
406
407    mode->eBands = compute_ebands(Fs, mode->shortMdctSize, res, &mode->nbEBands);
408    if (mode->eBands==NULL)
409       goto failure;
410
411    mode->effEBands = mode->nbEBands;
412    while (mode->eBands[mode->effEBands] > mode->shortMdctSize)
413       mode->effEBands--;
414    
415    /* Overlap must be divisible by 4 */
416    if (mode->nbShortMdcts > 1)
417       mode->overlap = (mode->shortMdctSize>>2)<<2;
418    else
419       mode->overlap = (frame_size>>3)<<2;
420
421
422    compute_allocation_table(mode, res);
423    if (mode->allocVectors==NULL)
424       goto failure;
425    
426    window = (celt_word16*)celt_alloc(mode->overlap*sizeof(celt_word16));
427    if (window==NULL)
428       goto failure;
429
430 #ifndef FIXED_POINT
431    for (i=0;i<mode->overlap;i++)
432       window[i] = Q15ONE*sin(.5*M_PI* sin(.5*M_PI*(i+.5)/mode->overlap) * sin(.5*M_PI*(i+.5)/mode->overlap));
433 #else
434    for (i=0;i<mode->overlap;i++)
435       window[i] = MIN32(32767,floor(.5+32768.*sin(.5*M_PI* sin(.5*M_PI*(i+.5)/mode->overlap) * sin(.5*M_PI*(i+.5)/mode->overlap))));
436 #endif
437    mode->window = window;
438
439    logN = (celt_int16*)celt_alloc(mode->nbEBands*sizeof(celt_int16));
440    if (logN==NULL)
441       goto failure;
442
443    for (i=0;i<mode->nbEBands;i++)
444       logN[i] = log2_frac(mode->eBands[i+1]-mode->eBands[i], BITRES);
445    mode->logN = logN;
446
447    compute_pulse_cache(mode, mode->maxLM);
448
449    clt_mdct_init(&mode->mdct, 2*mode->shortMdctSize*mode->nbShortMdcts, mode->maxLM);
450    if ((mode->mdct.trig==NULL)
451 #ifndef ENABLE_TI_DSPLIB55
452          || (mode->mdct.kfft==NULL)
453 #endif
454    )
455       goto failure;
456
457    if (error)
458       *error = CELT_OK;
459
460    return mode;
461 failure: 
462    if (error)
463       *error = CELT_INVALID_MODE;
464    if (mode!=NULL)
465       celt_mode_destroy(mode);
466    return NULL;
467 #endif /* !STATIC_MODES */
468 }
469
470 void celt_mode_destroy(CELTMode *mode)
471 {
472 #ifndef STATIC_MODES
473    if (mode == NULL)
474       return;
475
476    celt_free((celt_int16*)mode->eBands);
477    celt_free((celt_int16*)mode->allocVectors);
478    
479    celt_free((celt_word16*)mode->window);
480    celt_free((celt_int16*)mode->logN);
481
482    celt_free((celt_int16*)mode->cache.index);
483    celt_free((unsigned char*)mode->cache.bits);
484    clt_mdct_clear(&mode->mdct);
485
486    celt_free((CELTMode *)mode);
487 #endif
488 }