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