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