Allocation table stored with BITRES accuracy based on the short block size
[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_FRAME_SIZE:
65          *value = mode->mdctSize;
66          break;
67       case CELT_GET_LOOKAHEAD:
68          *value = mode->overlap;
69          break;
70       case CELT_GET_BITSTREAM_VERSION:
71          *value = CELT_BITSTREAM_VERSION;
72          break;
73       case CELT_GET_SAMPLE_RATE:
74          *value = mode->Fs;
75          break;
76       default:
77          return CELT_UNIMPLEMENTED;
78    }
79    return CELT_OK;
80 }
81
82 #ifndef STATIC_MODES
83
84 /* Defining 25 critical bands for the full 0-20 kHz audio bandwidth
85    Taken from http://ccrma.stanford.edu/~jos/bbt/Bark_Frequency_Scale.html */
86 #define BARK_BANDS 25
87 static const celt_int16 bark_freq[BARK_BANDS+1] = {
88       0,   100,   200,   300,   400,
89     510,   630,   770,   920,  1080,
90    1270,  1480,  1720,  2000,  2320,
91    2700,  3150,  3700,  4400,  5300,
92    6400,  7700,  9500, 12000, 15500,
93   20000};
94
95 /* This allocation table is per critical band. When creating a mode, the bits get added together 
96    into the codec bands, which are sometimes larger than one critical band at low frequency */
97
98 #ifdef STDIN_TUNING
99 int BITALLOC_SIZE;
100 int *band_allocation;
101 #else
102 #define BITALLOC_SIZE 12
103 static const int band_allocation[BARK_BANDS*BITALLOC_SIZE] = 
104    /* 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  */
105    {  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*/
106       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*/
107       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*/
108       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*/
109       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*/
110       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*/
111       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*/
112       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*/
113       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*/
114       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*/
115       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*/
116       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*/
117    };
118 #endif
119
120 static celt_int16 *compute_ebands(celt_int32 Fs, int frame_size, int res, int *nbEBands)
121 {
122    celt_int16 *eBands;
123    int i, lin, low, high, nBark, offset=0;
124
125    /* Find the number of critical bands supported by our sampling rate */
126    for (nBark=1;nBark<BARK_BANDS;nBark++)
127     if (bark_freq[nBark+1]*2 >= Fs)
128        break;
129
130    /* Find where the linear part ends (i.e. where the spacing is more than min_width */
131    for (lin=0;lin<nBark;lin++)
132       if (bark_freq[lin+1]-bark_freq[lin] >= res)
133          break;
134
135    low = (bark_freq[lin]+res/2)/res;
136    high = nBark-lin;
137    *nbEBands = low+high;
138    eBands = celt_alloc(sizeof(celt_int16)*(*nbEBands+2));
139    
140    if (eBands==NULL)
141       return NULL;
142    
143    /* Linear spacing (min_width) */
144    for (i=0;i<low;i++)
145       eBands[i] = i;
146    if (low>0)
147       offset = eBands[low-1]*res - bark_freq[lin-1];
148    /* Spacing follows critical bands */
149    for (i=0;i<high;i++)
150    {
151       int target = bark_freq[lin+i];
152       eBands[i+low] = (target+(offset+res)/2)/res;
153       offset = eBands[i+low]*res - target;
154    }
155    /* Enforce the minimum spacing at the boundary */
156    for (i=0;i<*nbEBands;i++)
157       if (eBands[i] < i)
158          eBands[i] = i;
159    eBands[*nbEBands] = (bark_freq[nBark]+res/2)/res;
160    eBands[*nbEBands+1] = frame_size;
161    if (eBands[*nbEBands] > eBands[*nbEBands+1])
162       eBands[*nbEBands] = eBands[*nbEBands+1];
163    for (i=1;i<*nbEBands-1;i++)
164    {
165       if (eBands[i+1]-eBands[i] < eBands[i]-eBands[i-1])
166       {
167          eBands[i] -= (2*eBands[i]-eBands[i-1]-eBands[i+1])/2;
168       }
169    }
170    /*for (i=0;i<=*nbEBands+1;i++)
171       printf ("%d ", eBands[i]);
172    printf ("\n");
173    exit(1);*/
174    /* FIXME: Remove last band if too small */
175    return eBands;
176 }
177
178 static void compute_allocation_table(CELTMode *mode, int res)
179 {
180    int i, j, nBark;
181    celt_int16 *allocVectors;
182
183    /* Find the number of critical bands supported by our sampling rate */
184    for (nBark=1;nBark<BARK_BANDS;nBark++)
185     if (bark_freq[nBark+1]*2 >= mode->Fs)
186        break;
187
188    mode->nbAllocVectors = BITALLOC_SIZE;
189    allocVectors = celt_alloc(sizeof(celt_int16)*(BITALLOC_SIZE*mode->nbEBands));
190    if (allocVectors==NULL)
191       return;
192    /* Compute per-codec-band allocation from per-critical-band matrix */
193    for (i=0;i<BITALLOC_SIZE;i++)
194    {
195       celt_int32 current = 0;
196       int eband = 0;
197       for (j=0;j<nBark;j++)
198       {
199          int edge, low, high;
200          celt_int32 alloc;
201
202          alloc = mode->shortMdctSize*band_allocation[i*BARK_BANDS+j];
203          low = bark_freq[j];
204          high = bark_freq[j+1];
205
206          edge = mode->eBands[eband+1]*res;
207          while (edge <= high && eband < mode->nbEBands)
208          {
209             celt_int32 num;
210             int den, bits;
211             num = alloc * (edge-low);
212             den = high-low;
213             /* Divide with rounding */
214             bits = (2*num+den)/(2*den);
215             allocVectors[i*mode->nbEBands+eband] = (current+bits+128)>>(8-BITRES);
216             /* Remove one bit from every band -- FIXME: this is just a temporary hack*/
217             allocVectors[i*mode->nbEBands+eband] -= 1<<BITRES;
218             if (allocVectors[i*mode->nbEBands+eband]<0)
219                allocVectors[i*mode->nbEBands+eband]=0;
220             /* Remove the part of the band we just allocated */
221             low = edge;
222             alloc -= bits;
223
224             /* Move to next eband */
225             current = 0;
226             eband++;
227             edge = mode->eBands[eband+1]*res;
228          }
229          current += alloc;
230       }
231       if (eband < mode->nbEBands)
232       {
233          allocVectors[i*mode->nbEBands+eband] = (current+128)>>(8-BITRES);
234          /* Same hack as above FIXME: again */
235          allocVectors[i*mode->nbEBands+eband] -= 1<<BITRES;
236          if (allocVectors[i*mode->nbEBands+eband]<0)
237             allocVectors[i*mode->nbEBands+eband]=0;
238       }
239    }
240    mode->allocVectors = allocVectors;
241 }
242
243 #endif /* STATIC_MODES */
244
245 CELTMode *celt_mode_create(celt_int32 Fs, int frame_size, int *error)
246 {
247    int i;
248 #ifdef STDIN_TUNING
249    scanf("%d ", &MIN_BINS);
250    scanf("%d ", &BITALLOC_SIZE);
251    band_allocation = celt_alloc(sizeof(int)*BARK_BANDS*BITALLOC_SIZE);
252    for (i=0;i<BARK_BANDS*BITALLOC_SIZE;i++)
253    {
254       scanf("%d ", band_allocation+i);
255    }
256 #endif
257 #ifdef STATIC_MODES
258    const CELTMode *m = NULL;
259    CELTMode *mode=NULL;
260    ALLOC_STACK;
261 #if !defined(VAR_ARRAYS) && !defined(USE_ALLOCA)
262    if (global_stack==NULL)
263    {
264       celt_free(global_stack);
265       goto failure;
266    }
267 #endif 
268    for (i=0;i<TOTAL_MODES;i++)
269    {
270       if (Fs == static_mode_list[i]->Fs &&
271           frame_size == static_mode_list[i]->mdctSize)
272       {
273          m = static_mode_list[i];
274          break;
275       }
276    }
277    if (m == NULL)
278    {
279       celt_warning("Mode not included as part of the static modes");
280       if (error)
281          *error = CELT_BAD_ARG;
282       return NULL;
283    }
284    mode = (CELTMode*)celt_alloc(sizeof(CELTMode));
285    if (mode==NULL)
286       goto failure;
287    CELT_COPY(mode, m, 1);
288    mode->marker_start = MODEPARTIAL;
289 #else
290    int res;
291    CELTMode *mode=NULL;
292    celt_word16 *window;
293    celt_int16 *logN;
294    ALLOC_STACK;
295 #if !defined(VAR_ARRAYS) && !defined(USE_ALLOCA)
296    if (global_stack==NULL)
297    {
298       celt_free(global_stack);
299       goto failure;
300    }
301 #endif 
302
303    /* The good thing here is that permutation of the arguments will automatically be invalid */
304    
305    if (Fs < 32000 || Fs > 96000)
306    {
307       celt_warning("Sampling rate must be between 32 kHz and 96 kHz");
308       if (error)
309          *error = CELT_BAD_ARG;
310       return NULL;
311    }
312    if (frame_size < 64 || frame_size > 1024 || frame_size%2!=0)
313    {
314       celt_warning("Only even frame sizes from 64 to 1024 are supported");
315       if (error)
316          *error = CELT_BAD_ARG;
317       return NULL;
318    }
319    
320    mode = celt_alloc(sizeof(CELTMode));
321    if (mode==NULL)
322       goto failure;
323    mode->marker_start = MODEPARTIAL;
324    mode->Fs = Fs;
325    mode->mdctSize = frame_size;
326    mode->ePredCoef = QCONST16(.8f,15);
327
328    if (frame_size > 640 && (frame_size%16)==0)
329    {
330      mode->nbShortMdcts = 8;
331    } else if (frame_size > 384 && (frame_size%8)==0)
332    {
333      mode->nbShortMdcts = 4;
334    } else if (frame_size > 384 && (frame_size%10)==0)
335    {
336      mode->nbShortMdcts = 5;
337    } else if (frame_size > 256 && (frame_size%6)==0)
338    {
339      mode->nbShortMdcts = 3;
340    } else if (frame_size > 256 && (frame_size%8)==0)
341    {
342      mode->nbShortMdcts = 4;
343    } else if (frame_size > 64 && (frame_size%4)==0)
344    {
345      mode->nbShortMdcts = 2;
346    } else if (frame_size > 128 && (frame_size%6)==0)
347    {
348      mode->nbShortMdcts = 3;
349    } else
350    {
351      mode->nbShortMdcts = 1;
352    }
353
354    mode->shortMdctSize = mode->mdctSize/mode->nbShortMdcts;
355    res = (mode->Fs+mode->shortMdctSize)/(2*mode->shortMdctSize);
356
357    mode->eBands = compute_ebands(Fs, mode->shortMdctSize, res, &mode->nbEBands);
358    if (mode->eBands==NULL)
359       goto failure;
360
361    mode->pitchEnd = 4000*(celt_int32)frame_size/Fs;
362    
363    /* Overlap must be divisible by 4 */
364    if (mode->nbShortMdcts > 1)
365       mode->overlap = ((frame_size/mode->nbShortMdcts)>>2)<<2; 
366    else
367       mode->overlap = (frame_size>>3)<<2;
368
369
370    compute_allocation_table(mode, res);
371    if (mode->allocVectors==NULL)
372       goto failure;
373    
374    window = (celt_word16*)celt_alloc(mode->overlap*sizeof(celt_word16));
375    if (window==NULL)
376       goto failure;
377
378 #ifndef FIXED_POINT
379    for (i=0;i<mode->overlap;i++)
380       window[i] = Q15ONE*sin(.5*M_PI* sin(.5*M_PI*(i+.5)/mode->overlap) * sin(.5*M_PI*(i+.5)/mode->overlap));
381 #else
382    for (i=0;i<mode->overlap;i++)
383       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))));
384 #endif
385    mode->window = window;
386
387    mode->bits = (const celt_int16 **)compute_alloc_cache(mode, 1, mode->nbShortMdcts);
388    if (mode->bits==NULL)
389       goto failure;
390
391    logN = (celt_int16*)celt_alloc(mode->nbEBands*sizeof(celt_int16));
392    if (logN==NULL)
393       goto failure;
394
395    for (i=0;i<mode->nbEBands;i++)
396       logN[i] = log2_frac(mode->nbShortMdcts*(mode->eBands[i+1]-mode->eBands[i]), BITRES);
397    mode->logN = logN;
398 #endif /* !STATIC_MODES */
399
400    clt_mdct_init(&mode->mdct, 2*mode->mdctSize);
401
402    clt_mdct_init(&mode->shortMdct, 2*mode->shortMdctSize);
403
404    mode->prob = quant_prob_alloc(mode);
405    if ((mode->mdct.trig==NULL) || (mode->shortMdct.trig==NULL)
406 #ifndef ENABLE_TI_DSPLIB55
407         || (mode->mdct.kfft==NULL) || (mode->shortMdct.kfft==NULL)
408 #endif
409         || (mode->prob==NULL))
410      goto failure;
411
412    mode->marker_start = MODEVALID;
413    mode->marker_end   = MODEVALID;
414    if (error)
415       *error = CELT_OK;
416    return mode;
417 failure: 
418    if (error)
419       *error = CELT_INVALID_MODE;
420    if (mode!=NULL)
421       celt_mode_destroy(mode);
422    return NULL;
423 }
424
425 void celt_mode_destroy(CELTMode *mode)
426 {
427    int i;
428    const celt_int16 *prevPtr = NULL;
429    if (mode == NULL)
430    {
431       celt_warning("NULL passed to celt_mode_destroy");
432       return;
433    }
434
435    if (mode->marker_start == MODEFREED || mode->marker_end == MODEFREED)
436    {
437       celt_warning("Freeing a mode which has already been freed"); 
438       return;
439    }
440
441    if (mode->marker_start != MODEVALID && mode->marker_start != MODEPARTIAL)
442    {
443       celt_warning("This is not a valid CELT mode structure");
444       return;  
445    }
446    mode->marker_start = MODEFREED;
447 #ifndef STATIC_MODES
448    if (mode->bits!=NULL)
449    {
450       for (i=0;i<mode->nbEBands;i++)
451       {
452          if (mode->bits[i] != prevPtr)
453          {
454             prevPtr = mode->bits[i];
455             celt_free((int*)mode->bits[i]);
456           }
457       }
458    }   
459    celt_free((celt_int16**)mode->bits);
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 #endif
467    clt_mdct_clear(&mode->mdct);
468    clt_mdct_clear(&mode->shortMdct);
469    quant_prob_free(mode->prob);
470    mode->marker_end = MODEFREED;
471    celt_free((CELTMode *)mode);
472 }
473
474 int check_mode(const CELTMode *mode)
475 {
476    if (mode==NULL)
477       return CELT_INVALID_MODE;
478    if (mode->marker_start == MODEVALID && mode->marker_end == MODEVALID)
479       return CELT_OK;
480    if (mode->marker_start == MODEFREED || mode->marker_end == MODEFREED)
481       celt_warning("Using a mode that has already been freed");
482    else
483       celt_warning("This is not a valid CELT mode");
484    return CELT_INVALID_MODE;
485 }