License for the kiss-fft headers
[opus.git] / libcelt / modes.c
1 /* (C) 2007-2008 Jean-Marc Valin, CSIRO
2    (C) 2008 Gregory Maxwell */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7    
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10    
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14    
15    - Neither the name of the Xiph.org Foundation nor the names of its
16    contributors may be used to endorse or promote products derived from
17    this software without specific prior written permission.
18    
19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32 #ifdef HAVE_CONFIG_H
33 #include "config.h"
34 #endif
35
36 #include "celt.h"
37 #include "modes.h"
38 #include "rate.h"
39 #include "os_support.h"
40 #include "stack_alloc.h"
41 #include "quant_bands.h"
42
43 #ifdef STATIC_MODES
44 #include "static_modes.c"
45 #endif
46
47 #define MODEVALID 0xa110ca7e
48 #define MODEFREED 0xb10cf8ee
49
50 #ifndef M_PI
51 #define M_PI 3.141592653
52 #endif
53
54
55 int celt_mode_info(const CELTMode *mode, int request, celt_int32_t *value)
56 {
57    switch (request)
58    {
59       case CELT_GET_FRAME_SIZE:
60          *value = mode->mdctSize;
61          break;
62       case CELT_GET_LOOKAHEAD:
63          *value = mode->overlap;
64          break;
65       case CELT_GET_NB_CHANNELS:
66          *value = mode->nbChannels;
67          break;
68       case CELT_GET_BITSTREAM_VERSION:
69          *value = CELT_BITSTREAM_VERSION;
70          break;
71       default:
72          return CELT_UNIMPLEMENTED;
73    }
74    return CELT_OK;
75 }
76
77 #ifndef STATIC_MODES
78
79 #define PBANDS 8
80
81 #ifdef STDIN_TUNING
82 int MIN_BINS;
83 #else
84 #define MIN_BINS 3
85 #endif
86
87 /* Defining 25 critical bands for the full 0-20 kHz audio bandwidth
88    Taken from http://ccrma.stanford.edu/~jos/bbt/Bark_Frequency_Scale.html */
89 #define BARK_BANDS 25
90 static const celt_int16_t bark_freq[BARK_BANDS+1] = {
91       0,   100,   200,   300,   400,
92     510,   630,   770,   920,  1080,
93    1270,  1480,  1720,  2000,  2320,
94    2700,  3150,  3700,  4400,  5300,
95    6400,  7700,  9500, 12000, 15500,
96   20000};
97
98 static const celt_int16_t pitch_freq[PBANDS+1] ={0, 345, 689, 1034, 1378, 2067, 3273, 5340, 6374};
99
100 /* This allocation table is per critical band. When creating a mode, the bits get added together 
101    into the codec bands, which are sometimes larger than one critical band at low frequency */
102
103 #ifdef STDIN_TUNING
104 int BITALLOC_SIZE;
105 int *band_allocation;
106 #else
107 #define BITALLOC_SIZE 12
108 static const int band_allocation[BARK_BANDS*BITALLOC_SIZE] = 
109    {  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,
110       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,
111       2,  2,  2,  1,  2,  2,  2,  2,  2,  1,  2,  2,  4,  5,  7,  7,  7,  5,  4,  0,  0,  0,  0,  0,  0,
112       2,  2,  2,  2,  3,  2,  2,  2,  2,  2,  3,  3,  5,  6,  8,  8,  8,  6,  5,  4,  0,  0,  0,  0,  0,
113       3,  2,  2,  2,  3,  3,  2,  3,  2,  3,  4,  4,  6,  7,  9,  9,  9,  7,  6,  5,  5,  5,  0,  0,  0,
114       3,  3,  2,  2,  3,  3,  3,  3,  3,  4,  4,  5,  7,  9, 10, 10, 10,  9,  6,  5,  5,  5,  5,  1,  0,
115       4,  3,  3,  3,  3,  3,  3,  3,  4,  4,  6,  7,  7,  9, 11, 10, 10,  9,  9,  8, 11, 10, 10,  1,  1,
116       5,  5,  5,  5,  5,  5,  5,  6,  6,  6,  8,  8, 10, 12, 12, 11, 11, 17, 12, 15, 15, 20, 18, 10,  1,
117       8,  7,  7,  7,  7,  7,  8,  8,  9, 10, 11, 12, 14, 17, 18, 21, 22, 27, 29, 39, 37, 38, 40, 35,  1,
118       7,  7,  7,  7,  7,  7, 10, 10, 10, 13, 14, 18, 20, 24, 28, 32, 32, 35, 38, 38, 42, 50, 59, 54, 31,
119       8,  8,  8,  8,  8,  9, 10, 12, 14, 20, 22, 25, 28, 30, 35, 42, 46, 50, 55, 60, 62, 62, 62, 62, 62,
120      12, 12, 12, 12, 12, 13, 15, 18, 22, 30, 32, 35, 40, 45, 55, 62, 66, 70, 85, 90, 92, 92, 92, 92, 92,
121    };
122 #endif
123
124 static celt_int16_t *compute_ebands(celt_int32_t Fs, int frame_size, int *nbEBands)
125 {
126    celt_int16_t *eBands;
127    int i, res, min_width, lin, low, high, nBark;
128    res = (Fs+frame_size)/(2*frame_size);
129    min_width = MIN_BINS*res;
130    /*printf ("min_width = %d\n", min_width);*/
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] >= min_width)
140          break;
141    
142    /*printf ("lin = %d (%d Hz)\n", lin, bark_freq[lin]);*/
143    low = ((bark_freq[lin]/res)+(MIN_BINS-1))/MIN_BINS;
144    high = nBark-lin;
145    *nbEBands = low+high;
146    eBands = celt_alloc(sizeof(celt_int16_t)*(*nbEBands+2));
147    
148    /* Linear spacing (min_width) */
149    for (i=0;i<low;i++)
150       eBands[i] = MIN_BINS*i;
151    /* Spacing follows critical bands */
152    for (i=0;i<high;i++)
153       eBands[i+low] = (bark_freq[lin+i]+res/2)/res;
154    /* Enforce the minimum spacing at the boundary */
155    for (i=0;i<*nbEBands;i++)
156       if (eBands[i] < MIN_BINS*i)
157          eBands[i] = MIN_BINS*i;
158    eBands[*nbEBands] = (bark_freq[nBark]+res/2)/res;
159    eBands[*nbEBands+1] = frame_size;
160    if (eBands[*nbEBands] > eBands[*nbEBands+1])
161       eBands[*nbEBands] = eBands[*nbEBands+1];
162    
163    /* FIXME: Remove last band if too small */
164    /*for (i=0;i<*nbEBands+2;i++)
165       printf("%d ", eBands[i]);
166    printf ("\n");
167    exit(1);*/
168    return eBands;
169 }
170
171 static void compute_pbands(CELTMode *mode, int res)
172 {
173    int i;
174    celt_int16_t *pBands;
175    pBands=celt_alloc(sizeof(celt_int16_t)*(PBANDS+2));
176    mode->nbPBands = PBANDS;
177    for (i=0;i<PBANDS+1;i++)
178    {
179       pBands[i] = (pitch_freq[i]+res/2)/res;
180       if (pBands[i] < mode->eBands[i])
181          pBands[i] = mode->eBands[i];
182    }
183    pBands[PBANDS+1] = mode->eBands[mode->nbEBands+1];
184    for (i=1;i<mode->nbPBands+1;i++)
185    {
186       int j;
187       for (j=0;j<mode->nbEBands;j++)
188          if (mode->eBands[j] <= pBands[i] && mode->eBands[j+1] > pBands[i])
189             break;
190       /*printf ("%d %d\n", i, j);*/
191       if (mode->eBands[j] != pBands[i])
192       {
193          if (pBands[i]-mode->eBands[j] < mode->eBands[j+1]-pBands[i] && 
194              mode->eBands[j] != pBands[i-1])
195             pBands[i] = mode->eBands[j];
196          else
197             pBands[i] = mode->eBands[j+1];
198       }
199    }
200    /*for (i=0;i<mode->nbPBands+2;i++)
201       printf("%d ", pBands[i]);
202    printf ("\n");*/
203    mode->pBands = pBands;
204    mode->pitchEnd = pBands[PBANDS];
205 }
206
207 static void compute_allocation_table(CELTMode *mode, int res)
208 {
209    int i, j, eband, nBark;
210    celt_int16_t *allocVectors, *allocEnergy;
211    const int C = CHANNELS(mode);
212
213    /* Find the number of critical bands supported by our sampling rate */
214    for (nBark=1;nBark<BARK_BANDS;nBark++)
215     if (bark_freq[nBark+1]*2 >= mode->Fs)
216        break;
217
218    mode->nbAllocVectors = BITALLOC_SIZE;
219    allocVectors = celt_alloc(sizeof(celt_int16_t)*(BITALLOC_SIZE*mode->nbEBands));
220    allocEnergy = celt_alloc(sizeof(celt_int16_t)*(mode->nbAllocVectors*(mode->nbEBands+1)));
221    /* Compute per-codec-band allocation from per-critical-band matrix */
222    for (i=0;i<BITALLOC_SIZE;i++)
223    {
224       eband = 0;
225       for (j=0;j<nBark;j++)
226       {
227          int edge, low;
228          celt_int32_t alloc;
229          edge = mode->eBands[eband+1]*res;
230          alloc = band_allocation[i*BARK_BANDS+j];
231          alloc = alloc*C*mode->mdctSize/256;
232          if (edge < bark_freq[j+1])
233          {
234             int num, den;
235             num = alloc * (edge-bark_freq[j]);
236             den = bark_freq[j+1]-bark_freq[j];
237             low = (num+den/2)/den;
238             allocVectors[i*mode->nbEBands+eband] += low;
239             eband++;
240             allocVectors[i*mode->nbEBands+eband] += alloc-low;
241          } else {
242             allocVectors[i*mode->nbEBands+eband] += alloc;
243          }
244       }
245    }
246    /* Compute fine energy resolution and update the pulse allocation table to subtract that */
247    for (i=0;i<mode->nbAllocVectors;i++)
248    {
249       int sum = 0;
250       for (j=0;j<mode->nbEBands;j++)
251       {
252          int ebits;
253          int min_bits=0;
254          if (allocVectors[i*mode->nbEBands+j] > 0)
255             min_bits = 1;
256          ebits = IMAX(min_bits , allocVectors[i*mode->nbEBands+j] / (C*(mode->eBands[j+1]-mode->eBands[j])));
257          if (ebits>7)
258             ebits=7;
259          /* The bits used for fine allocation can't be used for pulses */
260          /* However, we give two "free" bits to all modes to compensate for the fact that some energy
261             resolution is needed regardless of the frame size. */
262          if (ebits>1)
263             allocVectors[i*mode->nbEBands+j] -= C*(ebits-2);
264          if (allocVectors[i*mode->nbEBands+j] < 0)
265             allocVectors[i*mode->nbEBands+j] = 0;
266          sum += ebits;
267          allocEnergy[i*(mode->nbEBands+1)+j] = ebits;
268       }
269       allocEnergy[i*(mode->nbEBands+1)+mode->nbEBands] = sum;
270    }
271    mode->energy_alloc = allocEnergy;
272    mode->allocVectors = allocVectors;
273 }
274
275 #endif /* STATIC_MODES */
276
277 CELTMode *celt_mode_create(celt_int32_t Fs, int channels, int frame_size, int *error)
278 {
279    int i;
280 #ifdef STDIN_TUNING
281    scanf("%d ", &MIN_BINS);
282    scanf("%d ", &BITALLOC_SIZE);
283    band_allocation = celt_alloc(sizeof(int)*BARK_BANDS*BITALLOC_SIZE);
284    for (i=0;i<BARK_BANDS*BITALLOC_SIZE;i++)
285    {
286       scanf("%d ", band_allocation+i);
287    }
288 #endif
289 #ifdef STATIC_MODES
290    const CELTMode *m = NULL;
291    CELTMode *mode=NULL;
292    ALLOC_STACK;
293    for (i=0;i<TOTAL_MODES;i++)
294    {
295       if (Fs == static_mode_list[i]->Fs &&
296           channels == static_mode_list[i]->nbChannels &&
297           frame_size == static_mode_list[i]->mdctSize)
298       {
299          m = static_mode_list[i];
300          break;
301       }
302    }
303    if (m == NULL)
304    {
305       celt_warning("Mode not included as part of the static modes");
306       if (error)
307          *error = CELT_BAD_ARG;
308       return NULL;
309    }
310    mode = (CELTMode*)celt_alloc(sizeof(CELTMode));
311    CELT_COPY(mode, m, 1);
312 #else
313    int res;
314    CELTMode *mode;
315    celt_word16_t *window;
316    ALLOC_STACK;
317
318    /* The good thing here is that permutation of the arguments will automatically be invalid */
319    
320    if (Fs < 32000 || Fs > 96000)
321    {
322       celt_warning("Sampling rate must be between 32 kHz and 96 kHz");
323       if (error)
324          *error = CELT_BAD_ARG;
325       return NULL;
326    }
327    if (channels < 0 || channels > 2)
328    {
329       celt_warning("Only mono and stereo supported");
330       if (error)
331          *error = CELT_BAD_ARG;
332       return NULL;
333    }
334    if (frame_size < 64 || frame_size > 512 || frame_size%2!=0)
335    {
336       celt_warning("Only even frame sizes from 64 to 512 are supported");
337       if (error)
338          *error = CELT_BAD_ARG;
339       return NULL;
340    }
341    res = (Fs+frame_size)/(2*frame_size);
342    
343    mode = celt_alloc(sizeof(CELTMode));
344    mode->Fs = Fs;
345    mode->mdctSize = frame_size;
346    mode->nbChannels = channels;
347    mode->eBands = compute_ebands(Fs, frame_size, &mode->nbEBands);
348    compute_pbands(mode, res);
349    mode->ePredCoef = QCONST16(.8f,15);
350
351    if (frame_size > 384 && (frame_size%8)==0)
352    {
353      mode->nbShortMdcts = 4;
354    } else if (frame_size > 384 && (frame_size%10)==0)
355    {
356      mode->nbShortMdcts = 5;
357    } else if (frame_size > 256 && (frame_size%6)==0)
358    {
359      mode->nbShortMdcts = 3;
360    } else if (frame_size > 256 && (frame_size%8)==0)
361    {
362      mode->nbShortMdcts = 4;
363    } else if (frame_size > 64 && (frame_size%4)==0)
364    {
365      mode->nbShortMdcts = 2;
366    } else if (frame_size > 128 && (frame_size%6)==0)
367    {
368      mode->nbShortMdcts = 3;
369    } else
370    {
371      mode->nbShortMdcts = 1;
372    }
373
374    if (mode->nbShortMdcts > 1)
375       mode->overlap = ((frame_size/mode->nbShortMdcts)>>2)<<2; /* Overlap must be divisible by 4 */
376    else
377       mode->overlap = (frame_size>>3)<<2;
378
379    compute_allocation_table(mode, res);
380    /*printf ("%d bands\n", mode->nbEBands);*/
381    
382    window = (celt_word16_t*)celt_alloc(mode->overlap*sizeof(celt_word16_t));
383
384 #ifndef FIXED_POINT
385    for (i=0;i<mode->overlap;i++)
386       window[i] = Q15ONE*sin(.5*M_PI* sin(.5*M_PI*(i+.5)/mode->overlap) * sin(.5*M_PI*(i+.5)/mode->overlap));
387 #else
388    for (i=0;i<mode->overlap;i++)
389       window[i] = MIN32(32767,32768.*sin(.5*M_PI* sin(.5*M_PI*(i+.5)/mode->overlap) * sin(.5*M_PI*(i+.5)/mode->overlap)));
390 #endif
391    mode->window = window;
392
393    mode->bits = (const celt_int16_t **)compute_alloc_cache(mode, 1);
394    if (mode->nbChannels>=2)
395       mode->bits_stereo = (const celt_int16_t **)compute_alloc_cache(mode, mode->nbChannels);
396
397 #ifndef SHORTCUTS
398    psydecay_init(&mode->psy, MAX_PERIOD/2, mode->Fs);
399 #endif
400    
401    mode->marker_start = MODEVALID;
402    mode->marker_end = MODEVALID;
403 #endif /* !STATIC_MODES */
404    mdct_init(&mode->mdct, 2*mode->mdctSize);
405    mode->fft = pitch_state_alloc(MAX_PERIOD);
406
407    mode->shortMdctSize = mode->mdctSize/mode->nbShortMdcts;
408    mdct_init(&mode->shortMdct, 2*mode->shortMdctSize);
409    mode->shortWindow = mode->window;
410
411    mode->prob = quant_prob_alloc(mode);
412
413    if (error)
414       *error = CELT_OK;
415    return mode;
416 }
417
418 void celt_mode_destroy(CELTMode *mode)
419 {
420 #ifndef STATIC_MODES
421    int i;
422    const celt_int16_t *prevPtr = NULL;
423    for (i=0;i<mode->nbEBands;i++)
424    {
425       if (mode->bits[i] != prevPtr)
426       {
427          prevPtr = mode->bits[i];
428          celt_free((int*)mode->bits[i]);
429       }
430    }
431    celt_free((int**)mode->bits);
432    if (mode->bits_stereo != NULL)
433    {
434       for (i=0;i<mode->nbEBands;i++)
435       {
436          if (mode->bits_stereo[i] != prevPtr)
437          {
438             prevPtr = mode->bits_stereo[i];
439             celt_free((int*)mode->bits_stereo[i]);
440          }
441       }
442       celt_free((int**)mode->bits_stereo);
443    }
444    if (check_mode(mode) != CELT_OK)
445       return;
446    celt_free((int*)mode->eBands);
447    celt_free((int*)mode->pBands);
448    celt_free((int*)mode->allocVectors);
449    celt_free((celt_int16_t *)mode->energy_alloc);
450
451    celt_free((celt_word16_t*)mode->window);
452
453    mode->marker_start = MODEFREED;
454    mode->marker_end = MODEFREED;
455 #ifndef SHORTCUTS
456    psydecay_clear(&mode->psy);
457 #endif
458 #endif
459    mdct_clear(&mode->mdct);
460    mdct_clear(&mode->shortMdct);
461    pitch_state_free(mode->fft);
462    quant_prob_free(mode->prob);
463    celt_free((CELTMode *)mode);
464 }
465
466 int check_mode(const CELTMode *mode)
467 {
468    if (mode->marker_start == MODEVALID && mode->marker_end == MODEVALID)
469       return CELT_OK;
470    if (mode->marker_start == MODEFREED || mode->marker_end == MODEFREED)
471       celt_warning("Using a mode that has already been freed");
472    else
473       celt_warning("This is not a valid CELT mode");
474    return CELT_INVALID_MODE;
475 }