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