celt.c 76.6 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 clip;
108
   int disable_pf;
109
   int complexity;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
110
   int upsample;
111
112
   int start, end;

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

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

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

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

138
139
140
141
142
143
   /* 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
144
145
   celt_word32 preemph_memE[2];
   celt_word32 preemph_memD[2];
146

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

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

157
158
159
160
161
162
163
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)
164
165
166
{
   int size = sizeof(struct CELTEncoder)
         + (2*channels*mode->overlap-1)*sizeof(celt_sig)
167
         + channels*COMBFILTER_MAXPERIOD*sizeof(celt_sig)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
168
         + 3*channels*mode->nbEBands*sizeof(celt_word16);
169
170
171
   return size;
}

Jean-Marc Valin's avatar
Jean-Marc Valin committed
172
CELTEncoder *celt_encoder_create(int sampling_rate, int channels, int *error)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
173
{
174
   CELTEncoder *st;
175
   st = (CELTEncoder *)celt_alloc(celt_encoder_get_size(channels));
176
177
178
179
180
   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
181
   return st;
182
183
184
}

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

Jean-Marc Valin's avatar
Jean-Marc Valin committed
195
CELTEncoder *celt_encoder_init(CELTEncoder *st, int sampling_rate, int channels, int *error)
196
{
Jean-Marc Valin's avatar
Jean-Marc Valin committed
197
198
199
200
201
202
203
204
205
   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;
206
207
208
}

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

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

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

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

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

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

252
void celt_encoder_destroy(CELTEncoder *st)
253
254
255
256
{
   celt_free(st);
}

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

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

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

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

   /* 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;
313
      mem1 = x - .5f*y;
314
315
#endif
      tmp[i] = EXTRACT16(SHR(y,2));
316
   }
317
   /* First few samples are bad because we don't propagate the memory */
318
   for (i=0;i<12;i++)
319
320
      tmp[i] = 0;

321
   for (i=0;i<N;i++)
322
   {
323
324
325
326
327
      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;
328
   }
329
   for (i=0;i<N;i++)
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
357
358
      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;
359
   }
360
   RESTORE_STACK;
361
   return is_transient;
362
363
}

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

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

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

Jean-Marc Valin's avatar
Jean-Marc Valin committed
431
432
433
434
         for (b=0;b<B;b++)
         {
            /* De-interleaving the sub-frames */
            for (j=0;j<N2;j++)
435
               tmp[j] = X[(j*B+b)+c*N2*B];
Jean-Marc Valin's avatar
Jean-Marc Valin committed
436
            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
437
         }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
438

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

Jean-Marc Valin's avatar
Jean-Marc Valin committed
449
static void deemphasis(celt_sig *in[], celt_word16 *pcm, int N, int _C, int downsample, const celt_word16 *coef, celt_sig *mem)
450
{
451
   const int C = CHANNELS(_C);
452
   int c;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
453
   int count=0;
454
   c=0; do {
455
      int j;
456
457
458
      celt_sig * restrict x;
      celt_word16  * restrict y;
      celt_sig m = mem[c];
Jean-Marc Valin's avatar
Jean-Marc Valin committed
459
      x =in[c];
460
      y = pcm+c;
461
462
      for (j=0;j<N;j++)
      {
463
464
465
         celt_sig tmp = *x + m;
         m = MULT16_32_Q15(coef[0], tmp)
           - MULT16_32_Q15(coef[1], *x);
466
         tmp = SHL32(MULT16_32_Q15(coef[3], tmp), 2);
467
         x++;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
468
469
470
471
472
473
474
475
         /* 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;
         }
476
      }
477
      mem[c] = m;
478
   } while (++c<C);
479
480
}

481
482
#ifdef ENABLE_POSTFILTER
static void comb_filter(celt_word32 *y, celt_word32 *x, int T0, int T1, int N,
483
      celt_word16 g0, celt_word16 g1, int tapset0, int tapset1,
484
      const celt_word16 *window, int overlap)
485
486
487
488
{
   int i;
   /* printf ("%d %d %f %f\n", T0, T1, g0, g1); */
   celt_word16 g00, g01, g02, g10, g11, g12;
489
   static const celt_word16 gains[3][3] = {
490
491
492
         {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)}};
493
494
495
496
497
498
   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]);
499
500
501
502
503
   for (i=0;i<overlap;i++)
   {
      celt_word16 f;
      f = MULT16_16_Q15(window[i],window[i]);
      y[i] = x[i]
504
505
506
507
508
509
510
511
512
513
               + 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]);
514
515
516
517

   }
   for (i=overlap;i<N;i++)
      y[i] = x[i]
518
519
520
521
522
               + 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]);
523
524
525
}
#endif /* ENABLE_POSTFILTER */

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

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

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

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

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

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

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

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

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

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

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

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

662
663
      from0 = cost0 + lambda;
      from1 = cost1;
664
665
      if (from0 < from1)
      {
666
         curr1 = from0;
667
668
         path1[i]= 0;
      } else {
669
         curr1 = from1;
670
671
         path1[i]= 1;
      }
672
673
      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]);
674
   }
675
   tf_res[len-1] = cost0 < cost1 ? 0 : 1;
676
677
678
679
680
   /* 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];
681
      else
682
         tf_res[i] = path0[i+1];
683
   }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
684
   RESTORE_STACK;
685
   return tf_select;
686
687
}

688
static void tf_encode(int start, int end, int isTransient, int *tf_res, int LM, int tf_select, ec_enc *enc)
689
690
{
   int curr, i;
691
692
693
694
695
   int tf_select_rsv;
   int tf_changed;
   int logp;
   ec_uint32 budget;
   ec_uint32 tell;
696
697
   budget = enc->storage*8;
   tell = ec_tell(enc);
698
699
700
701
702
703
   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++)
704
   {
705
706
707
      if (tell+logp<=budget)
      {
         ec_enc_bit_logp(enc, tf_res[i] ^ curr, logp);
708
         tell = ec_tell(enc);
709
710
711
712
713
714
         curr = tf_res[i];
         tf_changed |= curr;
      }
      else
         tf_res[i] = curr;
      logp = isTransient ? 4 : 5;
715
   }
716
717
718
719
720
721
722
   /* 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;
723
   for (i=start;i<end;i++)
724
      tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
725
   /*printf("%d %d ", isTransient, tf_select); for(i=0;i<end;i++)printf("%d ", tf_res[i]);printf("\n");*/
726
727
}

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

737
738
   budget = dec->storage*8;
   tell = ec_tell(dec);
739
740
741
742
743
744
745
746
747
   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);
748
         tell = ec_tell(dec);
749
750
751
752
753
754
755
756
757
758
         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])
   {
759
      tf_select = ec_dec_bit_logp(dec, 1);
760
761
   }
   for (i=start;i<end;i++)
762
   {
763
      tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
764
   }
765
766
}

767
768
769
770
771
772
773
774
775
776
777
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;
   }
}

778
static int alloc_trim_analysis(const CELTMode *m, const celt_norm *X,
779
      const celt_word16 *bandLogE, int end, int LM, int C, int N0)
780
781
{
   int i;
782
783
   celt_word32 diff=0;
   int c;
784
   int trim_index = 5;
785
786
787
788
789
790
791
792
793
794
795
796
797
798
   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);*/
799
      if (sum > QCONST16(.995f,10))
800
         trim_index-=4;
801
      else if (sum > QCONST16(.92f,10))
802
803
         trim_index-=3;
      else if (sum > QCONST16(.85f,10))
804
         trim_index-=2;
805
      else if (sum > QCONST16(.8f,10))
806
         trim_index-=1;
807
   }
808
809

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

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

Jean-Marc Valin's avatar
Jean-Marc Valin committed
834
static int stereo_analysis(const CELTMode *m, const celt_norm *X,
835
      int LM, int N0)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
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
863
864
{
   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);
}

865
#ifdef FIXED_POINT
866
CELT_STATIC
867
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
868
{
869
#else
870
CELT_STATIC
871
int celt_encode_with_ec_float(CELTEncoder * restrict st, const celt_sig * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
872
873
{
#endif
874
   int i, c, N;
Timothy B. Terriberry's avatar
Timothy B. Terriberry committed
875
   celt_int32 bits;
876
   ec_enc _enc;
877
878
879
880
881
   VARDECL(celt_sig, in);
   VARDECL(celt_sig, freq);
   VARDECL(celt_norm, X);
   VARDECL(celt_ener, bandE);
   VARDECL(celt_word16, bandLogE);
882
   VARDECL(int, fine_quant);
883
   VARDECL(celt_word16, error);