celt.c 75.8 KB
Newer Older
Jean-Marc Valin's avatar
Jean-Marc Valin committed
1
/* Copyright (c) 2007-2008 CSIRO
Jean-Marc Valin's avatar
Jean-Marc Valin committed
2
   Copyright (c) 2007-2010 Xiph.Org Foundation
Jean-Marc Valin's avatar
Jean-Marc Valin committed
3
4
   Copyright (c) 2008 Gregory Maxwell 
   Written by Jean-Marc Valin and Gregory Maxwell */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
/*
   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:
   
   - Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
   
   - Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
   
   - Neither the name of the Xiph.org Foundation nor the names of its
   contributors may be used to endorse or promote products derived from
   this software without specific prior written permission.
   
   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

34
35
36
37
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

38
39
#define CELT_C

40
41
42
43
44
/* Always enable postfilter for Opus */
#if defined(OPUS_BUILD) && !defined(ENABLE_POSTFILTER)
#define ENABLE_POSTFILTER
#endif

Jean-Marc Valin's avatar
Jean-Marc Valin committed
45
#include "os_support.h"
46
47
#include "mdct.h"
#include <math.h>
48
#include "celt.h"
49
#include "pitch.h"
50
#include "bands.h"
51
#include "modes.h"
52
#include "entcode.h"
53
#include "quant_bands.h"
Jean-Marc Valin's avatar
Jean-Marc Valin committed
54
#include "rate.h"
55
#include "stack_alloc.h"
56
#include "mathops.h"
57
#include "float_cast.h"
58
#include <stdarg.h>
59
#include "plc.h"
60
#include "vq.h"
61

62
63
64
static const unsigned char trim_icdf[11] = {126, 124, 119, 109, 87, 41, 19, 9, 4, 2, 0};
/* Probs: NONE: 21.875%, LIGHT: 6.25%, NORMAL: 65.625%, AGGRESSIVE: 6.25% */
static const unsigned char spread_icdf[4] = {25, 23, 2, 0};
65

66
67
static const unsigned char tapset_icdf[3]={2,1,0};

68
#define COMBFILTER_MAXPERIOD 1024
Jean-Marc Valin's avatar
Jean-Marc Valin committed
69
#define COMBFILTER_MINPERIOD 15
70

Jean-Marc Valin's avatar
Jean-Marc Valin committed
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
static int resampling_factor(celt_int32 rate)
{
   int ret;
   switch (rate)
   {
   case 48000:
      ret = 1;
      break;
   case 24000:
      ret = 2;
      break;
   case 16000:
      ret = 3;
      break;
   case 12000:
      ret = 4;
      break;
   case 8000:
      ret = 6;
      break;
   default:
      ret = 0;
   }
   return ret;
}

Jean-Marc Valin's avatar
Jean-Marc Valin committed
97
98
99
/** Encoder state 
 @brief Encoder state
 */
100
struct CELTEncoder {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
101
   const CELTMode *mode;     /**< Mode used by the encoder */
102
   int overlap;
103
   int channels;
104
   int stream_channels;
105
   
106
   int force_intra;
107
   int disable_pf;
108
   int complexity;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
109
   int upsample;
110
111
   int start, end;

112
113
   celt_int32 bitrate;
   int vbr;
114
   int constrained_vbr;      /* If zero, VBR can do whatever it likes with the rate */
115
116

   /* Everything beyond this point gets cleared on a reset */
117
#define ENCODER_RESET_START rng
118

119
   ec_uint32 rng;
120
   int spread_decision;
121
   int delayedIntra;
122
   int tonal_average;
123
   int lastCodedBands;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
124
125
   int hf_average;
   int tapset_decision;
126

127
128
   int prefilter_period;
   celt_word16 prefilter_gain;
129
   int prefilter_tapset;
130
131
132
#ifdef RESYNTH
   int prefilter_period_old;
   celt_word16 prefilter_gain_old;
133
   int prefilter_tapset_old;
134
#endif
135
   int consec_transient;
136

137
138
139
140
141
142
   /* VBR-related parameters */
   celt_int32 vbr_reservoir;
   celt_int32 vbr_drift;
   celt_int32 vbr_offset;
   celt_int32 vbr_count;

Jean-Marc Valin's avatar
Jean-Marc Valin committed
143
144
   celt_word32 preemph_memE[2];
   celt_word32 preemph_memD[2];
145

146
147
148
149
#ifdef RESYNTH
   celt_sig syn_mem[2][2*MAX_PERIOD];
#endif

150
   celt_sig in_mem[1]; /* Size = channels*mode->overlap */
151
   /* celt_sig prefilter_mem[],  Size = channels*COMBFILTER_PERIOD */
152
   /* celt_sig overlap_mem[],  Size = channels*mode->overlap */
153
   /* celt_word16 oldEBands[], Size = 2*channels*mode->nbEBands */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
154
155
};

156
157
158
159
160
161
162
int celt_encoder_get_size(int channels)
{
   CELTMode *mode = celt_mode_create(48000, 960, NULL);
   return celt_encoder_get_size_custom(mode, channels);
}

int celt_encoder_get_size_custom(const CELTMode *mode, int channels)
163
164
165
{
   int size = sizeof(struct CELTEncoder)
         + (2*channels*mode->overlap-1)*sizeof(celt_sig)
166
         + channels*COMBFILTER_MAXPERIOD*sizeof(celt_sig)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
167
         + 3*channels*mode->nbEBands*sizeof(celt_word16);
168
169
170
   return size;
}

Jean-Marc Valin's avatar
Jean-Marc Valin committed
171
CELTEncoder *celt_encoder_create(int sampling_rate, int channels, int *error)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
172
{
173
   CELTEncoder *st;
174
   st = (CELTEncoder *)celt_alloc(celt_encoder_get_size(channels));
175
176
177
178
179
   if (st!=NULL && celt_encoder_init(st, sampling_rate, channels, error)==NULL)
   {
      celt_encoder_destroy(st);
      st = NULL;
   }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
180
   return st;
181
182
183
}

CELTEncoder *celt_encoder_create_custom(const CELTMode *mode, int channels, int *error)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
184
{
185
   CELTEncoder *st = (CELTEncoder *)celt_alloc(celt_encoder_get_size_custom(mode, channels));
186
187
188
189
190
191
   if (st!=NULL && celt_encoder_init_custom(st, mode, channels, error)==NULL)
   {
      celt_encoder_destroy(st);
      st = NULL;
   }
   return st;
192
}
193

Jean-Marc Valin's avatar
Jean-Marc Valin committed
194
CELTEncoder *celt_encoder_init(CELTEncoder *st, int sampling_rate, int channels, int *error)
195
{
Jean-Marc Valin's avatar
Jean-Marc Valin committed
196
197
198
199
200
201
202
203
204
   celt_encoder_init_custom(st, celt_mode_create(48000, 960, NULL), channels, error);
   st->upsample = resampling_factor(sampling_rate);
   if (st->upsample==0)
   {
      if (error)
         *error = CELT_BAD_ARG;
      return NULL;
   }
   return st;
205
206
207
}

CELTEncoder *celt_encoder_init_custom(CELTEncoder *st, const CELTMode *mode, int channels, int *error)
208
{
209
210
211
212
213
214
   if (channels < 0 || channels > 2)
   {
      if (error)
         *error = CELT_BAD_ARG;
      return NULL;
   }
215

216
   if (st==NULL || mode==NULL)
217
218
219
   {
      if (error)
         *error = CELT_ALLOC_FAIL;
220
      return NULL;
221
   }
222

223
   CELT_MEMSET((char*)st, 0, celt_encoder_get_size_custom(mode, channels));
224
   
225
   st->mode = mode;
226
   st->overlap = mode->overlap;
227
   st->stream_channels = st->channels = channels;
228

Jean-Marc Valin's avatar
Jean-Marc Valin committed
229
   st->upsample = 1;
230
   st->start = 0;
231
   st->end = st->mode->effEBands;
232
   st->constrained_vbr = 1;
233

234
235
   st->bitrate = 255000*channels;
   st->vbr = 0;
236
   st->vbr_offset = 0;
237
   st->force_intra  = 0;
238
   st->delayedIntra = 1;
239
   st->tonal_average = 256;
240
   st->spread_decision = SPREAD_NORMAL;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
241
242
   st->hf_average = 0;
   st->tapset_decision = 0;
243
   st->complexity = 5;
244

245
   if (error)
246
247
      *error = CELT_OK;
   return st;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
248
249
}

250
void celt_encoder_destroy(CELTEncoder *st)
251
252
253
254
{
   celt_free(st);
}

255
static inline celt_int16 FLOAT2INT16(float x)
256
{
257
   x = x*CELT_SIG_SCALE;
258
259
   x = MAX32(x, -32768);
   x = MIN32(x, 32767);
260
   return (celt_int16)float2int(x);
261
262
}

263
static inline celt_word16 SIG2WORD16(celt_sig x)
264
{
265
#ifdef FIXED_POINT
266
267
268
   x = PSHR32(x, SIG_SHIFT);
   x = MAX32(x, -32768);
   x = MIN32(x, 32767);
269
270
   return EXTRACT16(x);
#else
271
   return (celt_word16)x;
272
273
274
#endif
}

275
static int transient_analysis(const celt_word32 * restrict in, int len, int C,
276
                              int overlap)
277
{
278
   int i;
279
280
   VARDECL(celt_word16, tmp);
   celt_word32 mem0=0,mem1=0;
281
282
283
284
285
   int is_transient = 0;
   int block;
   int N;
   /* FIXME: Make that smaller */
   celt_word16 bins[50];
286
   SAVE_STACK;
287
288
   ALLOC(tmp, len, celt_word16);

289
290
   block = overlap/2;
   N=len/block;
291
   if (C==1)
292
293
   {
      for (i=0;i<len;i++)
294
         tmp[i] = SHR32(in[i],SIG_SHIFT);
295
296
   } else {
      for (i=0;i<len;i++)
297
         tmp[i] = SHR32(ADD32(in[i],in[i+len]), SIG_SHIFT+1);
298
299
300
301
302
303
304
305
306
307
308
309
310
   }

   /* High-pass filter: (1 - 2*z^-1 + z^-2) / (1 - z^-1 + .5*z^-2) */
   for (i=0;i<len;i++)
   {
      celt_word32 x,y;
      x = tmp[i];
      y = ADD32(mem0, x);
#ifdef FIXED_POINT
      mem0 = mem1 + y - SHL32(x,1);
      mem1 = x - SHR32(y,1);
#else
      mem0 = mem1 + y - 2*x;
311
      mem1 = x - .5f*y;
312
313
#endif
      tmp[i] = EXTRACT16(SHR(y,2));
314
   }
315
   /* First few samples are bad because we don't propagate the memory */
316
   for (i=0;i<12;i++)
317
318
      tmp[i] = 0;

319
   for (i=0;i<N;i++)
320
   {
321
322
323
324
325
      int j;
      float max_abs=0;
      for (j=0;j<block;j++)
         max_abs = MAX32(max_abs, tmp[i*block+j]);
      bins[i] = max_abs;
326
   }
327
   for (i=0;i<N;i++)
328
   {
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
      int j;
      int conseq=0;
      celt_word16 t1, t2, t3;

      t1 = MULT16_16_Q15(QCONST16(.15f, 15), bins[i]);
      t2 = MULT16_16_Q15(QCONST16(.4f, 15), bins[i]);
      t3 = MULT16_16_Q15(QCONST16(.15f, 15), bins[i]);
      for (j=0;j<i;j++)
      {
         if (bins[j] < t1)
            conseq++;
         if (bins[j] < t2)
            conseq++;
         else
            conseq = 0;
      }
      if (conseq>=3)
         is_transient=1;
      conseq = 0;
      for (j=i+1;j<N;j++)
      {
         if (bins[j] < t3)
            conseq++;
         else
            conseq = 0;
      }
      if (conseq>=7)
         is_transient=1;
357
   }
358
   RESTORE_STACK;
359
   return is_transient;
360
361
}

362
363
/** Apply window and compute the MDCT for all sub-frames and 
    all channels in a frame */
364
static void compute_mdcts(const CELTMode *mode, int shortBlocks, celt_sig * restrict in, celt_sig * restrict out, int _C, int LM)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
365
{
366
   const int C = CHANNELS(_C);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
367
   if (C==1 && !shortBlocks)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
368
   {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
369
      const int overlap = OVERLAP(mode);
370
      clt_mdct_forward(&mode->mdct, in, out, mode->window, overlap, mode->maxLM-LM);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
371
   } else {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
372
      const int overlap = OVERLAP(mode);
373
      int N = mode->shortMdctSize<<LM;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
374
375
      int B = 1;
      int b, c;
376
      VARDECL(celt_word32, tmp);
377
      SAVE_STACK;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
378
      if (shortBlocks)
379
      {
380
         /*lookup = &mode->mdct[0];*/
Jean-Marc Valin's avatar
Jean-Marc Valin committed
381
         N = mode->shortMdctSize;
382
         B = shortBlocks;
383
      }
384
      ALLOC(tmp, N, celt_word32);
385
      c=0; do {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
386
387
388
         for (b=0;b<B;b++)
         {
            int j;
389
            clt_mdct_forward(&mode->mdct, in+c*(B*N+overlap)+b*N, tmp, mode->window, overlap, shortBlocks ? mode->maxLM : mode->maxLM-LM);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
390
391
            /* Interleaving the sub-frames */
            for (j=0;j<N;j++)
392
               out[(j*B+b)+c*N*B] = tmp[j];
Jean-Marc Valin's avatar
Jean-Marc Valin committed
393
         }
394
      } while (++c<C);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
395
      RESTORE_STACK;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
396
397
398
   }
}

399
400
/** Compute the IMDCT and apply window for all sub-frames and 
    all channels in a frame */
401
static void compute_inv_mdcts(const CELTMode *mode, int shortBlocks, celt_sig *X,
402
      celt_sig * restrict out_mem[],
403
      celt_sig * restrict overlap_mem[], int _C, int LM)
404
{
405
   int c;
406
   const int C = CHANNELS(_C);
407
   const int N = mode->shortMdctSize<<LM;
408
   const int overlap = OVERLAP(mode);
409
   c=0; do {
410
      int j;
411
412
         VARDECL(celt_word32, x);
         VARDECL(celt_word32, tmp);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
413
414
415
         int b;
         int N2 = N;
         int B = 1;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
416
         SAVE_STACK;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
417
         
Jean-Marc Valin's avatar
Jean-Marc Valin committed
418
         ALLOC(x, N+overlap, celt_word32);
419
         ALLOC(tmp, N, celt_word32);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
420
421
422
423

         if (shortBlocks)
         {
            N2 = mode->shortMdctSize;
424
            B = shortBlocks;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
425
         }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
426
         /* Prevents problems from the imdct doing the overlap-add */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
427
         CELT_MEMSET(x, 0, overlap);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
428

Jean-Marc Valin's avatar
Jean-Marc Valin committed
429
430
431
432
         for (b=0;b<B;b++)
         {
            /* De-interleaving the sub-frames */
            for (j=0;j<N2;j++)
433
               tmp[j] = X[(j*B+b)+c*N2*B];
Jean-Marc Valin's avatar
Jean-Marc Valin committed
434
            clt_mdct_backward(&mode->mdct, tmp, x+N2*b, mode->window, overlap, shortBlocks ? mode->maxLM : mode->maxLM-LM);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
435
         }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
436

437
         for (j=0;j<overlap;j++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
438
            out_mem[c][j] = x[j] + overlap_mem[c][j];
439
         for (;j<N;j++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
440
            out_mem[c][j] = x[j];
441
         for (j=0;j<overlap;j++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
442
            overlap_mem[c][j] = x[N+j];
443
         RESTORE_STACK;
444
   } while (++c<C);
445
446
}

Jean-Marc Valin's avatar
Jean-Marc Valin committed
447
static void deemphasis(celt_sig *in[], celt_word16 *pcm, int N, int _C, int downsample, const celt_word16 *coef, celt_sig *mem)
448
{
449
   const int C = CHANNELS(_C);
450
   int c;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
451
   int count=0;
452
   c=0; do {
453
      int j;
454
455
456
      celt_sig * restrict x;
      celt_word16  * restrict y;
      celt_sig m = mem[c];
Jean-Marc Valin's avatar
Jean-Marc Valin committed
457
      x =in[c];
458
      y = pcm+c;
459
460
      for (j=0;j<N;j++)
      {
461
462
463
         celt_sig tmp = *x + m;
         m = MULT16_32_Q15(coef[0], tmp)
           - MULT16_32_Q15(coef[1], *x);
464
         tmp = SHL32(MULT16_32_Q15(coef[3], tmp), 2);
465
         x++;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
466
467
468
469
470
471
472
473
         /* Technically the store could be moved outside of the if because
            the stores we don't want will just be overwritten */
         if (++count==downsample)
         {
            *y = SCALEOUT(SIG2WORD16(tmp));
            y+=C;
            count=0;
         }
474
      }
475
      mem[c] = m;
476
   } while (++c<C);
477
478
}

479
480
#ifdef ENABLE_POSTFILTER
static void comb_filter(celt_word32 *y, celt_word32 *x, int T0, int T1, int N,
481
      celt_word16 g0, celt_word16 g1, int tapset0, int tapset1,
482
      const celt_word16 *window, int overlap)
483
484
485
486
{
   int i;
   /* printf ("%d %d %f %f\n", T0, T1, g0, g1); */
   celt_word16 g00, g01, g02, g10, g11, g12;
487
   static const celt_word16 gains[3][3] = {
488
489
490
         {QCONST16(0.3066406250f, 15), QCONST16(0.2170410156f, 15), QCONST16(0.1296386719f, 15)},
         {QCONST16(0.4638671875f, 15), QCONST16(0.2680664062f, 15), QCONST16(0.f, 15)},
         {QCONST16(0.7998046875f, 15), QCONST16(0.1000976562f, 15), QCONST16(0.f, 15)}};
491
492
493
494
495
496
   g00 = MULT16_16_Q15(g0, gains[tapset0][0]);
   g01 = MULT16_16_Q15(g0, gains[tapset0][1]);
   g02 = MULT16_16_Q15(g0, gains[tapset0][2]);
   g10 = MULT16_16_Q15(g1, gains[tapset1][0]);
   g11 = MULT16_16_Q15(g1, gains[tapset1][1]);
   g12 = MULT16_16_Q15(g1, gains[tapset1][2]);
497
498
499
500
501
   for (i=0;i<overlap;i++)
   {
      celt_word16 f;
      f = MULT16_16_Q15(window[i],window[i]);
      y[i] = x[i]
502
503
504
505
506
507
508
509
510
511
               + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g00),x[i-T0])
               + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g01),x[i-T0-1])
               + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g01),x[i-T0+1])
               + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g02),x[i-T0-2])
               + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g02),x[i-T0+2])
               + MULT16_32_Q15(MULT16_16_Q15(f,g10),x[i-T1])
               + MULT16_32_Q15(MULT16_16_Q15(f,g11),x[i-T1-1])
               + MULT16_32_Q15(MULT16_16_Q15(f,g11),x[i-T1+1])
               + MULT16_32_Q15(MULT16_16_Q15(f,g12),x[i-T1-2])
               + MULT16_32_Q15(MULT16_16_Q15(f,g12),x[i-T1+2]);
512
513
514
515

   }
   for (i=overlap;i<N;i++)
      y[i] = x[i]
516
517
518
519
520
               + MULT16_32_Q15(g10,x[i-T1])
               + MULT16_32_Q15(g11,x[i-T1-1])
               + MULT16_32_Q15(g11,x[i-T1+1])
               + MULT16_32_Q15(g12,x[i-T1-2])
               + MULT16_32_Q15(g12,x[i-T1+2]);
521
522
523
}
#endif /* ENABLE_POSTFILTER */

524
static const signed char tf_select_table[4][8] = {
525
      {0, -1, 0, -1,    0,-1, 0,-1},
526
527
      {0, -1, 0, -2,    1, 0, 1,-1},
      {0, -2, 0, -3,    2, 0, 1,-1},
528
      {0, -2, 0, -3,    3, 0, 1,-1},
529
530
};

531
static celt_word32 l1_metric(const celt_norm *tmp, int N, int LM, int width)
532
533
{
   int i, j;
534
   static const celt_word16 sqrtM_1[4] = {Q15ONE, QCONST16(.70710678f,15), QCONST16(0.5f,15), QCONST16(0.35355339f,15)};
535
536
537
538
539
540
541
542
543
544
545
   celt_word32 L1;
   celt_word16 bias;
   L1=0;
   for (i=0;i<1<<LM;i++)
   {
      celt_word32 L2 = 0;
      for (j=0;j<N>>LM;j++)
         L2 = MAC16_16(L2, tmp[(j<<LM)+i], tmp[(j<<LM)+i]);
      L1 += celt_sqrt(L2);
   }
   L1 = MULT16_32_Q15(sqrtM_1[LM], L1);
546
547
548
549
550
551
   if (width==1)
      bias = QCONST16(.12f,15)*LM;
   else if (width==2)
      bias = QCONST16(.05f,15)*LM;
   else
      bias = QCONST16(.02f,15)*LM;
552
553
554
555
   L1 = MAC16_32_Q15(L1, bias, L1);
   return L1;
}

556
static int tf_analysis(const CELTMode *m, celt_word16 *bandLogE, celt_word16 *oldBandE,
557
558
      int len, int C, int isTransient, int *tf_res, int nbCompressedBytes, celt_norm *X,
      int N0, int LM, int *tf_sum)
559
{
560
   int i;
561
562
563
   VARDECL(int, metric);
   int cost0;
   int cost1;
564
565
   VARDECL(int, path0);
   VARDECL(int, path1);
566
567
   VARDECL(celt_norm, tmp);
   int lambda;
568
   int tf_select=0;
569
570
   SAVE_STACK;

571
572
   if (nbCompressedBytes<15*C)
   {
573
      *tf_sum = 0;
574
      for (i=0;i<len;i++)
575
         tf_res[i] = isTransient;
576
577
      return 0;
   }
578
   if (nbCompressedBytes<40)
579
      lambda = 12;
580
   else if (nbCompressedBytes<60)
581
      lambda = 6;
582
   else if (nbCompressedBytes<100)
583
      lambda = 4;
584
   else
585
      lambda = 3;
586

587
588
   ALLOC(metric, len, int);
   ALLOC(tmp, (m->eBands[len]-m->eBands[len-1])<<LM, celt_norm);
589
590
   ALLOC(path0, len, int);
   ALLOC(path1, len, int);
591

592
   *tf_sum = 0;
593
   for (i=0;i<len;i++)
594
   {
595
596
597
598
599
600
      int j, k, N;
      celt_word32 L1, best_L1;
      int best_level=0;
      N = (m->eBands[i+1]-m->eBands[i])<<LM;
      for (j=0;j<N;j++)
         tmp[j] = X[j+(m->eBands[i]<<LM)];
601
      /* Just add the right channel if we're in stereo */
602
      if (C==2)
603
         for (j=0;j<N;j++)
604
            tmp[j] = ADD16(tmp[j],X[N0+j+(m->eBands[i]<<LM)]);
605
      L1 = l1_metric(tmp, N, isTransient ? LM : 0, N>>LM);
606
607
608
      best_L1 = L1;
      /*printf ("%f ", L1);*/
      for (k=0;k<LM;k++)
609
      {
610
611
612
613
614
615
616
         int B;

         if (isTransient)
            B = (LM-k-1);
         else
            B = k+1;

617
618
619
620
621
         if (isTransient)
            haar1(tmp, N>>(LM-k), 1<<(LM-k));
         else
            haar1(tmp, N>>k, 1<<k);

622
         L1 = l1_metric(tmp, N, B, N>>LM);
623
624
625
626
627
628

         if (L1 < best_L1)
         {
            best_L1 = L1;
            best_level = k+1;
         }
629
      }
630
631
632
633
634
      /*printf ("%d ", isTransient ? LM-best_level : best_level);*/
      if (isTransient)
         metric[i] = best_level;
      else
         metric[i] = -best_level;
635
      *tf_sum += metric[i];
636
   }
637
638
   /*printf("\n");*/
   /* FIXME: Figure out how to set this */
639
   tf_select = 0;
640

641
   cost0 = 0;
642
   cost1 = isTransient ? 0 : lambda;
643
644
   /* Viterbi forward pass */
   for (i=1;i<len;i++)
645
   {
646
647
      int curr0, curr1;
      int from0, from1;
648

649
650
      from0 = cost0;
      from1 = cost1 + lambda;
651
652
      if (from0 < from1)
      {
653
         curr0 = from0;
654
655
         path0[i]= 0;
      } else {
656
         curr0 = from1;
657
658
659
         path0[i]= 1;
      }

660
661
      from0 = cost0 + lambda;
      from1 = cost1;
662
663
      if (from0 < from1)
      {
664
         curr1 = from0;
665
666
         path1[i]= 0;
      } else {
667
         curr1 = from1;
668
669
         path1[i]= 1;
      }
670
671
      cost0 = curr0 + abs(metric[i]-tf_select_table[LM][4*isTransient+2*tf_select+0]);
      cost1 = curr1 + abs(metric[i]-tf_select_table[LM][4*isTransient+2*tf_select+1]);
672
   }
673
   tf_res[len-1] = cost0 < cost1 ? 0 : 1;
674
675
676
677
678
   /* Viterbi backward pass to check the decisions */
   for (i=len-2;i>=0;i--)
   {
      if (tf_res[i+1] == 1)
         tf_res[i] = path1[i+1];
679
      else
680
         tf_res[i] = path0[i+1];
681
   }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
682
   RESTORE_STACK;
683
   return tf_select;
684
685
}

686
static void tf_encode(int start, int end, int isTransient, int *tf_res, int LM, int tf_select, ec_enc *enc)
687
688
{
   int curr, i;
689
690
691
692
693
   int tf_select_rsv;
   int tf_changed;
   int logp;
   ec_uint32 budget;
   ec_uint32 tell;
694
695
   budget = enc->storage*8;
   tell = ec_tell(enc);
696
697
698
699
700
701
   logp = isTransient ? 2 : 4;
   /* Reserve space to code the tf_select decision. */
   tf_select_rsv = LM>0 && tell+logp+1 <= budget;
   budget -= tf_select_rsv;
   curr = tf_changed = 0;
   for (i=start;i<end;i++)
702
   {
703
704
705
      if (tell+logp<=budget)
      {
         ec_enc_bit_logp(enc, tf_res[i] ^ curr, logp);
706
         tell = ec_tell(enc);
707
708
709
710
711
712
         curr = tf_res[i];
         tf_changed |= curr;
      }
      else
         tf_res[i] = curr;
      logp = isTransient ? 4 : 5;
713
   }
714
715
716
717
718
719
720
   /* Only code tf_select if it would actually make a difference. */
   if (tf_select_rsv &&
         tf_select_table[LM][4*isTransient+0+tf_changed]!=
         tf_select_table[LM][4*isTransient+2+tf_changed])
      ec_enc_bit_logp(enc, tf_select, 1);
   else
      tf_select = 0;
721
   for (i=start;i<end;i++)
722
      tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
723
   /*printf("%d %d ", isTransient, tf_select); for(i=0;i<end;i++)printf("%d ", tf_res[i]);printf("\n");*/
724
725
}

726
static void tf_decode(int start, int end, int isTransient, int *tf_res, int LM, ec_dec *dec)
727
{
728
   int i, curr, tf_select;
729
730
731
732
733
734
   int tf_select_rsv;
   int tf_changed;
   int logp;
   ec_uint32 budget;
   ec_uint32 tell;

735
736
   budget = dec->storage*8;
   tell = ec_tell(dec);
737
738
739
740
741
742
743
744
745
   logp = isTransient ? 2 : 4;
   tf_select_rsv = LM>0 && tell+logp+1<=budget;
   budget -= tf_select_rsv;
   tf_changed = curr = 0;
   for (i=start;i<end;i++)
   {
      if (tell+logp<=budget)
      {
         curr ^= ec_dec_bit_logp(dec, logp);
746
         tell = ec_tell(dec);
747
748
749
750
751
752
753
754
755
756
         tf_changed |= curr;
      }
      tf_res[i] = curr;
      logp = isTransient ? 4 : 5;
   }
   tf_select = 0;
   if (tf_select_rsv &&
     tf_select_table[LM][4*isTransient+0+tf_changed] !=
     tf_select_table[LM][4*isTransient+2+tf_changed])
   {
757
      tf_select = ec_dec_bit_logp(dec, 1);
758
759
   }
   for (i=start;i<end;i++)
760
   {
761
      tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
762
   }
763
764
}

765
766
767
768
769
770
771
772
773
774
775
static void init_caps(const CELTMode *m,int *cap,int LM,int C)
{
   int i;
   for (i=0;i<m->nbEBands;i++)
   {
      int N;
      N=(m->eBands[i+1]-m->eBands[i])<<LM;
      cap[i] = (m->cache.caps[m->nbEBands*(2*LM+C-1)+i]+64)*C*N>>2;
   }
}

776
777
778
779
static int alloc_trim_analysis(const CELTMode *m, const celt_norm *X,
      const celt_word16 *bandLogE, int nbEBands, int LM, int C, int N0)
{
   int i;
780
781
   celt_word32 diff=0;
   int c;
782
   int trim_index = 5;
783
784
785
786
787
788
789
790
791
792
793
794
795
796
   if (C==2)
   {
      celt_word16 sum = 0; /* Q10 */
      /* Compute inter-channel correlation for low frequencies */
      for (i=0;i<8;i++)
      {
         int j;
         celt_word32 partial = 0;
         for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
            partial = MAC16_16(partial, X[j], X[N0+j]);
         sum = ADD16(sum, EXTRACT16(SHR32(partial, 18)));
      }
      sum = MULT16_16_Q15(QCONST16(1.f/8, 15), sum);
      /*printf ("%f\n", sum);*/
797
      if (sum > QCONST16(.995f,10))
798
         trim_index-=4;
799
      else if (sum > QCONST16(.92f,10))
800
801
         trim_index-=3;
      else if (sum > QCONST16(.85f,10))
802
         trim_index-=2;
803
      else if (sum > QCONST16(.8f,10))
804
         trim_index-=1;
805
   }
806
807

   /* Estimate spectral tilt */
808
   c=0; do {
809
810
      for (i=0;i<nbEBands-1;i++)
      {
811
         diff += bandLogE[i+c*nbEBands]*(celt_int32)(2+2*i-nbEBands);
812
      }
813
   } while (++c<0);
814
815
   diff /= C*(nbEBands-1);
   /*printf("%f\n", diff);*/
816
   if (diff > QCONST16(2.f, DB_SHIFT))
817
      trim_index--;
818
   if (diff > QCONST16(8.f, DB_SHIFT))
819
      trim_index--;
820
   if (diff < -QCONST16(4.f, DB_SHIFT))
Jean-Marc Valin's avatar
Jean-Marc Valin committed
821
      trim_index++;
822
823
824
   if (diff < -QCONST16(10.f, DB_SHIFT))
      trim_index++;

825
826
   if (trim_index<0)
      trim_index = 0;
827
828
   if (trim_index>10)
      trim_index = 10;
829
830
831
   return trim_index;
}

Jean-Marc Valin's avatar
Jean-Marc Valin committed
832
static int stereo_analysis(const CELTMode *m, const celt_norm *X,
833
      int LM, int N0)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
{
   int i;
   int thetas;
   celt_word32 sumLR = EPSILON, sumMS = EPSILON;

   /* Use the L1 norm to model the entropy of the L/R signal vs the M/S signal */
   for (i=0;i<13;i++)
   {
      int j;
      for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
      {
         celt_word16 L, R, M, S;
         L = X[j];
         R = X[N0+j];
         M = L+R;
         S = L-R;
         sumLR += EXTEND32(ABS16(L)) + EXTEND32(ABS16(R));
         sumMS += EXTEND32(ABS16(M)) + EXTEND32(ABS16(S));
      }
   }
   sumMS = MULT16_32_Q15(QCONST16(0.707107f, 15), sumMS);
   thetas = 13;
   /* We don't need thetas for lower bands with LM<=1 */
   if (LM<=1)
      thetas -= 8;
   return MULT16_32_Q15((m->eBands[13]<<(LM+1))+thetas, sumMS)
         > MULT16_32_Q15(m->eBands[13]<<(LM+1), sumLR);
}

863
#ifdef FIXED_POINT
864
CELT_STATIC
865
int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
866
{
867
#else
868
CELT_STATIC
869
int celt_encode_with_ec_float(CELTEncoder * restrict st, const celt_sig * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
870
871
{
#endif
872
   int i, c, N;
873
   int bits;
874
   ec_enc _enc;
875
876
877
878
879
   VARDECL(celt_sig, in);
   VARDECL(celt_sig, freq);
   VARDECL(celt_norm, X);
   VARDECL(celt_ener, bandE);
   VARDECL(celt_word16, bandLogE);
880
   VARDECL(int, fine_quant);
881
   VARDECL(celt_word16, error);
882
   VARDECL(int, pulses);
883
   VARDECL(int, cap);
884
   VARDECL(int, offsets);
885
   VARDECL(int, fine_priority);
886
   VARDECL(int, tf_res);
Timothy B. Terriberry's avatar