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