Making the allocation matrix a bit smoother
[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, 56, 44, 40, 32, 24, 20, 40, 48, 52, 44, 32, 16,  8,  4,  0,  0,  0,  0,  0,
100  52, 40, 68, 64, 56, 48, 40, 32, 48, 56, 56, 48, 36, 20, 12,  8,  8,  4,  0,  0,  0,
101  68, 84, 92,104, 96, 80, 68, 64, 68, 72, 64, 56, 44, 24, 12,  8,  8,  4,  4,  0,  0,
102 100,100,110,110,100, 96, 90, 88, 84, 76, 68, 60, 44, 28, 20, 20, 20, 12, 12,  0,  0,
103 130,130,130,128,120,112,104,104,100, 96, 76, 68, 60, 52, 36, 36, 32, 28, 20,  8,  0,
104 168,168,168,148,140,136,135,130,136,140,128,124,108, 96, 88, 83, 72, 56, 44, 28,  2,
105 184,196,184,184,168,172,176,188,200,208,204,192,156,128,108, 96, 88, 76, 68, 44, 20,
106 212,212,196,192,220,255,255,255,255,255,255,255,224,188,164,148,124, 96, 80, 64, 40,
107 255,255,255,255,255,255,255,255,255,255,255,255,255,255,252,220,188,144,104, 84, 60,
108 };
109
110 static celt_int16 *compute_ebands(celt_int32 Fs, int frame_size, int res, int *nbEBands)
111 {
112    celt_int16 *eBands;
113    int i, lin, low, high, nBark, offset=0;
114
115    /* All modes that have 2.5 ms short blocks use the same definition */
116    if (Fs == 400*(celt_int32)frame_size)
117    {
118       *nbEBands = sizeof(eband5ms)/sizeof(eband5ms[0])-1;
119       eBands = celt_alloc(sizeof(celt_int16)*(*nbEBands+1));
120       for (i=0;i<*nbEBands+1;i++)
121          eBands[i] = eband5ms[i];
122       return eBands;
123    }
124    /* Find the number of critical bands supported by our sampling rate */
125    for (nBark=1;nBark<BARK_BANDS;nBark++)
126     if (bark_freq[nBark+1]*2 >= Fs)
127        break;
128
129    /* Find where the linear part ends (i.e. where the spacing is more than min_width */
130    for (lin=0;lin<nBark;lin++)
131       if (bark_freq[lin+1]-bark_freq[lin] >= res)
132          break;
133
134    low = (bark_freq[lin]+res/2)/res;
135    high = nBark-lin;
136    *nbEBands = low+high;
137    eBands = celt_alloc(sizeof(celt_int16)*(*nbEBands+2));
138    
139    if (eBands==NULL)
140       return NULL;
141    
142    /* Linear spacing (min_width) */
143    for (i=0;i<low;i++)
144       eBands[i] = i;
145    if (low>0)
146       offset = eBands[low-1]*res - bark_freq[lin-1];
147    /* Spacing follows critical bands */
148    for (i=0;i<high;i++)
149    {
150       int target = bark_freq[lin+i];
151       eBands[i+low] = (target+(offset+res)/2)/res;
152       offset = eBands[i+low]*res - target;
153    }
154    /* Enforce the minimum spacing at the boundary */
155    for (i=0;i<*nbEBands;i++)
156       if (eBands[i] < i)
157          eBands[i] = i;
158    eBands[*nbEBands] = (bark_freq[nBark]+res/2)/res;
159    if (eBands[*nbEBands] > frame_size)
160       eBands[*nbEBands] = frame_size;
161    for (i=1;i<*nbEBands-1;i++)
162    {
163       if (eBands[i+1]-eBands[i] < eBands[i]-eBands[i-1])
164       {
165          eBands[i] -= (2*eBands[i]-eBands[i-1]-eBands[i+1])/2;
166       }
167    }
168    /*for (i=0;i<=*nbEBands+1;i++)
169       printf ("%d ", eBands[i]);
170    printf ("\n");
171    exit(1);*/
172    /* FIXME: Remove last band if too small */
173    return eBands;
174 }
175
176 static void compute_allocation_table(CELTMode *mode, int res)
177 {
178    int i, j;
179    unsigned char *allocVectors;
180    int maxBands = sizeof(eband5ms)/sizeof(eband5ms[0])-1;
181
182    mode->nbAllocVectors = BITALLOC_SIZE;
183    allocVectors = celt_alloc(sizeof(unsigned char)*(BITALLOC_SIZE*mode->nbEBands));
184    if (allocVectors==NULL)
185       return;
186
187    /* Check for standard mode */
188    if (mode->Fs == 400*(celt_int32)mode->shortMdctSize && mode->Fs >= 40000)
189    {
190       for (i=0;i<BITALLOC_SIZE*mode->nbEBands;i++)
191          allocVectors[i] = band_allocation[i];
192       mode->allocVectors = allocVectors;
193       return;
194    }
195    /* If not the standard mode, interpolate */
196
197    /* Compute per-codec-band allocation from per-critical-band matrix */
198    for (i=0;i<BITALLOC_SIZE;i++)
199    {
200       celt_int32 current = 0;
201       int eband = 0;
202       /* We may be looping over too many bands, but eband will stop being
203          incremented once we reach the last band */
204       for (j=0;j<maxBands;j++)
205       {
206          int edge, low, high;
207          celt_int32 alloc;
208          alloc = band_allocation[i*maxBands + j]*(mode->eBands[eband+1]-mode->eBands[eband])<<4;
209          low = eband5ms[j]*200;
210          high = eband5ms[j+1]*200;
211          edge = mode->eBands[eband+1]*res;
212          while (edge <= high && eband < mode->nbEBands)
213          {
214             celt_int32 num;
215             int den, bits;
216             int N = (mode->eBands[eband+1]-mode->eBands[eband]);
217             num = alloc * (edge-low);
218             den = high-low;
219             /* Divide with rounding */
220             bits = (2*num+den)/(2*den);
221             allocVectors[i*mode->nbEBands+eband] = (2*(current+bits)+(N<<4))/(2*N<<4);
222             /* Remove the part of the band we just allocated */
223             low = edge;
224             alloc -= bits;
225
226             /* Move to next eband */
227             current = 0;
228             eband++;
229             if (eband < mode->nbEBands)
230                edge = mode->eBands[eband+1]*res;
231          }
232          current += alloc;
233       }
234       if (eband < mode->nbEBands)
235       {
236          int N = (mode->eBands[eband+1]-mode->eBands[eband]);
237          allocVectors[i*mode->nbEBands+eband] = (2*current+(N<<4))/(2*N<<4);
238       }
239    }
240    /*printf ("\n");
241    for (i=0;i<BITALLOC_SIZE;i++)
242    {
243       for (j=0;j<mode->nbEBands;j++)
244          printf ("%d ", allocVectors[i*mode->nbEBands+j]);
245       printf ("\n");
246    }
247    exit(0);*/
248
249    mode->allocVectors = allocVectors;
250 }
251
252 #endif /* STATIC_MODES */
253
254 CELTMode *celt_mode_create(celt_int32 Fs, int frame_size, int *error)
255 {
256    int i;
257 #ifdef STATIC_MODES
258    for (i=0;i<TOTAL_MODES;i++)
259    {
260       if (Fs == static_mode_list[i]->Fs &&
261           frame_size == static_mode_list[i]->shortMdctSize*static_mode_list[i]->nbShortMdcts)
262       {
263          return (CELTMode*)static_mode_list[i];
264       }
265    }
266    if (error)
267       *error = CELT_BAD_ARG;
268    return NULL;
269 #else
270    int res;
271    CELTMode *mode=NULL;
272    celt_word16 *window;
273    celt_int16 *logN;
274    int LM;
275 #ifdef STDIN_TUNING
276    scanf("%d ", &MIN_BINS);
277    scanf("%d ", &BITALLOC_SIZE);
278    band_allocation = celt_alloc(sizeof(int)*BARK_BANDS*BITALLOC_SIZE);
279    for (i=0;i<BARK_BANDS*BITALLOC_SIZE;i++)
280    {
281       scanf("%d ", band_allocation+i);
282    }
283 #endif
284    ALLOC_STACK;
285 #if !defined(VAR_ARRAYS) && !defined(USE_ALLOCA)
286    if (global_stack==NULL)
287       goto failure;
288 #endif 
289
290    /* The good thing here is that permutation of the arguments will automatically be invalid */
291    
292    if (Fs < 8000 || Fs > 96000)
293    {
294       if (error)
295          *error = CELT_BAD_ARG;
296       return NULL;
297    }
298    if (frame_size < 40 || frame_size > 1024 || frame_size%2!=0)
299    {
300       if (error)
301          *error = CELT_BAD_ARG;
302       return NULL;
303    }
304    
305    mode = celt_alloc(sizeof(CELTMode));
306    if (mode==NULL)
307       goto failure;
308    mode->Fs = Fs;
309
310    /* Pre/de-emphasis depends on sampling rate. The "standard" pre-emphasis
311       is defined as A(z) = 1 - 0.85*z^-1 at 48 kHz. Other rates should
312       approximate that. */
313    if(Fs < 12000) /* 8 kHz */
314    {
315       mode->preemph[0] =  QCONST16(.35f, 15);
316       mode->preemph[1] = -QCONST16(.18f, 15);
317       mode->preemph[2] =  QCONST16(.272f, SIG_SHIFT);
318       mode->preemph[3] =  QCONST16(3.6765f, 13);
319    } else if(Fs < 24000) /* 16 kHz */
320    {
321       mode->preemph[0] =  QCONST16(.6f, 15);
322       mode->preemph[1] = -QCONST16(.18f, 15);
323       mode->preemph[2] =  QCONST16(.4425f, SIG_SHIFT);
324       mode->preemph[3] =  QCONST16(2.259887f, 13);
325    } else if(Fs < 40000) /* 32 kHz */
326    {
327       mode->preemph[0] =  QCONST16(.78f, 15);
328       mode->preemph[1] = -QCONST16(.1f, 15);
329       mode->preemph[2] =  QCONST16(.75f, SIG_SHIFT);
330       mode->preemph[3] =  QCONST16(1.33333333f, 13);
331    } else /* 48 kHz */
332    {
333       mode->preemph[0] =  QCONST16(.85f, 15);
334       mode->preemph[1] =  QCONST16(.0f, 15);
335       mode->preemph[2] =  QCONST16(1.f, SIG_SHIFT);
336       mode->preemph[3] =  QCONST16(1.f, 13);
337    }
338
339    if ((celt_int32)frame_size*75 >= Fs && (frame_size%16)==0)
340    {
341      LM = 3;
342    } else if ((celt_int32)frame_size*150 >= Fs && (frame_size%8)==0)
343    {
344      LM = 2;
345    } else if ((celt_int32)frame_size*300 >= Fs && (frame_size%4)==0)
346    {
347      LM = 1;
348    } else
349    {
350      LM = 0;
351    }
352
353    mode->maxLM = LM;
354    mode->nbShortMdcts = 1<<LM;
355    mode->shortMdctSize = frame_size/mode->nbShortMdcts;
356    res = (mode->Fs+mode->shortMdctSize)/(2*mode->shortMdctSize);
357
358    mode->eBands = compute_ebands(Fs, mode->shortMdctSize, res, &mode->nbEBands);
359    if (mode->eBands==NULL)
360       goto failure;
361
362    mode->effEBands = mode->nbEBands;
363    while (mode->eBands[mode->effEBands] > mode->shortMdctSize)
364       mode->effEBands--;
365    
366    /* Overlap must be divisible by 4 */
367    if (mode->nbShortMdcts > 1)
368       mode->overlap = (mode->shortMdctSize>>2)<<2;
369    else
370       mode->overlap = (frame_size>>3)<<2;
371
372
373    compute_allocation_table(mode, res);
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    mode->prob = quant_prob_alloc(mode);
409    if (mode->prob==NULL)
410      goto failure;
411
412    if (error)
413       *error = CELT_OK;
414
415    return mode;
416 failure: 
417    if (error)
418       *error = CELT_INVALID_MODE;
419    if (mode!=NULL)
420       celt_mode_destroy(mode);
421    return NULL;
422 #endif /* !STATIC_MODES */
423 }
424
425 void celt_mode_destroy(CELTMode *mode)
426 {
427 #ifndef STATIC_MODES
428    if (mode == NULL)
429       return;
430
431    celt_free((celt_int16*)mode->eBands);
432    celt_free((celt_int16*)mode->allocVectors);
433    
434    celt_free((celt_word16*)mode->window);
435    celt_free((celt_int16*)mode->logN);
436
437    celt_free((celt_int16*)mode->cache.index);
438    celt_free((unsigned char*)mode->cache.bits);
439    clt_mdct_clear(&mode->mdct);
440    quant_prob_free(mode->prob);
441
442    celt_free((CELTMode *)mode);
443 #endif
444 }