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