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