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