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