Max delta: +/- 16384
[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       eBands[i+low] = (target+(offset+res)/2)/res;
175       offset = eBands[i+low]*res - target;
176    }
177    /* Enforce the minimum spacing at the boundary */
178    for (i=0;i<*nbEBands;i++)
179       if (eBands[i] < i)
180          eBands[i] = i;
181    eBands[*nbEBands] = (bark_freq[nBark]+res/2)/res;
182    if (eBands[*nbEBands] > frame_size)
183       eBands[*nbEBands] = frame_size;
184    for (i=1;i<*nbEBands-1;i++)
185    {
186       if (eBands[i+1]-eBands[i] < eBands[i]-eBands[i-1])
187       {
188          eBands[i] -= (2*eBands[i]-eBands[i-1]-eBands[i+1])/2;
189       }
190    }
191    /*for (i=0;i<=*nbEBands+1;i++)
192       printf ("%d ", eBands[i]);
193    printf ("\n");
194    exit(1);*/
195    /* FIXME: Remove last band if too small */
196    return eBands;
197 }
198
199 static void compute_allocation_table(CELTMode *mode, int res)
200 {
201    int i, j;
202    unsigned char *allocVectors;
203    int maxBands = sizeof(eband5ms)/sizeof(eband5ms[0])-1;
204
205    mode->nbAllocVectors = BITALLOC_SIZE;
206    allocVectors = celt_alloc(sizeof(unsigned char)*(BITALLOC_SIZE*mode->nbEBands));
207    if (allocVectors==NULL)
208       return;
209
210    /* Check for standard mode */
211    if (mode->Fs == 400*(celt_int32)mode->shortMdctSize)
212    {
213       for (i=0;i<BITALLOC_SIZE*mode->nbEBands;i++)
214          allocVectors[i] = band_allocation[i];
215       mode->allocVectors = allocVectors;
216       return;
217    }
218    /* If not the standard mode, interpolate */
219    /* Compute per-codec-band allocation from per-critical-band matrix */
220    for (i=0;i<BITALLOC_SIZE;i++)
221    {
222       for (j=0;j<mode->nbEBands;j++)
223       {
224          int k;
225          for (k=0;k<maxBands;k++)
226          {
227             if (400*(celt_int32)eband5ms[k] > mode->eBands[j]*(celt_int32)mode->Fs/mode->shortMdctSize)
228                break;
229          }
230          if (k>mode->nbEBands-1)
231             allocVectors[i*mode->nbEBands+j] = band_allocation[i*maxBands + maxBands-1];
232          else {
233             celt_int32 a0, a1;
234             a1 = mode->eBands[j]*(celt_int32)mode->Fs/mode->shortMdctSize - 400*(celt_int32)eband5ms[k-1];
235             a0 = 400*(celt_int32)eband5ms[k] - mode->eBands[j]*(celt_int32)mode->Fs/mode->shortMdctSize;
236             allocVectors[i*mode->nbEBands+j] = (a0*band_allocation[i*maxBands+k-1]
237                                              + a1*band_allocation[i*maxBands+k])/(a0+a1);
238          }
239       }
240    }
241
242    /*printf ("\n");
243    for (i=0;i<BITALLOC_SIZE;i++)
244    {
245       for (j=0;j<mode->nbEBands;j++)
246          printf ("%d ", allocVectors[i*mode->nbEBands+j]);
247       printf ("\n");
248    }
249    exit(0);*/
250
251    mode->allocVectors = allocVectors;
252 }
253
254 #endif /* STATIC_MODES */
255
256 CELTMode *celt_mode_create(celt_int32 Fs, int frame_size, int *error)
257 {
258    int i;
259 #ifdef STATIC_MODES
260    for (i=0;i<TOTAL_MODES;i++)
261    {
262       if (Fs == static_mode_list[i]->Fs &&
263           frame_size == static_mode_list[i]->shortMdctSize*static_mode_list[i]->nbShortMdcts)
264       {
265          return (CELTMode*)static_mode_list[i];
266       }
267    }
268    if (error)
269       *error = CELT_BAD_ARG;
270    return NULL;
271 #else
272    int res;
273    CELTMode *mode=NULL;
274    celt_word16 *window;
275    celt_int16 *logN;
276    int LM;
277 #ifdef STDIN_TUNING
278    scanf("%d ", &MIN_BINS);
279    scanf("%d ", &BITALLOC_SIZE);
280    band_allocation = celt_alloc(sizeof(int)*BARK_BANDS*BITALLOC_SIZE);
281    for (i=0;i<BARK_BANDS*BITALLOC_SIZE;i++)
282    {
283       scanf("%d ", band_allocation+i);
284    }
285 #endif
286    ALLOC_STACK;
287 #if !defined(VAR_ARRAYS) && !defined(USE_ALLOCA)
288    if (global_stack==NULL)
289       goto failure;
290 #endif 
291
292    /* The good thing here is that permutation of the arguments will automatically be invalid */
293    
294    if (Fs < 8000 || Fs > 96000)
295    {
296       if (error)
297          *error = CELT_BAD_ARG;
298       return NULL;
299    }
300    if (frame_size < 40 || frame_size > 1024 || frame_size%2!=0)
301    {
302       if (error)
303          *error = CELT_BAD_ARG;
304       return NULL;
305    }
306    
307    mode = celt_alloc(sizeof(CELTMode));
308    if (mode==NULL)
309       goto failure;
310    mode->Fs = Fs;
311
312    /* Pre/de-emphasis depends on sampling rate. The "standard" pre-emphasis
313       is defined as A(z) = 1 - 0.85*z^-1 at 48 kHz. Other rates should
314       approximate that. */
315    if(Fs < 12000) /* 8 kHz */
316    {
317       mode->preemph[0] =  QCONST16(.35f, 15);
318       mode->preemph[1] = -QCONST16(.18f, 15);
319       mode->preemph[2] =  QCONST16(.272f, SIG_SHIFT);
320       mode->preemph[3] =  QCONST16(3.6765f, 13);
321    } else if(Fs < 24000) /* 16 kHz */
322    {
323       mode->preemph[0] =  QCONST16(.6f, 15);
324       mode->preemph[1] = -QCONST16(.18f, 15);
325       mode->preemph[2] =  QCONST16(.4425f, SIG_SHIFT);
326       mode->preemph[3] =  QCONST16(2.259887f, 13);
327    } else if(Fs < 40000) /* 32 kHz */
328    {
329       mode->preemph[0] =  QCONST16(.78f, 15);
330       mode->preemph[1] = -QCONST16(.1f, 15);
331       mode->preemph[2] =  QCONST16(.75f, SIG_SHIFT);
332       mode->preemph[3] =  QCONST16(1.33333333f, 13);
333    } else /* 48 kHz */
334    {
335       mode->preemph[0] =  QCONST16(.85f, 15);
336       mode->preemph[1] =  QCONST16(.0f, 15);
337       mode->preemph[2] =  QCONST16(1.f, SIG_SHIFT);
338       mode->preemph[3] =  QCONST16(1.f, 13);
339    }
340
341    if ((celt_int32)frame_size*75 >= Fs && (frame_size%16)==0)
342    {
343      LM = 3;
344    } else if ((celt_int32)frame_size*150 >= Fs && (frame_size%8)==0)
345    {
346      LM = 2;
347    } else if ((celt_int32)frame_size*300 >= Fs && (frame_size%4)==0)
348    {
349      LM = 1;
350    } else
351    {
352      LM = 0;
353    }
354
355    mode->maxLM = LM;
356    mode->nbShortMdcts = 1<<LM;
357    mode->shortMdctSize = frame_size/mode->nbShortMdcts;
358    res = (mode->Fs+mode->shortMdctSize)/(2*mode->shortMdctSize);
359
360    mode->eBands = compute_ebands(Fs, mode->shortMdctSize, res, &mode->nbEBands);
361    if (mode->eBands==NULL)
362       goto failure;
363
364    mode->effEBands = mode->nbEBands;
365    while (mode->eBands[mode->effEBands] > mode->shortMdctSize)
366       mode->effEBands--;
367    
368    /* Overlap must be divisible by 4 */
369    mode->overlap = ((mode->shortMdctSize>>2)<<2);
370
371    compute_allocation_table(mode, res);
372    if (mode->allocVectors==NULL)
373       goto failure;
374    
375    window = (celt_word16*)celt_alloc(mode->overlap*sizeof(celt_word16));
376    if (window==NULL)
377       goto failure;
378
379 #ifndef FIXED_POINT
380    for (i=0;i<mode->overlap;i++)
381       window[i] = Q15ONE*sin(.5*M_PI* sin(.5*M_PI*(i+.5)/mode->overlap) * sin(.5*M_PI*(i+.5)/mode->overlap));
382 #else
383    for (i=0;i<mode->overlap;i++)
384       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))));
385 #endif
386    mode->window = window;
387
388    logN = (celt_int16*)celt_alloc(mode->nbEBands*sizeof(celt_int16));
389    if (logN==NULL)
390       goto failure;
391
392    for (i=0;i<mode->nbEBands;i++)
393       logN[i] = log2_frac(mode->eBands[i+1]-mode->eBands[i], BITRES);
394    mode->logN = logN;
395
396    compute_pulse_cache(mode, mode->maxLM);
397
398    clt_mdct_init(&mode->mdct, 2*mode->shortMdctSize*mode->nbShortMdcts, mode->maxLM);
399    if ((mode->mdct.trig==NULL)
400 #ifndef ENABLE_TI_DSPLIB55
401          || (mode->mdct.kfft==NULL)
402 #endif
403    )
404       goto failure;
405
406    if (error)
407       *error = CELT_OK;
408
409    return mode;
410 failure: 
411    if (error)
412       *error = CELT_INVALID_MODE;
413    if (mode!=NULL)
414       celt_mode_destroy(mode);
415    return NULL;
416 #endif /* !STATIC_MODES */
417 }
418
419 void celt_mode_destroy(CELTMode *mode)
420 {
421 #ifndef STATIC_MODES
422    if (mode == NULL)
423       return;
424
425    celt_free((celt_int16*)mode->eBands);
426    celt_free((celt_int16*)mode->allocVectors);
427    
428    celt_free((celt_word16*)mode->window);
429    celt_free((celt_int16*)mode->logN);
430
431    celt_free((celt_int16*)mode->cache.index);
432    celt_free((unsigned char*)mode->cache.bits);
433    clt_mdct_clear(&mode->mdct);
434
435    celt_free((CELTMode *)mode);
436 #endif
437 }