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