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