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