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