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