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