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