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