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