Add assertions for band size restrictions.
[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=1;i<*nbEBands;i++)
200    {
201       /* Every band must be smaller than the last band. */
202       celt_assert(eBands[i]-eBands[i-1]<=eBands[*nbEBands]-eBands[*nbEBands-1]);
203       /* Each band must be no larger than twice the size of the previous one. */
204       celt_assert(eBands[i+1]-eBands[i]<=2*(eBands[i]-eBands[i-1]));
205    }
206
207    /*for (i=0;i<=*nbEBands+1;i++)
208       printf ("%d ", eBands[i]);
209    printf ("\n");
210    exit(1);*/
211    /* FIXME: Remove last band if too small */
212    return eBands;
213 }
214
215 static void compute_allocation_table(CELTMode *mode)
216 {
217    int i, j;
218    unsigned char *allocVectors;
219    int maxBands = sizeof(eband5ms)/sizeof(eband5ms[0])-1;
220
221    mode->nbAllocVectors = BITALLOC_SIZE;
222    allocVectors = celt_alloc(sizeof(unsigned char)*(BITALLOC_SIZE*mode->nbEBands));
223    if (allocVectors==NULL)
224       return;
225
226    /* Check for standard mode */
227    if (mode->Fs == 400*(celt_int32)mode->shortMdctSize)
228    {
229       for (i=0;i<BITALLOC_SIZE*mode->nbEBands;i++)
230          allocVectors[i] = band_allocation[i];
231       mode->allocVectors = allocVectors;
232       return;
233    }
234    /* If not the standard mode, interpolate */
235    /* Compute per-codec-band allocation from per-critical-band matrix */
236    for (i=0;i<BITALLOC_SIZE;i++)
237    {
238       for (j=0;j<mode->nbEBands;j++)
239       {
240          int k;
241          for (k=0;k<maxBands;k++)
242          {
243             if (400*(celt_int32)eband5ms[k] > mode->eBands[j]*(celt_int32)mode->Fs/mode->shortMdctSize)
244                break;
245          }
246          if (k>mode->nbEBands-1)
247             allocVectors[i*mode->nbEBands+j] = band_allocation[i*maxBands + maxBands-1];
248          else {
249             celt_int32 a0, a1;
250             a1 = mode->eBands[j]*(celt_int32)mode->Fs/mode->shortMdctSize - 400*(celt_int32)eband5ms[k-1];
251             a0 = 400*(celt_int32)eband5ms[k] - mode->eBands[j]*(celt_int32)mode->Fs/mode->shortMdctSize;
252             allocVectors[i*mode->nbEBands+j] = (a0*band_allocation[i*maxBands+k-1]
253                                              + a1*band_allocation[i*maxBands+k])/(a0+a1);
254          }
255       }
256    }
257
258    /*printf ("\n");
259    for (i=0;i<BITALLOC_SIZE;i++)
260    {
261       for (j=0;j<mode->nbEBands;j++)
262          printf ("%d ", allocVectors[i*mode->nbEBands+j]);
263       printf ("\n");
264    }
265    exit(0);*/
266
267    mode->allocVectors = allocVectors;
268 }
269
270 #endif /* CUSTOM_MODES */
271
272 CELTMode *celt_mode_create(celt_int32 Fs, int frame_size, int *error)
273 {
274    int i;
275    int res;
276    CELTMode *mode=NULL;
277    celt_word16 *window;
278    celt_int16 *logN;
279    int LM;
280 #ifdef STDIN_TUNING
281    scanf("%d ", &MIN_BINS);
282    scanf("%d ", &BITALLOC_SIZE);
283    band_allocation = celt_alloc(sizeof(int)*BARK_BANDS*BITALLOC_SIZE);
284    for (i=0;i<BARK_BANDS*BITALLOC_SIZE;i++)
285    {
286       scanf("%d ", band_allocation+i);
287    }
288 #endif
289    ALLOC_STACK;
290 #if !defined(VAR_ARRAYS) && !defined(USE_ALLOCA)
291    if (global_stack==NULL)
292       goto failure;
293 #endif 
294
295 #ifndef CUSTOM_MODES_ONLY
296    for (i=0;i<TOTAL_MODES;i++)
297    {
298       int j;
299       for (j=0;j<4;j++)
300       {
301          if (Fs == static_mode_list[i]->Fs &&
302                (frame_size<<j) == static_mode_list[i]->shortMdctSize*static_mode_list[i]->nbShortMdcts)
303          {
304             if (error)
305                *error = CELT_OK;
306             return (CELTMode*)static_mode_list[i];
307          }
308       }
309    }
310 #endif /* CUSTOM_MODES_ONLY */
311
312 #ifndef CUSTOM_MODES
313    if (error)
314       *error = CELT_BAD_ARG;
315    return NULL;
316 #else
317
318    /* The good thing here is that permutation of the arguments will automatically be invalid */
319    
320    if (Fs < 8000 || Fs > 96000)
321    {
322       if (error)
323          *error = CELT_BAD_ARG;
324       return NULL;
325    }
326    if (frame_size < 40 || frame_size > 1024 || frame_size%2!=0)
327    {
328       if (error)
329          *error = CELT_BAD_ARG;
330       return NULL;
331    }
332    
333    mode = celt_alloc(sizeof(CELTMode));
334    if (mode==NULL)
335       goto failure;
336    mode->Fs = Fs;
337
338    /* Pre/de-emphasis depends on sampling rate. The "standard" pre-emphasis
339       is defined as A(z) = 1 - 0.85*z^-1 at 48 kHz. Other rates should
340       approximate that. */
341    if(Fs < 12000) /* 8 kHz */
342    {
343       mode->preemph[0] =  QCONST16(0.3500061035f, 15);
344       mode->preemph[1] = -QCONST16(0.1799926758f, 15);
345       mode->preemph[2] =  QCONST16(0.2719968125f, SIG_SHIFT); /* exact 1/preemph[3] */
346       mode->preemph[3] =  QCONST16(3.6765136719f, 13);
347    } else if(Fs < 24000) /* 16 kHz */
348    {
349       mode->preemph[0] =  QCONST16(0.6000061035f, 15);
350       mode->preemph[1] = -QCONST16(0.1799926758f, 15);
351       mode->preemph[2] =  QCONST16(0.4424998650f, SIG_SHIFT); /* exact 1/preemph[3] */
352       mode->preemph[3] =  QCONST16(2.2598876953f, 13);
353    } else if(Fs < 40000) /* 32 kHz */
354    {
355       mode->preemph[0] =  QCONST16(0.7799987793f, 15);
356       mode->preemph[1] = -QCONST16(0.1000061035f, 15);
357       mode->preemph[2] =  QCONST16(0.7499771125f, SIG_SHIFT); /* exact 1/preemph[3] */
358       mode->preemph[3] =  QCONST16(1.3333740234f, 13);
359    } else /* 48 kHz */
360    {
361       mode->preemph[0] =  QCONST16(0.8500061035f, 15);
362       mode->preemph[1] =  QCONST16(0.0f, 15);
363       mode->preemph[2] =  QCONST16(1.f, SIG_SHIFT);
364       mode->preemph[3] =  QCONST16(1.f, 13);
365    }
366
367    if ((celt_int32)frame_size*75 >= Fs && (frame_size%16)==0)
368    {
369      LM = 3;
370    } else if ((celt_int32)frame_size*150 >= Fs && (frame_size%8)==0)
371    {
372      LM = 2;
373    } else if ((celt_int32)frame_size*300 >= Fs && (frame_size%4)==0)
374    {
375      LM = 1;
376    } else
377    {
378      LM = 0;
379    }
380
381    mode->maxLM = LM;
382    mode->nbShortMdcts = 1<<LM;
383    mode->shortMdctSize = frame_size/mode->nbShortMdcts;
384    res = (mode->Fs+mode->shortMdctSize)/(2*mode->shortMdctSize);
385
386    mode->eBands = compute_ebands(Fs, mode->shortMdctSize, res, &mode->nbEBands);
387    if (mode->eBands==NULL)
388       goto failure;
389
390    mode->effEBands = mode->nbEBands;
391    while (mode->eBands[mode->effEBands] > mode->shortMdctSize)
392       mode->effEBands--;
393    
394    /* Overlap must be divisible by 4 */
395    mode->overlap = ((mode->shortMdctSize>>2)<<2);
396
397    compute_allocation_table(mode);
398    if (mode->allocVectors==NULL)
399       goto failure;
400    
401    window = (celt_word16*)celt_alloc(mode->overlap*sizeof(celt_word16));
402    if (window==NULL)
403       goto failure;
404
405 #ifndef FIXED_POINT
406    for (i=0;i<mode->overlap;i++)
407       window[i] = Q15ONE*sin(.5*M_PI* sin(.5*M_PI*(i+.5)/mode->overlap) * sin(.5*M_PI*(i+.5)/mode->overlap));
408 #else
409    for (i=0;i<mode->overlap;i++)
410       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))));
411 #endif
412    mode->window = window;
413
414    logN = (celt_int16*)celt_alloc(mode->nbEBands*sizeof(celt_int16));
415    if (logN==NULL)
416       goto failure;
417
418    for (i=0;i<mode->nbEBands;i++)
419       logN[i] = log2_frac(mode->eBands[i+1]-mode->eBands[i], BITRES);
420    mode->logN = logN;
421
422    compute_pulse_cache(mode, mode->maxLM);
423
424    clt_mdct_init(&mode->mdct, 2*mode->shortMdctSize*mode->nbShortMdcts, mode->maxLM);
425    if ((mode->mdct.trig==NULL)
426 #ifndef ENABLE_TI_DSPLIB55
427          || (mode->mdct.kfft==NULL)
428 #endif
429    )
430       goto failure;
431
432    if (error)
433       *error = CELT_OK;
434
435    return mode;
436 failure: 
437    if (error)
438       *error = CELT_INVALID_MODE;
439    if (mode!=NULL)
440       celt_mode_destroy(mode);
441    return NULL;
442 #endif /* !CUSTOM_MODES */
443 }
444
445 void celt_mode_destroy(CELTMode *mode)
446 {
447 #ifdef CUSTOM_MODES
448    int i;
449    if (mode == NULL)
450       return;
451 #ifndef CUSTOM_MODES_ONLY
452    for (i=0;i<TOTAL_MODES;i++)
453    {
454       if (mode == static_mode_list[i])
455       {
456          return;
457       }
458    }
459 #endif /* CUSTOM_MODES_ONLY */
460    celt_free((celt_int16*)mode->eBands);
461    celt_free((celt_int16*)mode->allocVectors);
462    
463    celt_free((celt_word16*)mode->window);
464    celt_free((celt_int16*)mode->logN);
465
466    celt_free((celt_int16*)mode->cache.index);
467    celt_free((unsigned char*)mode->cache.bits);
468    celt_free((unsigned char*)mode->cache.caps);
469    clt_mdct_clear(&mode->mdct);
470
471    celt_free((CELTMode *)mode);
472 #endif
473 }