git.sesse.net Git - ffmpeg/blob - libavcodec/opus_celt.c

   1 /*
   2  * Copyright (c) 2012 Andrew D'Addesio
   3  * Copyright (c) 2013-2014 Mozilla Corporation
   4  *
   5  * This file is part of FFmpeg.
   6  *
   7  * FFmpeg is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Lesser General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2.1 of the License, or (at your option) any later version.
  11  *
  12  * FFmpeg is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public
  18  * License along with FFmpeg; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20  */
  21
  22 /**
  23  * @file
  24  * Opus CELT decoder
  25  */
  26
  27 #include <stdint.h>
  28
  29 #include "libavutil/float_dsp.h"
  30 #include "libavutil/libm.h"
  31
  32 #include "imdct15.h"
  33 #include "opus.h"
  34
  35 enum CeltSpread {
  36     CELT_SPREAD_NONE,
  37     CELT_SPREAD_LIGHT,
  38     CELT_SPREAD_NORMAL,
  39     CELT_SPREAD_AGGRESSIVE
  40 };
  41
  42 typedef struct CeltFrame {
  43     float energy[CELT_MAX_BANDS];
  44     float prev_energy[2][CELT_MAX_BANDS];
  45
  46     uint8_t collapse_masks[CELT_MAX_BANDS];
  47
  48     /* buffer for mdct output + postfilter */
  49     DECLARE_ALIGNED(32, float, buf)[2048];
  50
  51     /* postfilter parameters */
  52     int pf_period_new;
  53     float pf_gains_new[3];
  54     int pf_period;
  55     float pf_gains[3];
  56     int pf_period_old;
  57     float pf_gains_old[3];
  58
  59     float deemph_coeff;
  60 } CeltFrame;
  61
  62 struct CeltContext {
  63     // constant values that do not change during context lifetime
  64     AVCodecContext    *avctx;
  65     IMDCT15Context    *imdct[4];
  66     AVFloatDSPContext  *dsp;
  67     int output_channels;
  68
  69     // values that have inter-frame effect and must be reset on flush
  70     CeltFrame frame[2];
  71     uint32_t seed;
  72     int flushed;
  73
  74     // values that only affect a single frame
  75     int coded_channels;
  76     int framebits;
  77     int duration;
  78
  79     /* number of iMDCT blocks in the frame */
  80     int blocks;
  81     /* size of each block */
  82     int blocksize;
  83
  84     int startband;
  85     int endband;
  86     int codedbands;
  87
  88     int anticollapse_bit;
  89
  90     int intensitystereo;
  91     int dualstereo;
  92     enum CeltSpread spread;
  93
  94     int remaining;
  95     int remaining2;
  96     int fine_bits    [CELT_MAX_BANDS];
  97     int fine_priority[CELT_MAX_BANDS];
  98     int pulses       [CELT_MAX_BANDS];
  99     int tf_change    [CELT_MAX_BANDS];
 100
 101     DECLARE_ALIGNED(32, float, coeffs)[2][CELT_MAX_FRAME_SIZE];
 102     DECLARE_ALIGNED(32, float, scratch)[22 * 8]; // MAX(celt_freq_range) * 1<<CELT_MAX_LOG_BLOCKS
 103 };
 104
 105 static const uint16_t celt_model_tapset[] = { 4, 2, 3, 4 };
 106
 107 static const uint16_t celt_model_spread[] = { 32, 7, 9, 30, 32 };
 108
 109 static const uint16_t celt_model_alloc_trim[] = {
 110     128,   2,   4,   9,  19,  41,  87, 109, 119, 124, 126, 128
 111 };
 112
 113 static const uint16_t celt_model_energy_small[] = { 4, 2, 3, 4 };
 114
 115 static const uint8_t celt_freq_bands[] = { /* in steps of 200Hz */
 116     0,  1,  2,  3,  4,  5,  6,  7,  8, 10, 12, 14, 16, 20, 24, 28, 34, 40, 48, 60, 78, 100
 117 };
 118
 119 static const uint8_t celt_freq_range[] = {
 120     1,  1,  1,  1,  1,  1,  1,  1,  2,  2,  2,  2,  4,  4,  4,  6,  6,  8, 12, 18, 22
 121 };
 122
 123 static const uint8_t celt_log_freq_range[] = {
 124     0,  0,  0,  0,  0,  0,  0,  0,  8,  8,  8,  8, 16, 16, 16, 21, 21, 24, 29, 34, 36
 125 };
 126
 127 static const int8_t celt_tf_select[4][2][2][2] = {
 128     { { { 0, -1 }, { 0, -1 } }, { { 0, -1 }, { 0, -1 } } },
 129     { { { 0, -1 }, { 0, -2 } }, { { 1,  0 }, { 1, -1 } } },
 130     { { { 0, -2 }, { 0, -3 } }, { { 2,  0 }, { 1, -1 } } },
 131     { { { 0, -2 }, { 0, -3 } }, { { 3,  0 }, { 1, -1 } } }
 132 };
 133
 134 static const float celt_mean_energy[] = {
 135     6.437500f, 6.250000f, 5.750000f, 5.312500f, 5.062500f,
 136     4.812500f, 4.500000f, 4.375000f, 4.875000f, 4.687500f,
 137     4.562500f, 4.437500f, 4.875000f, 4.625000f, 4.312500f,
 138     4.500000f, 4.375000f, 4.625000f, 4.750000f, 4.437500f,
 139     3.750000f, 3.750000f, 3.750000f, 3.750000f, 3.750000f
 140 };
 141
 142 static const float celt_alpha_coef[] = {
 143     29440.0f/32768.0f,    26112.0f/32768.0f,    21248.0f/32768.0f,    16384.0f/32768.0f
 144 };
 145
 146 static const float celt_beta_coef[] = { /* TODO: precompute 1 minus this if the code ends up neater */
 147     30147.0f/32768.0f,    22282.0f/32768.0f,    12124.0f/32768.0f,     6554.0f/32768.0f
 148 };
 149
 150 static const uint8_t celt_coarse_energy_dist[4][2][42] = {
 151     {
 152         {       // 120-sample inter
 153              72, 127,  65, 129,  66, 128,  65, 128,  64, 128,  62, 128,  64, 128,
 154              64, 128,  92,  78,  92,  79,  92,  78,  90,  79, 116,  41, 115,  40,
 155             114,  40, 132,  26, 132,  26, 145,  17, 161,  12, 176,  10, 177,  11
 156         }, {    // 120-sample intra
 157              24, 179,  48, 138,  54, 135,  54, 132,  53, 134,  56, 133,  55, 132,
 158              55, 132,  61, 114,  70,  96,  74,  88,  75,  88,  87,  74,  89,  66,
 159              91,  67, 100,  59, 108,  50, 120,  40, 122,  37,  97,  43,  78,  50
 160         }
 161     }, {
 162         {       // 240-sample inter
 163              83,  78,  84,  81,  88,  75,  86,  74,  87,  71,  90,  73,  93,  74,
 164              93,  74, 109,  40, 114,  36, 117,  34, 117,  34, 143,  17, 145,  18,
 165             146,  19, 162,  12, 165,  10, 178,   7, 189,   6, 190,   8, 177,   9
 166         }, {    // 240-sample intra
 167              23, 178,  54, 115,  63, 102,  66,  98,  69,  99,  74,  89,  71,  91,
 168              73,  91,  78,  89,  86,  80,  92,  66,  93,  64, 102,  59, 103,  60,
 169             104,  60, 117,  52, 123,  44, 138,  35, 133,  31,  97,  38,  77,  45
 170         }
 171     }, {
 172         {       // 480-sample inter
 173              61,  90,  93,  60, 105,  42, 107,  41, 110,  45, 116,  38, 113,  38,
 174             112,  38, 124,  26, 132,  27, 136,  19, 140,  20, 155,  14, 159,  16,
 175             158,  18, 170,  13, 177,  10, 187,   8, 192,   6, 175,   9, 159,  10
 176         }, {    // 480-sample intra
 177              21, 178,  59, 110,  71,  86,  75,  85,  84,  83,  91,  66,  88,  73,
 178              87,  72,  92,  75,  98,  72, 105,  58, 107,  54, 115,  52, 114,  55,
 179             112,  56, 129,  51, 132,  40, 150,  33, 140,  29,  98,  35,  77,  42
 180         }
 181     }, {
 182         {       // 960-sample inter
 183              42, 121,  96,  66, 108,  43, 111,  40, 117,  44, 123,  32, 120,  36,
 184             119,  33, 127,  33, 134,  34, 139,  21, 147,  23, 152,  20, 158,  25,
 185             154,  26, 166,  21, 173,  16, 184,  13, 184,  10, 150,  13, 139,  15
 186         }, {    // 960-sample intra
 187              22, 178,  63, 114,  74,  82,  84,  83,  92,  82, 103,  62,  96,  72,
 188              96,  67, 101,  73, 107,  72, 113,  55, 118,  52, 125,  52, 118,  52,
 189             117,  55, 135,  49, 137,  39, 157,  32, 145,  29,  97,  33,  77,  40
 190         }
 191     }
 192 };
 193
 194 static const uint8_t celt_static_alloc[11][21] = {  /* 1/32 bit/sample */
 195     {   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0 },
 196     {  90,  80,  75,  69,  63,  56,  49,  40,  34,  29,  20,  18,  10,   0,   0,   0,   0,   0,   0,   0,   0 },
 197     { 110, 100,  90,  84,  78,  71,  65,  58,  51,  45,  39,  32,  26,  20,  12,   0,   0,   0,   0,   0,   0 },
 198     { 118, 110, 103,  93,  86,  80,  75,  70,  65,  59,  53,  47,  40,  31,  23,  15,   4,   0,   0,   0,   0 },
 199     { 126, 119, 112, 104,  95,  89,  83,  78,  72,  66,  60,  54,  47,  39,  32,  25,  17,  12,   1,   0,   0 },
 200     { 134, 127, 120, 114, 103,  97,  91,  85,  78,  72,  66,  60,  54,  47,  41,  35,  29,  23,  16,  10,   1 },
 201     { 144, 137, 130, 124, 113, 107, 101,  95,  88,  82,  76,  70,  64,  57,  51,  45,  39,  33,  26,  15,   1 },
 202     { 152, 145, 138, 132, 123, 117, 111, 105,  98,  92,  86,  80,  74,  67,  61,  55,  49,  43,  36,  20,   1 },
 203     { 162, 155, 148, 142, 133, 127, 121, 115, 108, 102,  96,  90,  84,  77,  71,  65,  59,  53,  46,  30,   1 },
 204     { 172, 165, 158, 152, 143, 137, 131, 125, 118, 112, 106, 100,  94,  87,  81,  75,  69,  63,  56,  45,  20 },
 205     { 200, 200, 200, 200, 200, 200, 200, 200, 198, 193, 188, 183, 178, 173, 168, 163, 158, 153, 148, 129, 104 }
 206 };
 207
 208 static const uint8_t celt_static_caps[4][2][21] = {
 209     {       // 120-sample
 210         {224, 224, 224, 224, 224, 224, 224, 224, 160, 160,
 211          160, 160, 185, 185, 185, 178, 178, 168, 134,  61,  37},
 212         {224, 224, 224, 224, 224, 224, 224, 224, 240, 240,
 213          240, 240, 207, 207, 207, 198, 198, 183, 144,  66,  40},
 214     }, {    // 240-sample
 215         {160, 160, 160, 160, 160, 160, 160, 160, 185, 185,
 216          185, 185, 193, 193, 193, 183, 183, 172, 138,  64,  38},
 217         {240, 240, 240, 240, 240, 240, 240, 240, 207, 207,
 218          207, 207, 204, 204, 204, 193, 193, 180, 143,  66,  40},
 219     }, {    // 480-sample
 220         {185, 185, 185, 185, 185, 185, 185, 185, 193, 193,
 221          193, 193, 193, 193, 193, 183, 183, 172, 138,  65,  39},
 222         {207, 207, 207, 207, 207, 207, 207, 207, 204, 204,
 223          204, 204, 201, 201, 201, 188, 188, 176, 141,  66,  40},
 224     }, {    // 960-sample
 225         {193, 193, 193, 193, 193, 193, 193, 193, 193, 193,
 226          193, 193, 194, 194, 194, 184, 184, 173, 139,  65,  39},
 227         {204, 204, 204, 204, 204, 204, 204, 204, 201, 201,
 228          201, 201, 198, 198, 198, 187, 187, 175, 140,  66,  40}
 229     }
 230 };
 231
 232 static const uint8_t celt_cache_bits[392] = {
 233     40, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 234     7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 235     7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 40, 15, 23, 28,
 236     31, 34, 36, 38, 39, 41, 42, 43, 44, 45, 46, 47, 47, 49, 50,
 237     51, 52, 53, 54, 55, 55, 57, 58, 59, 60, 61, 62, 63, 63, 65,
 238     66, 67, 68, 69, 70, 71, 71, 40, 20, 33, 41, 48, 53, 57, 61,
 239     64, 66, 69, 71, 73, 75, 76, 78, 80, 82, 85, 87, 89, 91, 92,
 240     94, 96, 98, 101, 103, 105, 107, 108, 110, 112, 114, 117, 119, 121, 123,
 241     124, 126, 128, 40, 23, 39, 51, 60, 67, 73, 79, 83, 87, 91, 94,
 242     97, 100, 102, 105, 107, 111, 115, 118, 121, 124, 126, 129, 131, 135, 139,
 243     142, 145, 148, 150, 153, 155, 159, 163, 166, 169, 172, 174, 177, 179, 35,
 244     28, 49, 65, 78, 89, 99, 107, 114, 120, 126, 132, 136, 141, 145, 149,
 245     153, 159, 165, 171, 176, 180, 185, 189, 192, 199, 205, 211, 216, 220, 225,
 246     229, 232, 239, 245, 251, 21, 33, 58, 79, 97, 112, 125, 137, 148, 157,
 247     166, 174, 182, 189, 195, 201, 207, 217, 227, 235, 243, 251, 17, 35, 63,
 248     86, 106, 123, 139, 152, 165, 177, 187, 197, 206, 214, 222, 230, 237, 250,
 249     25, 31, 55, 75, 91, 105, 117, 128, 138, 146, 154, 161, 168, 174, 180,
 250     185, 190, 200, 208, 215, 222, 229, 235, 240, 245, 255, 16, 36, 65, 89,
 251     110, 128, 144, 159, 173, 185, 196, 207, 217, 226, 234, 242, 250, 11, 41,
 252     74, 103, 128, 151, 172, 191, 209, 225, 241, 255, 9, 43, 79, 110, 138,
 253     163, 186, 207, 227, 246, 12, 39, 71, 99, 123, 144, 164, 182, 198, 214,
 254     228, 241, 253, 9, 44, 81, 113, 142, 168, 192, 214, 235, 255, 7, 49,
 255     90, 127, 160, 191, 220, 247, 6, 51, 95, 134, 170, 203, 234, 7, 47,
 256     87, 123, 155, 184, 212, 237, 6, 52, 97, 137, 174, 208, 240, 5, 57,
 257     106, 151, 192, 231, 5, 59, 111, 158, 202, 243, 5, 55, 103, 147, 187,
 258     224, 5, 60, 113, 161, 206, 248, 4, 65, 122, 175, 224, 4, 67, 127,
 259     182, 234
 260 };
 261
 262 static const int16_t celt_cache_index[105] = {
 263     -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 41, 41, 41,
 264     82, 82, 123, 164, 200, 222, 0, 0, 0, 0, 0, 0, 0, 0, 41,
 265     41, 41, 41, 123, 123, 123, 164, 164, 240, 266, 283, 295, 41, 41, 41,
 266     41, 41, 41, 41, 41, 123, 123, 123, 123, 240, 240, 240, 266, 266, 305,
 267     318, 328, 336, 123, 123, 123, 123, 123, 123, 123, 123, 240, 240, 240, 240,
 268     305, 305, 305, 318, 318, 343, 351, 358, 364, 240, 240, 240, 240, 240, 240,
 269     240, 240, 305, 305, 305, 305, 343, 343, 343, 351, 351, 370, 376, 382, 387,
 270 };
 271
 272 static const uint8_t celt_log2_frac[] = {
 273     0, 8, 13, 16, 19, 21, 23, 24, 26, 27, 28, 29, 30, 31, 32, 32, 33, 34, 34, 35, 36, 36, 37, 37
 274 };
 275
 276 static const uint8_t celt_bit_interleave[] = {
 277     0, 1, 1, 1, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3
 278 };
 279
 280 static const uint8_t celt_bit_deinterleave[] = {
 281     0x00, 0x03, 0x0C, 0x0F, 0x30, 0x33, 0x3C, 0x3F,
 282     0xC0, 0xC3, 0xCC, 0xCF, 0xF0, 0xF3, 0xFC, 0xFF
 283 };
 284
 285 static const uint8_t celt_hadamard_ordery[] = {
 286     1,   0,
 287     3,   0,  2,  1,
 288     7,   0,  4,  3,  6,  1,  5,  2,
 289     15,  0,  8,  7, 12,  3, 11,  4, 14,  1,  9,  6, 13,  2, 10,  5
 290 };
 291
 292 static const uint16_t celt_qn_exp2[] = {
 293     16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048
 294 };
 295
 296 static const uint32_t celt_pvq_u[1272] = {
 297     /* N = 0, K = 0...176 */
 298     1, 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,
 299     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, 0,
 300     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, 0,
 301     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, 0,
 302     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, 0,
 303     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, 0,
 304     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 305     /* N = 1, K = 1...176 */
 306     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 307     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 308     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 309     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 310     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 311     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 312     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 313     /* N = 2, K = 2...176 */
 314     3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41,
 315     43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 65, 67, 69, 71, 73, 75, 77, 79,
 316     81, 83, 85, 87, 89, 91, 93, 95, 97, 99, 101, 103, 105, 107, 109, 111, 113,
 317     115, 117, 119, 121, 123, 125, 127, 129, 131, 133, 135, 137, 139, 141, 143,
 318     145, 147, 149, 151, 153, 155, 157, 159, 161, 163, 165, 167, 169, 171, 173,
 319     175, 177, 179, 181, 183, 185, 187, 189, 191, 193, 195, 197, 199, 201, 203,
 320     205, 207, 209, 211, 213, 215, 217, 219, 221, 223, 225, 227, 229, 231, 233,
 321     235, 237, 239, 241, 243, 245, 247, 249, 251, 253, 255, 257, 259, 261, 263,
 322     265, 267, 269, 271, 273, 275, 277, 279, 281, 283, 285, 287, 289, 291, 293,
 323     295, 297, 299, 301, 303, 305, 307, 309, 311, 313, 315, 317, 319, 321, 323,
 324     325, 327, 329, 331, 333, 335, 337, 339, 341, 343, 345, 347, 349, 351,
 325     /* N = 3, K = 3...176 */
 326     13, 25, 41, 61, 85, 113, 145, 181, 221, 265, 313, 365, 421, 481, 545, 613,
 327     685, 761, 841, 925, 1013, 1105, 1201, 1301, 1405, 1513, 1625, 1741, 1861,
 328     1985, 2113, 2245, 2381, 2521, 2665, 2813, 2965, 3121, 3281, 3445, 3613, 3785,
 329     3961, 4141, 4325, 4513, 4705, 4901, 5101, 5305, 5513, 5725, 5941, 6161, 6385,
 330     6613, 6845, 7081, 7321, 7565, 7813, 8065, 8321, 8581, 8845, 9113, 9385, 9661,
 331     9941, 10225, 10513, 10805, 11101, 11401, 11705, 12013, 12325, 12641, 12961,
 332     13285, 13613, 13945, 14281, 14621, 14965, 15313, 15665, 16021, 16381, 16745,
 333     17113, 17485, 17861, 18241, 18625, 19013, 19405, 19801, 20201, 20605, 21013,
 334     21425, 21841, 22261, 22685, 23113, 23545, 23981, 24421, 24865, 25313, 25765,
 335     26221, 26681, 27145, 27613, 28085, 28561, 29041, 29525, 30013, 30505, 31001,
 336     31501, 32005, 32513, 33025, 33541, 34061, 34585, 35113, 35645, 36181, 36721,
 337     37265, 37813, 38365, 38921, 39481, 40045, 40613, 41185, 41761, 42341, 42925,
 338     43513, 44105, 44701, 45301, 45905, 46513, 47125, 47741, 48361, 48985, 49613,
 339     50245, 50881, 51521, 52165, 52813, 53465, 54121, 54781, 55445, 56113, 56785,
 340     57461, 58141, 58825, 59513, 60205, 60901, 61601,
 341     /* N = 4, K = 4...176 */
 342     63, 129, 231, 377, 575, 833, 1159, 1561, 2047, 2625, 3303, 4089, 4991, 6017,
 343     7175, 8473, 9919, 11521, 13287, 15225, 17343, 19649, 22151, 24857, 27775,
 344     30913, 34279, 37881, 41727, 45825, 50183, 54809, 59711, 64897, 70375, 76153,
 345     82239, 88641, 95367, 102425, 109823, 117569, 125671, 134137, 142975, 152193,
 346     161799, 171801, 182207, 193025, 204263, 215929, 228031, 240577, 253575,
 347     267033, 280959, 295361, 310247, 325625, 341503, 357889, 374791, 392217,
 348     410175, 428673, 447719, 467321, 487487, 508225, 529543, 551449, 573951,
 349     597057, 620775, 645113, 670079, 695681, 721927, 748825, 776383, 804609,
 350     833511, 863097, 893375, 924353, 956039, 988441, 1021567, 1055425, 1090023,
 351     1125369, 1161471, 1198337, 1235975, 1274393, 1313599, 1353601, 1394407,
 352     1436025, 1478463, 1521729, 1565831, 1610777, 1656575, 1703233, 1750759,
 353     1799161, 1848447, 1898625, 1949703, 2001689, 2054591, 2108417, 2163175,
 354     2218873, 2275519, 2333121, 2391687, 2451225, 2511743, 2573249, 2635751,
 355     2699257, 2763775, 2829313, 2895879, 2963481, 3032127, 3101825, 3172583,
 356     3244409, 3317311, 3391297, 3466375, 3542553, 3619839, 3698241, 3777767,
 357     3858425, 3940223, 4023169, 4107271, 4192537, 4278975, 4366593, 4455399,
 358     4545401, 4636607, 4729025, 4822663, 4917529, 5013631, 5110977, 5209575,
 359     5309433, 5410559, 5512961, 5616647, 5721625, 5827903, 5935489, 6044391,
 360     6154617, 6266175, 6379073, 6493319, 6608921, 6725887, 6844225, 6963943,
 361     7085049, 7207551,
 362     /* N = 5, K = 5...176 */
 363     321, 681, 1289, 2241, 3649, 5641, 8361, 11969, 16641, 22569, 29961, 39041,
 364     50049, 63241, 78889, 97281, 118721, 143529, 172041, 204609, 241601, 283401,
 365     330409, 383041, 441729, 506921, 579081, 658689, 746241, 842249, 947241,
 366     1061761, 1186369, 1321641, 1468169, 1626561, 1797441, 1981449, 2179241,
 367     2391489, 2618881, 2862121, 3121929, 3399041, 3694209, 4008201, 4341801,
 368     4695809, 5071041, 5468329, 5888521, 6332481, 6801089, 7295241, 7815849,
 369     8363841, 8940161, 9545769, 10181641, 10848769, 11548161, 12280841, 13047849,
 370     13850241, 14689089, 15565481, 16480521, 17435329, 18431041, 19468809,
 371     20549801, 21675201, 22846209, 24064041, 25329929, 26645121, 28010881,
 372     29428489, 30899241, 32424449, 34005441, 35643561, 37340169, 39096641,
 373     40914369, 42794761, 44739241, 46749249, 48826241, 50971689, 53187081,
 374     55473921, 57833729, 60268041, 62778409, 65366401, 68033601, 70781609,
 375     73612041, 76526529, 79526721, 82614281, 85790889, 89058241, 92418049,
 376     95872041, 99421961, 103069569, 106816641, 110664969, 114616361, 118672641,
 377     122835649, 127107241, 131489289, 135983681, 140592321, 145317129, 150160041,
 378     155123009, 160208001, 165417001, 170752009, 176215041, 181808129, 187533321,
 379     193392681, 199388289, 205522241, 211796649, 218213641, 224775361, 231483969,
 380     238341641, 245350569, 252512961, 259831041, 267307049, 274943241, 282741889,
 381     290705281, 298835721, 307135529, 315607041, 324252609, 333074601, 342075401,
 382     351257409, 360623041, 370174729, 379914921, 389846081, 399970689, 410291241,
 383     420810249, 431530241, 442453761, 453583369, 464921641, 476471169, 488234561,
 384     500214441, 512413449, 524834241, 537479489, 550351881, 563454121, 576788929,
 385     590359041, 604167209, 618216201, 632508801,
 386     /* N = 6, K = 6...96 (technically V(109,5) fits in 32 bits, but that can't be
 387      achieved by splitting an Opus band) */
 388     1683, 3653, 7183, 13073, 22363, 36365, 56695, 85305, 124515, 177045, 246047,
 389     335137, 448427, 590557, 766727, 982729, 1244979, 1560549, 1937199, 2383409,
 390     2908411, 3522221, 4235671, 5060441, 6009091, 7095093, 8332863, 9737793,
 391     11326283, 13115773, 15124775, 17372905, 19880915, 22670725, 25765455,
 392     29189457, 32968347, 37129037, 41699767, 46710137, 52191139, 58175189,
 393     64696159, 71789409, 79491819, 87841821, 96879431, 106646281, 117185651,
 394     128542501, 140763503, 153897073, 167993403, 183104493, 199284183, 216588185,
 395     235074115, 254801525, 275831935, 298228865, 322057867, 347386557, 374284647,
 396     402823977, 433078547, 465124549, 499040399, 534906769, 572806619, 612825229,
 397     655050231, 699571641, 746481891, 795875861, 847850911, 902506913, 959946283,
 398     1020274013, 1083597703, 1150027593, 1219676595, 1292660325, 1369097135,
 399     1449108145, 1532817275, 1620351277, 1711839767, 1807415257, 1907213187,
 400     2011371957, 2120032959,
 401     /* N = 7, K = 7...54 (technically V(60,6) fits in 32 bits, but that can't be
 402      achieved by splitting an Opus band) */
 403     8989, 19825, 40081, 75517, 134245, 227305, 369305, 579125, 880685, 1303777,
 404     1884961, 2668525, 3707509, 5064793, 6814249, 9041957, 11847485, 15345233,
 405     19665841, 24957661, 31388293, 39146185, 48442297, 59511829, 72616013,
 406     88043969, 106114625, 127178701, 151620757, 179861305, 212358985, 249612805,
 407     292164445, 340600625, 395555537, 457713341, 527810725, 606639529, 695049433,
 408     793950709, 904317037, 1027188385, 1163673953, 1314955181, 1482288821,
 409     1667010073, 1870535785, 2094367717,
 410     /* N = 8, K = 8...37 (technically V(40,7) fits in 32 bits, but that can't be
 411      achieved by splitting an Opus band) */
 412     48639, 108545, 224143, 433905, 795455, 1392065, 2340495, 3800305, 5984767,
 413     9173505, 13726991, 20103025, 28875327, 40754369, 56610575, 77500017,
 414     104692735, 139703809, 184327311, 240673265, 311207743, 398796225, 506750351,
 415     638878193, 799538175, 993696769, 1226990095, 1505789553, 1837271615,
 416     2229491905,
 417     /* N = 9, K = 9...28 (technically V(29,8) fits in 32 bits, but that can't be
 418      achieved by splitting an Opus band) */
 419     265729, 598417, 1256465, 2485825, 4673345, 8405905, 14546705, 24331777,
 420     39490049, 62390545, 96220561, 145198913, 214828609, 312193553, 446304145,
 421     628496897, 872893441, 1196924561, 1621925137, 2173806145,
 422     /* N = 10, K = 10...24 */
 423     1462563, 3317445, 7059735, 14218905, 27298155, 50250765, 89129247, 152951073,
 424     254831667, 413442773, 654862247, 1014889769, 1541911931, 2300409629,
 425     3375210671,
 426     /* N = 11, K = 11...19 (technically V(20,10) fits in 32 bits, but that can't be
 427      achieved by splitting an Opus band) */
 428     8097453, 18474633, 39753273, 81270333, 158819253, 298199265, 540279585,
 429     948062325, 1616336765,
 430     /* N = 12, K = 12...18 */
 431     45046719, 103274625, 224298231, 464387817, 921406335, 1759885185,
 432     3248227095,
 433     /* N = 13, K = 13...16 */
 434     251595969, 579168825, 1267854873, 2653649025,
 435     /* N = 14, K = 14 */
 436     1409933619
 437 };
 438
 439 DECLARE_ALIGNED(32, static const float, celt_window)[120] = {
 440     6.7286966e-05f, 0.00060551348f, 0.0016815970f, 0.0032947962f, 0.0054439943f,
 441     0.0081276923f, 0.011344001f, 0.015090633f, 0.019364886f, 0.024163635f,
 442     0.029483315f, 0.035319905f, 0.041668911f, 0.048525347f, 0.055883718f,
 443     0.063737999f, 0.072081616f, 0.080907428f, 0.090207705f, 0.099974111f,
 444     0.11019769f, 0.12086883f, 0.13197729f, 0.14351214f, 0.15546177f,
 445     0.16781389f, 0.18055550f, 0.19367290f, 0.20715171f, 0.22097682f,
 446     0.23513243f, 0.24960208f, 0.26436860f, 0.27941419f, 0.29472040f,
 447     0.31026818f, 0.32603788f, 0.34200931f, 0.35816177f, 0.37447407f,
 448     0.39092462f, 0.40749142f, 0.42415215f, 0.44088423f, 0.45766484f,
 449     0.47447104f, 0.49127978f, 0.50806798f, 0.52481261f, 0.54149077f,
 450     0.55807973f, 0.57455701f, 0.59090049f, 0.60708841f, 0.62309951f,
 451     0.63891306f, 0.65450896f, 0.66986776f, 0.68497077f, 0.69980010f,
 452     0.71433873f, 0.72857055f, 0.74248043f, 0.75605424f, 0.76927895f,
 453     0.78214257f, 0.79463430f, 0.80674445f, 0.81846456f, 0.82978733f,
 454     0.84070669f, 0.85121779f, 0.86131698f, 0.87100183f, 0.88027111f,
 455     0.88912479f, 0.89756398f, 0.90559094f, 0.91320904f, 0.92042270f,
 456     0.92723738f, 0.93365955f, 0.93969656f, 0.94535671f, 0.95064907f,
 457     0.95558353f, 0.96017067f, 0.96442171f, 0.96834849f, 0.97196334f,
 458     0.97527906f, 0.97830883f, 0.98106616f, 0.98356480f, 0.98581869f,
 459     0.98784191f, 0.98964856f, 0.99125274f, 0.99266849f, 0.99390969f,
 460     0.99499004f, 0.99592297f, 0.99672162f, 0.99739874f, 0.99796667f,
 461     0.99843728f, 0.99882195f, 0.99913147f, 0.99937606f, 0.99956527f,
 462     0.99970802f, 0.99981248f, 0.99988613f, 0.99993565f, 0.99996697f,
 463     0.99998518f, 0.99999457f, 0.99999859f, 0.99999982f, 1.0000000f,
 464 };
 465
 466 /* square of the window, used for the postfilter */
 467 const float ff_celt_window2[120] = {
 468     4.5275357e-09f, 3.66647e-07f, 2.82777e-06f, 1.08557e-05f, 2.96371e-05f, 6.60594e-05f,
 469     0.000128686f, 0.000227727f, 0.000374999f, 0.000583881f, 0.000869266f, 0.0012475f,
 470     0.0017363f, 0.00235471f, 0.00312299f, 0.00406253f, 0.00519576f, 0.00654601f,
 471     0.00813743f, 0.00999482f, 0.0121435f, 0.0146093f, 0.017418f, 0.0205957f, 0.0241684f,
 472     0.0281615f, 0.0326003f, 0.0375092f, 0.0429118f, 0.0488308f, 0.0552873f, 0.0623012f,
 473     0.0698908f, 0.0780723f, 0.0868601f, 0.0962664f, 0.106301f, 0.11697f, 0.12828f,
 474     0.140231f, 0.152822f, 0.166049f, 0.179905f, 0.194379f, 0.209457f, 0.225123f, 0.241356f,
 475     0.258133f, 0.275428f, 0.293212f, 0.311453f, 0.330116f, 0.349163f, 0.368556f, 0.388253f,
 476     0.40821f, 0.428382f, 0.448723f, 0.469185f, 0.48972f, 0.51028f, 0.530815f, 0.551277f,
 477     0.571618f, 0.59179f, 0.611747f, 0.631444f, 0.650837f, 0.669884f, 0.688547f, 0.706788f,
 478     0.724572f, 0.741867f, 0.758644f, 0.774877f, 0.790543f, 0.805621f, 0.820095f, 0.833951f,
 479     0.847178f, 0.859769f, 0.87172f, 0.88303f, 0.893699f, 0.903734f, 0.91314f, 0.921928f,
 480     0.930109f, 0.937699f, 0.944713f, 0.951169f, 0.957088f, 0.962491f, 0.9674f, 0.971838f,
 481     0.975832f, 0.979404f, 0.982582f, 0.985391f, 0.987857f, 0.990005f, 0.991863f, 0.993454f,
 482     0.994804f, 0.995937f, 0.996877f, 0.997645f, 0.998264f, 0.998753f, 0.999131f, 0.999416f,
 483     0.999625f, 0.999772f, 0.999871f, 0.999934f, 0.99997f, 0.999989f, 0.999997f, 0.99999964f, 1.0f,
 484 };
 485
 486 static const uint32_t * const celt_pvq_u_row[15] = {
 487     celt_pvq_u +    0, celt_pvq_u +  176, celt_pvq_u +  351,
 488     celt_pvq_u +  525, celt_pvq_u +  698, celt_pvq_u +  870,
 489     celt_pvq_u + 1041, celt_pvq_u + 1131, celt_pvq_u + 1178,
 490     celt_pvq_u + 1207, celt_pvq_u + 1226, celt_pvq_u + 1240,
 491     celt_pvq_u + 1248, celt_pvq_u + 1254, celt_pvq_u + 1257
 492 };
 493
 494 static inline int16_t celt_cos(int16_t x)
 495 {
 496     x = (MUL16(x, x) + 4096) >> 13;
 497     x = (32767-x) + ROUND_MUL16(x, (-7651 + ROUND_MUL16(x, (8277 + ROUND_MUL16(-626, x)))));
 498     return 1+x;
 499 }
 500
 501 static inline int celt_log2tan(int isin, int icos)
 502 {
 503     int lc, ls;
 504     lc = opus_ilog(icos);
 505     ls = opus_ilog(isin);
 506     icos <<= 15 - lc;
 507     isin <<= 15 - ls;
 508     return (ls << 11) - (lc << 11) +
 509            ROUND_MUL16(isin, ROUND_MUL16(isin, -2597) + 7932) -
 510            ROUND_MUL16(icos, ROUND_MUL16(icos, -2597) + 7932);
 511 }
 512
 513 static inline uint32_t celt_rng(CeltContext *s)
 514 {
 515     s->seed = 1664525 * s->seed + 1013904223;
 516     return s->seed;
 517 }
 518
 519 static void celt_decode_coarse_energy(CeltContext *s, OpusRangeCoder *rc)
 520 {
 521     int i, j;
 522     float prev[2] = {0};
 523     float alpha, beta;
 524     const uint8_t *model;
 525
 526     /* use the 2D z-transform to apply prediction in both */
 527     /* the time domain (alpha) and the frequency domain (beta) */
 528
 529     if (opus_rc_tell(rc)+3 <= s->framebits && opus_rc_p2model(rc, 3)) {
 530         /* intra frame */
 531         alpha = 0;
 532         beta  = 1.0f - 4915.0f/32768.0f;
 533         model = celt_coarse_energy_dist[s->duration][1];
 534     } else {
 535         alpha = celt_alpha_coef[s->duration];
 536         beta  = 1.0f - celt_beta_coef[s->duration];
 537         model = celt_coarse_energy_dist[s->duration][0];
 538     }
 539
 540     for (i = 0; i < CELT_MAX_BANDS; i++) {
 541         for (j = 0; j < s->coded_channels; j++) {
 542             CeltFrame *frame = &s->frame[j];
 543             float value;
 544             int available;
 545
 546             if (i < s->startband || i >= s->endband) {
 547                 frame->energy[i] = 0.0;
 548                 continue;
 549             }
 550
 551             available = s->framebits - opus_rc_tell(rc);
 552             if (available >= 15) {
 553                 /* decode using a Laplace distribution */
 554                 int k = FFMIN(i, 20) << 1;
 555                 value = opus_rc_laplace(rc, model[k] << 7, model[k+1] << 6);
 556             } else if (available >= 2) {
 557                 int x = opus_rc_getsymbol(rc, celt_model_energy_small);
 558                 value = (x>>1) ^ -(x&1);
 559             } else if (available >= 1) {
 560                 value = -(float)opus_rc_p2model(rc, 1);
 561             } else value = -1;
 562
 563             frame->energy[i] = FFMAX(-9.0f, frame->energy[i]) * alpha + prev[j] + value;
 564             prev[j] += beta * value;
 565         }
 566     }
 567 }
 568
 569 static void celt_decode_fine_energy(CeltContext *s, OpusRangeCoder *rc)
 570 {
 571     int i;
 572     for (i = s->startband; i < s->endband; i++) {
 573         int j;
 574         if (!s->fine_bits[i])
 575             continue;
 576
 577         for (j = 0; j < s->coded_channels; j++) {
 578             CeltFrame *frame = &s->frame[j];
 579             int q2;
 580             float offset;
 581             q2 = opus_getrawbits(rc, s->fine_bits[i]);
 582             offset = (q2 + 0.5f) * (1 << (14 - s->fine_bits[i])) / 16384.0f - 0.5f;
 583             frame->energy[i] += offset;
 584         }
 585     }
 586 }
 587
 588 static void celt_decode_final_energy(CeltContext *s, OpusRangeCoder *rc,
 589                                      int bits_left)
 590 {
 591     int priority, i, j;
 592
 593     for (priority = 0; priority < 2; priority++) {
 594         for (i = s->startband; i < s->endband && bits_left >= s->coded_channels; i++) {
 595             if (s->fine_priority[i] != priority || s->fine_bits[i] >= CELT_MAX_FINE_BITS)
 596                 continue;
 597
 598             for (j = 0; j < s->coded_channels; j++) {
 599                 int q2;
 600                 float offset;
 601                 q2 = opus_getrawbits(rc, 1);
 602                 offset = (q2 - 0.5f) * (1 << (14 - s->fine_bits[i] - 1)) / 16384.0f;
 603                 s->frame[j].energy[i] += offset;
 604                 bits_left--;
 605             }
 606         }
 607     }
 608 }
 609
 610 static void celt_decode_tf_changes(CeltContext *s, OpusRangeCoder *rc,
 611                                    int transient)
 612 {
 613     int i, diff = 0, tf_select = 0, tf_changed = 0, tf_select_bit;
 614     int consumed, bits = transient ? 2 : 4;
 615
 616     consumed = opus_rc_tell(rc);
 617     tf_select_bit = (s->duration != 0 && consumed+bits+1 <= s->framebits);
 618
 619     for (i = s->startband; i < s->endband; i++) {
 620         if (consumed+bits+tf_select_bit <= s->framebits) {
 621             diff ^= opus_rc_p2model(rc, bits);
 622             consumed = opus_rc_tell(rc);
 623             tf_changed |= diff;
 624         }
 625         s->tf_change[i] = diff;
 626         bits = transient ? 4 : 5;
 627     }
 628
 629     if (tf_select_bit && celt_tf_select[s->duration][transient][0][tf_changed] !=
 630                          celt_tf_select[s->duration][transient][1][tf_changed])
 631         tf_select = opus_rc_p2model(rc, 1);
 632
 633     for (i = s->startband; i < s->endband; i++) {
 634         s->tf_change[i] = celt_tf_select[s->duration][transient][tf_select][s->tf_change[i]];
 635     }
 636 }
 637
 638 static void celt_decode_allocation(CeltContext *s, OpusRangeCoder *rc)
 639 {
 640     // approx. maximum bit allocation for each band before boost/trim
 641     int cap[CELT_MAX_BANDS];
 642     int boost[CELT_MAX_BANDS];
 643     int threshold[CELT_MAX_BANDS];
 644     int bits1[CELT_MAX_BANDS];
 645     int bits2[CELT_MAX_BANDS];
 646     int trim_offset[CELT_MAX_BANDS];
 647
 648     int skip_startband = s->startband;
 649     int dynalloc       = 6;
 650     int alloctrim      = 5;
 651     int extrabits      = 0;
 652
 653     int skip_bit            = 0;
 654     int intensitystereo_bit = 0;
 655     int dualstereo_bit      = 0;
 656
 657     int remaining, bandbits;
 658     int low, high, total, done;
 659     int totalbits;
 660     int consumed;
 661     int i, j;
 662
 663     consumed = opus_rc_tell(rc);
 664
 665     /* obtain spread flag */
 666     s->spread = CELT_SPREAD_NORMAL;
 667     if (consumed + 4 <= s->framebits)
 668         s->spread = opus_rc_getsymbol(rc, celt_model_spread);
 669
 670     /* generate static allocation caps */
 671     for (i = 0; i < CELT_MAX_BANDS; i++) {
 672         cap[i] = (celt_static_caps[s->duration][s->coded_channels - 1][i] + 64)
 673                  * celt_freq_range[i] << (s->coded_channels - 1) << s->duration >> 2;
 674     }
 675
 676     /* obtain band boost */
 677     totalbits = s->framebits << 3; // convert to 1/8 bits
 678     consumed = opus_rc_tell_frac(rc);
 679     for (i = s->startband; i < s->endband; i++) {
 680         int quanta, band_dynalloc;
 681
 682         boost[i] = 0;
 683
 684         quanta = celt_freq_range[i] << (s->coded_channels - 1) << s->duration;
 685         quanta = FFMIN(quanta << 3, FFMAX(6 << 3, quanta));
 686         band_dynalloc = dynalloc;
 687         while (consumed + (band_dynalloc<<3) < totalbits && boost[i] < cap[i]) {
 688             int add = opus_rc_p2model(rc, band_dynalloc);
 689             consumed = opus_rc_tell_frac(rc);
 690             if (!add)
 691                 break;
 692
 693             boost[i]     += quanta;
 694             totalbits    -= quanta;
 695             band_dynalloc = 1;
 696         }
 697         /* dynalloc is more likely to occur if it's already been used for earlier bands */
 698         if (boost[i])
 699             dynalloc = FFMAX(2, dynalloc - 1);
 700     }
 701
 702     /* obtain allocation trim */
 703     if (consumed + (6 << 3) <= totalbits)
 704         alloctrim = opus_rc_getsymbol(rc, celt_model_alloc_trim);
 705
 706     /* anti-collapse bit reservation */
 707     totalbits = (s->framebits << 3) - opus_rc_tell_frac(rc) - 1;
 708     s->anticollapse_bit = 0;
 709     if (s->blocks > 1 && s->duration >= 2 &&
 710         totalbits >= ((s->duration + 2) << 3))
 711         s->anticollapse_bit = 1 << 3;
 712     totalbits -= s->anticollapse_bit;
 713
 714     /* band skip bit reservation */
 715     if (totalbits >= 1 << 3)
 716         skip_bit = 1 << 3;
 717     totalbits -= skip_bit;
 718
 719     /* intensity/dual stereo bit reservation */
 720     if (s->coded_channels == 2) {
 721         intensitystereo_bit = celt_log2_frac[s->endband - s->startband];
 722         if (intensitystereo_bit <= totalbits) {
 723             totalbits -= intensitystereo_bit;
 724             if (totalbits >= 1 << 3) {
 725                 dualstereo_bit = 1 << 3;
 726                 totalbits -= 1 << 3;
 727             }
 728         } else
 729             intensitystereo_bit = 0;
 730     }
 731
 732     for (i = s->startband; i < s->endband; i++) {
 733         int trim     = alloctrim - 5 - s->duration;
 734         int band     = celt_freq_range[i] * (s->endband - i - 1);
 735         int duration = s->duration + 3;
 736         int scale    = duration + s->coded_channels - 1;
 737
 738         /* PVQ minimum allocation threshold, below this value the band is
 739          * skipped */
 740         threshold[i] = FFMAX(3 * celt_freq_range[i] << duration >> 4,
 741                              s->coded_channels << 3);
 742
 743         trim_offset[i] = trim * (band << scale) >> 6;
 744
 745         if (celt_freq_range[i] << s->duration == 1)
 746             trim_offset[i] -= s->coded_channels << 3;
 747     }
 748
 749     /* bisection */
 750     low  = 1;
 751     high = CELT_VECTORS - 1;
 752     while (low <= high) {
 753         int center = (low + high) >> 1;
 754         done = total = 0;
 755
 756         for (i = s->endband - 1; i >= s->startband; i--) {
 757             bandbits = celt_freq_range[i] * celt_static_alloc[center][i]
 758                        << (s->coded_channels - 1) << s->duration >> 2;
 759
 760             if (bandbits)
 761                 bandbits = FFMAX(0, bandbits + trim_offset[i]);
 762             bandbits += boost[i];
 763
 764             if (bandbits >= threshold[i] || done) {
 765                 done = 1;
 766                 total += FFMIN(bandbits, cap[i]);
 767             } else if (bandbits >= s->coded_channels << 3)
 768                 total += s->coded_channels << 3;
 769         }
 770
 771         if (total > totalbits)
 772             high = center - 1;
 773         else
 774             low = center + 1;
 775     }
 776     high = low--;
 777
 778     for (i = s->startband; i < s->endband; i++) {
 779         bits1[i] = celt_freq_range[i] * celt_static_alloc[low][i]
 780                    << (s->coded_channels - 1) << s->duration >> 2;
 781         bits2[i] = high >= CELT_VECTORS ? cap[i] :
 782                    celt_freq_range[i] * celt_static_alloc[high][i]
 783                    << (s->coded_channels - 1) << s->duration >> 2;
 784
 785         if (bits1[i])
 786             bits1[i] = FFMAX(0, bits1[i] + trim_offset[i]);
 787         if (bits2[i])
 788             bits2[i] = FFMAX(0, bits2[i] + trim_offset[i]);
 789         if (low)
 790             bits1[i] += boost[i];
 791         bits2[i] += boost[i];
 792
 793         if (boost[i])
 794             skip_startband = i;
 795         bits2[i] = FFMAX(0, bits2[i] - bits1[i]);
 796     }
 797
 798     /* bisection */
 799     low  = 0;
 800     high = 1 << CELT_ALLOC_STEPS;
 801     for (i = 0; i < CELT_ALLOC_STEPS; i++) {
 802         int center = (low + high) >> 1;
 803         done = total = 0;
 804
 805         for (j = s->endband - 1; j >= s->startband; j--) {
 806             bandbits = bits1[j] + (center * bits2[j] >> CELT_ALLOC_STEPS);
 807
 808             if (bandbits >= threshold[j] || done) {
 809                 done = 1;
 810                 total += FFMIN(bandbits, cap[j]);
 811             } else if (bandbits >= s->coded_channels << 3)
 812                 total += s->coded_channels << 3;
 813         }
 814         if (total > totalbits)
 815             high = center;
 816         else
 817             low = center;
 818     }
 819
 820     done = total = 0;
 821     for (i = s->endband - 1; i >= s->startband; i--) {
 822         bandbits = bits1[i] + (low * bits2[i] >> CELT_ALLOC_STEPS);
 823
 824         if (bandbits >= threshold[i] || done)
 825             done = 1;
 826         else
 827             bandbits = (bandbits >= s->coded_channels << 3) ?
 828                        s->coded_channels << 3 : 0;
 829
 830         bandbits     = FFMIN(bandbits, cap[i]);
 831         s->pulses[i] = bandbits;
 832         total      += bandbits;
 833     }
 834
 835     /* band skipping */
 836     for (s->codedbands = s->endband; ; s->codedbands--) {
 837         int allocation;
 838         j = s->codedbands - 1;
 839
 840         if (j == skip_startband) {
 841             /* all remaining bands are not skipped */
 842             totalbits += skip_bit;
 843             break;
 844         }
 845
 846         /* determine the number of bits available for coding "do not skip" markers */
 847         remaining   = totalbits - total;
 848         bandbits    = remaining / (celt_freq_bands[j+1] - celt_freq_bands[s->startband]);
 849         remaining  -= bandbits  * (celt_freq_bands[j+1] - celt_freq_bands[s->startband]);
 850         allocation  = s->pulses[j] + bandbits * celt_freq_range[j]
 851                       + FFMAX(0, remaining - (celt_freq_bands[j] - celt_freq_bands[s->startband]));
 852
 853         /* a "do not skip" marker is only coded if the allocation is
 854            above the chosen threshold */
 855         if (allocation >= FFMAX(threshold[j], (s->coded_channels + 1) <<3 )) {
 856             if (opus_rc_p2model(rc, 1))
 857                 break;
 858
 859             total      += 1 << 3;
 860             allocation -= 1 << 3;
 861         }
 862
 863         /* the band is skipped, so reclaim its bits */
 864         total -= s->pulses[j];
 865         if (intensitystereo_bit) {
 866             total -= intensitystereo_bit;
 867             intensitystereo_bit = celt_log2_frac[j - s->startband];
 868             total += intensitystereo_bit;
 869         }
 870
 871         total += s->pulses[j] = (allocation >= s->coded_channels << 3) ?
 872                               s->coded_channels << 3 : 0;
 873     }
 874
 875     /* obtain stereo flags */
 876     s->intensitystereo = 0;
 877     s->dualstereo      = 0;
 878     if (intensitystereo_bit)
 879         s->intensitystereo = s->startband +
 880                           opus_rc_unimodel(rc, s->codedbands + 1 - s->startband);
 881     if (s->intensitystereo <= s->startband)
 882         totalbits += dualstereo_bit; /* no intensity stereo means no dual stereo */
 883     else if (dualstereo_bit)
 884         s->dualstereo = opus_rc_p2model(rc, 1);
 885
 886     /* supply the remaining bits in this frame to lower bands */
 887     remaining = totalbits - total;
 888     bandbits  = remaining / (celt_freq_bands[s->codedbands] - celt_freq_bands[s->startband]);
 889     remaining -= bandbits * (celt_freq_bands[s->codedbands] - celt_freq_bands[s->startband]);
 890     for (i = s->startband; i < s->codedbands; i++) {
 891         int bits = FFMIN(remaining, celt_freq_range[i]);
 892
 893         s->pulses[i] += bits + bandbits * celt_freq_range[i];
 894         remaining    -= bits;
 895     }
 896
 897     for (i = s->startband; i < s->codedbands; i++) {
 898         int N = celt_freq_range[i] << s->duration;
 899         int prev_extra = extrabits;
 900         s->pulses[i] += extrabits;
 901
 902         if (N > 1) {
 903             int dof;        // degrees of freedom
 904             int temp;       // dof * channels * log(dof)
 905             int offset;     // fine energy quantization offset, i.e.
 906                             // extra bits assigned over the standard
 907                             // totalbits/dof
 908             int fine_bits, max_bits;
 909
 910             extrabits = FFMAX(0, s->pulses[i] - cap[i]);
 911             s->pulses[i] -= extrabits;
 912
 913             /* intensity stereo makes use of an extra degree of freedom */
 914             dof = N * s->coded_channels
 915                   + (s->coded_channels == 2 && N > 2 && !s->dualstereo && i < s->intensitystereo);
 916             temp = dof * (celt_log_freq_range[i] + (s->duration<<3));
 917             offset = (temp >> 1) - dof * CELT_FINE_OFFSET;
 918             if (N == 2) /* dof=2 is the only case that doesn't fit the model */
 919                 offset += dof<<1;
 920
 921             /* grant an additional bias for the first and second pulses */
 922             if (s->pulses[i] + offset < 2 * (dof << 3))
 923                 offset += temp >> 2;
 924             else if (s->pulses[i] + offset < 3 * (dof << 3))
 925                 offset += temp >> 3;
 926
 927             fine_bits = (s->pulses[i] + offset + (dof << 2)) / (dof << 3);
 928             max_bits  = FFMIN((s->pulses[i]>>3) >> (s->coded_channels - 1),
 929                               CELT_MAX_FINE_BITS);
 930
 931             max_bits  = FFMAX(max_bits, 0);
 932
 933             s->fine_bits[i] = av_clip(fine_bits, 0, max_bits);
 934
 935             /* if fine_bits was rounded down or capped,
 936                give priority for the final fine energy pass */
 937             s->fine_priority[i] = (s->fine_bits[i] * (dof<<3) >= s->pulses[i] + offset);
 938
 939             /* the remaining bits are assigned to PVQ */
 940             s->pulses[i] -= s->fine_bits[i] << (s->coded_channels - 1) << 3;
 941         } else {
 942             /* all bits go to fine energy except for the sign bit */
 943             extrabits = FFMAX(0, s->pulses[i] - (s->coded_channels << 3));
 944             s->pulses[i] -= extrabits;
 945             s->fine_bits[i] = 0;
 946             s->fine_priority[i] = 1;
 947         }
 948
 949         /* hand back a limited number of extra fine energy bits to this band */
 950         if (extrabits > 0) {
 951             int fineextra = FFMIN(extrabits >> (s->coded_channels + 2),
 952                                   CELT_MAX_FINE_BITS - s->fine_bits[i]);
 953             s->fine_bits[i] += fineextra;
 954
 955             fineextra <<= s->coded_channels + 2;
 956             s->fine_priority[i] = (fineextra >= extrabits - prev_extra);
 957             extrabits -= fineextra;
 958         }
 959     }
 960     s->remaining = extrabits;
 961
 962     /* skipped bands dedicate all of their bits for fine energy */
 963     for (; i < s->endband; i++) {
 964         s->fine_bits[i]     = s->pulses[i] >> (s->coded_channels - 1) >> 3;
 965         s->pulses[i]        = 0;
 966         s->fine_priority[i] = s->fine_bits[i] < 1;
 967     }
 968 }
 969
 970 static inline int celt_bits2pulses(const uint8_t *cache, int bits)
 971 {
 972     // TODO: Find the size of cache and make it into an array in the parameters list
 973     int i, low = 0, high;
 974
 975     high = cache[0];
 976     bits--;
 977
 978     for (i = 0; i < 6; i++) {
 979         int center = (low + high + 1) >> 1;
 980         if (cache[center] >= bits)
 981             high = center;
 982         else
 983             low = center;
 984     }
 985
 986     return (bits - (low == 0 ? -1 : cache[low]) <= cache[high] - bits) ? low : high;
 987 }
 988
 989 static inline int celt_pulses2bits(const uint8_t *cache, int pulses)
 990 {
 991     // TODO: Find the size of cache and make it into an array in the parameters list
 992    return (pulses == 0) ? 0 : cache[pulses] + 1;
 993 }
 994
 995 static inline void celt_normalize_residual(const int * av_restrict iy, float * av_restrict X,
 996                                            int N, float g)
 997 {
 998     int i;
 999     for (i = 0; i < N; i++)
1000         X[i] = g * iy[i];
1001 }
1002
1003 static void celt_exp_rotation1(float *X, unsigned int len, unsigned int stride,
1004                                float c, float s)
1005 {
1006     float *Xptr;
1007     int i;
1008
1009     Xptr = X;
1010     for (i = 0; i < len - stride; i++) {
1011         float x1, x2;
1012         x1           = Xptr[0];
1013         x2           = Xptr[stride];
1014         Xptr[stride] = c * x2 + s * x1;
1015         *Xptr++      = c * x1 - s * x2;
1016     }
1017
1018     Xptr = &X[len - 2 * stride - 1];
1019     for (i = len - 2 * stride - 1; i >= 0; i--) {
1020         float x1, x2;
1021         x1           = Xptr[0];
1022         x2           = Xptr[stride];
1023         Xptr[stride] = c * x2 + s * x1;
1024         *Xptr--      = c * x1 - s * x2;
1025     }
1026 }
1027
1028 static inline void celt_exp_rotation(float *X, unsigned int len,
1029                                      unsigned int stride, unsigned int K,
1030                                      enum CeltSpread spread)
1031 {
1032     unsigned int stride2 = 0;
1033     float c, s;
1034     float gain, theta;
1035     int i;
1036
1037     if (2*K >= len || spread == CELT_SPREAD_NONE)
1038         return;
1039
1040     gain = (float)len / (len + (20 - 5*spread) * K);
1041     theta = M_PI * gain * gain / 4;
1042
1043     c = cos(theta);
1044     s = sin(theta);
1045
1046     if (len >= stride << 3) {
1047         stride2 = 1;
1048         /* This is just a simple (equivalent) way of computing sqrt(len/stride) with rounding.
1049         It's basically incrementing long as (stride2+0.5)^2 < len/stride. */
1050         while ((stride2 * stride2 + stride2) * stride + (stride >> 2) < len)
1051             stride2++;
1052     }
1053
1054     /*NOTE: As a minor optimization, we could be passing around log2(B), not B, for both this and for
1055     extract_collapse_mask().*/
1056     len /= stride;
1057     for (i = 0; i < stride; i++) {
1058         if (stride2)
1059             celt_exp_rotation1(X + i * len, len, stride2, s, c);
1060         celt_exp_rotation1(X + i * len, len, 1, c, s);
1061     }
1062 }
1063
1064 static inline unsigned int celt_extract_collapse_mask(const int *iy,
1065                                                       unsigned int N,
1066                                                       unsigned int B)
1067 {
1068     unsigned int collapse_mask;
1069     int N0;
1070     int i, j;
1071
1072     if (B <= 1)
1073         return 1;
1074
1075     /*NOTE: As a minor optimization, we could be passing around log2(B), not B, for both this and for
1076     exp_rotation().*/
1077     N0 = N/B;
1078     collapse_mask = 0;
1079     for (i = 0; i < B; i++)
1080         for (j = 0; j < N0; j++)
1081             collapse_mask |= (iy[i*N0+j]!=0)<<i;
1082     return collapse_mask;
1083 }
1084
1085 static inline void celt_renormalize_vector(float *X, int N, float gain)
1086 {
1087     int i;
1088     float g = 1e-15f;
1089     for (i = 0; i < N; i++)
1090         g += X[i] * X[i];
1091     g = gain / sqrtf(g);
1092
1093     for (i = 0; i < N; i++)
1094         X[i] *= g;
1095 }
1096
1097 static inline void celt_stereo_merge(float *X, float *Y, float mid, int N)
1098 {
1099     int i;
1100     float xp = 0, side = 0;
1101     float E[2];
1102     float mid2;
1103     float t, gain[2];
1104
1105     /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
1106     for (i = 0; i < N; i++) {
1107         xp   += X[i] * Y[i];
1108         side += Y[i] * Y[i];
1109     }
1110
1111     /* Compensating for the mid normalization */
1112     xp *= mid;
1113     mid2 = mid;
1114     E[0] = mid2 * mid2 + side - 2 * xp;
1115     E[1] = mid2 * mid2 + side + 2 * xp;
1116     if (E[0] < 6e-4f || E[1] < 6e-4f) {
1117         for (i = 0; i < N; i++)
1118             Y[i] = X[i];
1119         return;
1120     }
1121
1122     t = E[0];
1123     gain[0] = 1.0f / sqrtf(t);
1124     t = E[1];
1125     gain[1] = 1.0f / sqrtf(t);
1126
1127     for (i = 0; i < N; i++) {
1128         float value[2];
1129         /* Apply mid scaling (side is already scaled) */
1130         value[0] = mid * X[i];
1131         value[1] = Y[i];
1132         X[i] = gain[0] * (value[0] - value[1]);
1133         Y[i] = gain[1] * (value[0] + value[1]);
1134     }
1135 }
1136
1137 static void celt_interleave_hadamard(float *tmp, float *X, int N0,
1138                                      int stride, int hadamard)
1139 {
1140     int i, j;
1141     int N = N0*stride;
1142
1143     if (hadamard) {
1144         const uint8_t *ordery = celt_hadamard_ordery + stride - 2;
1145         for (i = 0; i < stride; i++)
1146             for (j = 0; j < N0; j++)
1147                 tmp[j*stride+i] = X[ordery[i]*N0+j];
1148     } else {
1149         for (i = 0; i < stride; i++)
1150             for (j = 0; j < N0; j++)
1151                 tmp[j*stride+i] = X[i*N0+j];
1152     }
1153
1154     for (i = 0; i < N; i++)
1155         X[i] = tmp[i];
1156 }
1157
1158 static void celt_deinterleave_hadamard(float *tmp, float *X, int N0,
1159                                        int stride, int hadamard)
1160 {
1161     int i, j;
1162     int N = N0*stride;
1163
1164     if (hadamard) {
1165         const uint8_t *ordery = celt_hadamard_ordery + stride - 2;
1166         for (i = 0; i < stride; i++)
1167             for (j = 0; j < N0; j++)
1168                 tmp[ordery[i]*N0+j] = X[j*stride+i];
1169     } else {
1170         for (i = 0; i < stride; i++)
1171             for (j = 0; j < N0; j++)
1172                 tmp[i*N0+j] = X[j*stride+i];
1173     }
1174
1175     for (i = 0; i < N; i++)
1176         X[i] = tmp[i];
1177 }
1178
1179 static void celt_haar1(float *X, int N0, int stride)
1180 {
1181     int i, j;
1182     N0 >>= 1;
1183     for (i = 0; i < stride; i++) {
1184         for (j = 0; j < N0; j++) {
1185             float x0 = X[stride * (2 * j + 0) + i];
1186             float x1 = X[stride * (2 * j + 1) + i];
1187             X[stride * (2 * j + 0) + i] = (x0 + x1) * M_SQRT1_2;
1188             X[stride * (2 * j + 1) + i] = (x0 - x1) * M_SQRT1_2;
1189         }
1190     }
1191 }
1192
1193 static inline int celt_compute_qn(int N, int b, int offset, int pulse_cap,
1194                                   int dualstereo)
1195 {
1196     int qn, qb;
1197     int N2 = 2 * N - 1;
1198     if (dualstereo && N == 2)
1199         N2--;
1200
1201     /* The upper limit ensures that in a stereo split with itheta==16384, we'll
1202      * always have enough bits left over to code at least one pulse in the
1203      * side; otherwise it would collapse, since it doesn't get folded. */
1204     qb = FFMIN3(b - pulse_cap - (4 << 3), (b + N2 * offset) / N2, 8 << 3);
1205     qn = (qb < (1 << 3 >> 1)) ? 1 : ((celt_qn_exp2[qb & 0x7] >> (14 - (qb >> 3))) + 1) >> 1 << 1;
1206     return qn;
1207 }
1208
1209 // this code was adapted from libopus
1210 static inline uint64_t celt_cwrsi(unsigned int N, unsigned int K, unsigned int i, int *y)
1211 {
1212     uint64_t norm = 0;
1213     uint32_t p;
1214     int s, val;
1215     int k0;
1216
1217     while (N > 2) {
1218         uint32_t q;
1219
1220         /*Lots of pulses case:*/
1221         if (K >= N) {
1222             const uint32_t *row = celt_pvq_u_row[N];
1223
1224             /* Are the pulses in this dimension negative? */
1225             p  = row[K + 1];
1226             s  = -(i >= p);
1227             i -= p & s;
1228
1229             /*Count how many pulses were placed in this dimension.*/
1230             k0 = K;
1231             q = row[N];
1232             if (q > i) {
1233                 K = N;
1234                 do {
1235                     p = celt_pvq_u_row[--K][N];
1236                 } while (p > i);
1237             } else
1238                 for (p = row[K]; p > i; p = row[K])
1239                     K--;
1240
1241             i    -= p;
1242             val   = (k0 - K + s) ^ s;
1243             norm += val * val;
1244             *y++  = val;
1245         } else { /*Lots of dimensions case:*/
1246             /*Are there any pulses in this dimension at all?*/
1247             p = celt_pvq_u_row[K    ][N];
1248             q = celt_pvq_u_row[K + 1][N];
1249
1250             if (p <= i && i < q) {
1251                 i -= p;
1252                 *y++ = 0;
1253             } else {
1254                 /*Are the pulses in this dimension negative?*/
1255                 s  = -(i >= q);
1256                 i -= q & s;
1257
1258                 /*Count how many pulses were placed in this dimension.*/
1259                 k0 = K;
1260                 do p = celt_pvq_u_row[--K][N];
1261                 while (p > i);
1262
1263                 i    -= p;
1264                 val   = (k0 - K + s) ^ s;
1265                 norm += val * val;
1266                 *y++  = val;
1267             }
1268         }
1269         N--;
1270     }
1271
1272     /* N == 2 */
1273     p  = 2 * K + 1;
1274     s  = -(i >= p);
1275     i -= p & s;
1276     k0 = K;
1277     K  = (i + 1) / 2;
1278
1279     if (K)
1280         i -= 2 * K - 1;
1281
1282     val   = (k0 - K + s) ^ s;
1283     norm += val * val;
1284     *y++  = val;
1285
1286     /* N==1 */
1287     s     = -i;
1288     val   = (K + s) ^ s;
1289     norm += val * val;
1290     *y    = val;
1291
1292     return norm;
1293 }
1294
1295 static inline float celt_decode_pulses(OpusRangeCoder *rc, int *y, unsigned int N, unsigned int K)
1296 {
1297     unsigned int idx;
1298 #define CELT_PVQ_U(n, k) (celt_pvq_u_row[FFMIN(n, k)][FFMAX(n, k)])
1299 #define CELT_PVQ_V(n, k) (CELT_PVQ_U(n, k) + CELT_PVQ_U(n, (k) + 1))
1300     idx = opus_rc_unimodel(rc, CELT_PVQ_V(N, K));
1301     return celt_cwrsi(N, K, idx, y);
1302 }
1303
1304 /** Decode pulse vector and combine the result with the pitch vector to produce
1305     the final normalised signal in the current band. */
1306 static inline unsigned int celt_alg_unquant(OpusRangeCoder *rc, float *X,
1307                                             unsigned int N, unsigned int K,
1308                                             enum CeltSpread spread,
1309                                             unsigned int blocks, float gain)
1310 {
1311     int y[176];
1312
1313     gain /= sqrtf(celt_decode_pulses(rc, y, N, K));
1314     celt_normalize_residual(y, X, N, gain);
1315     celt_exp_rotation(X, N, blocks, K, spread);
1316     return celt_extract_collapse_mask(y, N, blocks);
1317 }
1318
1319 static unsigned int celt_decode_band(CeltContext *s, OpusRangeCoder *rc,
1320                                      const int band, float *X, float *Y,
1321                                      int N, int b, unsigned int blocks,
1322                                      float *lowband, int duration,
1323                                      float *lowband_out, int level,
1324                                      float gain, float *lowband_scratch,
1325                                      int fill)
1326 {
1327     const uint8_t *cache;
1328     int dualstereo, split;
1329     int imid = 0, iside = 0;
1330     unsigned int N0 = N;
1331     int N_B;
1332     int N_B0;
1333     int B0 = blocks;
1334     int time_divide = 0;
1335     int recombine = 0;
1336     int inv = 0;
1337     float mid = 0, side = 0;
1338     int longblocks = (B0 == 1);
1339     unsigned int cm = 0;
1340
1341     N_B0 = N_B = N / blocks;
1342     split = dualstereo = (Y != NULL);
1343
1344     if (N == 1) {
1345         /* special case for one sample */
1346         int i;
1347         float *x = X;
1348         for (i = 0; i <= dualstereo; i++) {
1349             int sign = 0;
1350             if (s->remaining2 >= 1<<3) {
1351                 sign           = opus_getrawbits(rc, 1);
1352                 s->remaining2 -= 1 << 3;
1353                 b             -= 1 << 3;
1354             }
1355             x[0] = sign ? -1.0f : 1.0f;
1356             x = Y;
1357         }
1358         if (lowband_out)
1359             lowband_out[0] = X[0];
1360         return 1;
1361     }
1362
1363     if (!dualstereo && level == 0) {
1364         int tf_change = s->tf_change[band];
1365         int k;
1366         if (tf_change > 0)
1367             recombine = tf_change;
1368         /* Band recombining to increase frequency resolution */
1369
1370         if (lowband &&
1371             (recombine || ((N_B & 1) == 0 && tf_change < 0) || B0 > 1)) {
1372             int j;
1373             for (j = 0; j < N; j++)
1374                 lowband_scratch[j] = lowband[j];
1375             lowband = lowband_scratch;
1376         }
1377
1378         for (k = 0; k < recombine; k++) {
1379             if (lowband)
1380                 celt_haar1(lowband, N >> k, 1 << k);
1381             fill = celt_bit_interleave[fill & 0xF] | celt_bit_interleave[fill >> 4] << 2;
1382         }
1383         blocks >>= recombine;
1384         N_B <<= recombine;
1385
1386         /* Increasing the time resolution */
1387         while ((N_B & 1) == 0 && tf_change < 0) {
1388             if (lowband)
1389                 celt_haar1(lowband, N_B, blocks);
1390             fill |= fill << blocks;
1391             blocks <<= 1;
1392             N_B >>= 1;
1393             time_divide++;
1394             tf_change++;
1395         }
1396         B0 = blocks;
1397         N_B0 = N_B;
1398
1399         /* Reorganize the samples in time order instead of frequency order */
1400         if (B0 > 1 && lowband)
1401             celt_deinterleave_hadamard(s->scratch, lowband, N_B >> recombine,
1402                                        B0 << recombine, longblocks);
1403     }
1404
1405     /* If we need 1.5 more bit than we can produce, split the band in two. */
1406     cache = celt_cache_bits +
1407             celt_cache_index[(duration + 1) * CELT_MAX_BANDS + band];
1408     if (!dualstereo && duration >= 0 && b > cache[cache[0]] + 12 && N > 2) {
1409         N >>= 1;
1410         Y = X + N;
1411         split = 1;
1412         duration -= 1;
1413         if (blocks == 1)
1414             fill = (fill & 1) | (fill << 1);
1415         blocks = (blocks + 1) >> 1;
1416     }
1417
1418     if (split) {
1419         int qn;
1420         int itheta = 0;
1421         int mbits, sbits, delta;
1422         int qalloc;
1423         int pulse_cap;
1424         int offset;
1425         int orig_fill;
1426         int tell;
1427
1428         /* Decide on the resolution to give to the split parameter theta */
1429         pulse_cap = celt_log_freq_range[band] + duration * 8;
1430         offset = (pulse_cap >> 1) - (dualstereo && N == 2 ? CELT_QTHETA_OFFSET_TWOPHASE :
1431                                                           CELT_QTHETA_OFFSET);
1432         qn = (dualstereo && band >= s->intensitystereo) ? 1 :
1433              celt_compute_qn(N, b, offset, pulse_cap, dualstereo);
1434         tell = opus_rc_tell_frac(rc);
1435         if (qn != 1) {
1436             /* Entropy coding of the angle. We use a uniform pdf for the
1437             time split, a step for stereo, and a triangular one for the rest. */
1438             if (dualstereo && N > 2)
1439                 itheta = opus_rc_stepmodel(rc, qn/2);
1440             else if (dualstereo || B0 > 1)
1441                 itheta = opus_rc_unimodel(rc, qn+1);
1442             else
1443                 itheta = opus_rc_trimodel(rc, qn);
1444             itheta = itheta * 16384 / qn;
1445             /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.
1446             Let's do that at higher complexity */
1447         } else if (dualstereo) {
1448             inv = (b > 2 << 3 && s->remaining2 > 2 << 3) ? opus_rc_p2model(rc, 2) : 0;
1449             itheta = 0;
1450         }
1451         qalloc = opus_rc_tell_frac(rc) - tell;
1452         b -= qalloc;
1453
1454         orig_fill = fill;
1455         if (itheta == 0) {
1456             imid = 32767;
1457             iside = 0;
1458             fill = av_mod_uintp2(fill, blocks);
1459             delta = -16384;
1460         } else if (itheta == 16384) {
1461             imid = 0;
1462             iside = 32767;
1463             fill &= ((1 << blocks) - 1) << blocks;
1464             delta = 16384;
1465         } else {
1466             imid = celt_cos(itheta);
1467             iside = celt_cos(16384-itheta);
1468             /* This is the mid vs side allocation that minimizes squared error
1469             in that band. */
1470             delta = ROUND_MUL16((N - 1) << 7, celt_log2tan(iside, imid));
1471         }
1472
1473         mid  = imid  / 32768.0f;
1474         side = iside / 32768.0f;
1475
1476         /* This is a special case for N=2 that only works for stereo and takes
1477         advantage of the fact that mid and side are orthogonal to encode
1478         the side with just one bit. */
1479         if (N == 2 && dualstereo) {
1480             int c;
1481             int sign = 0;
1482             float tmp;
1483             float *x2, *y2;
1484             mbits = b;
1485             /* Only need one bit for the side */
1486             sbits = (itheta != 0 && itheta != 16384) ? 1 << 3 : 0;
1487             mbits -= sbits;
1488             c = (itheta > 8192);
1489             s->remaining2 -= qalloc+sbits;
1490
1491             x2 = c ? Y : X;
1492             y2 = c ? X : Y;
1493             if (sbits)
1494                 sign = opus_getrawbits(rc, 1);
1495             sign = 1 - 2 * sign;
1496             /* We use orig_fill here because we want to fold the side, but if
1497             itheta==16384, we'll have cleared the low bits of fill. */
1498             cm = celt_decode_band(s, rc, band, x2, NULL, N, mbits, blocks,
1499                                   lowband, duration, lowband_out, level, gain,
1500                                   lowband_scratch, orig_fill);
1501             /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
1502             and there's no need to worry about mixing with the other channel. */
1503             y2[0] = -sign * x2[1];
1504             y2[1] =  sign * x2[0];
1505             X[0] *= mid;
1506             X[1] *= mid;
1507             Y[0] *= side;
1508             Y[1] *= side;
1509             tmp = X[0];
1510             X[0] = tmp - Y[0];
1511             Y[0] = tmp + Y[0];
1512             tmp = X[1];
1513             X[1] = tmp - Y[1];
1514             Y[1] = tmp + Y[1];
1515         } else {
1516             /* "Normal" split code */
1517             float *next_lowband2     = NULL;
1518             float *next_lowband_out1 = NULL;
1519             int next_level = 0;
1520             int rebalance;
1521
1522             /* Give more bits to low-energy MDCTs than they would
1523              * otherwise deserve */
1524             if (B0 > 1 && !dualstereo && (itheta & 0x3fff)) {
1525                 if (itheta > 8192)
1526                     /* Rough approximation for pre-echo masking */
1527                     delta -= delta >> (4 - duration);
1528                 else
1529                     /* Corresponds to a forward-masking slope of
1530                      * 1.5 dB per 10 ms */
1531                     delta = FFMIN(0, delta + (N << 3 >> (5 - duration)));
1532             }
1533             mbits = av_clip((b - delta) / 2, 0, b);
1534             sbits = b - mbits;
1535             s->remaining2 -= qalloc;
1536
1537             if (lowband && !dualstereo)
1538                 next_lowband2 = lowband + N; /* >32-bit split case */
1539
1540             /* Only stereo needs to pass on lowband_out.
1541              * Otherwise, it's handled at the end */
1542             if (dualstereo)
1543                 next_lowband_out1 = lowband_out;
1544             else
1545                 next_level = level + 1;
1546
1547             rebalance = s->remaining2;
1548             if (mbits >= sbits) {
1549                 /* In stereo mode, we do not apply a scaling to the mid
1550                  * because we need the normalized mid for folding later */
1551                 cm = celt_decode_band(s, rc, band, X, NULL, N, mbits, blocks,
1552                                       lowband, duration, next_lowband_out1,
1553                                       next_level, dualstereo ? 1.0f : (gain * mid),
1554                                       lowband_scratch, fill);
1555
1556                 rebalance = mbits - (rebalance - s->remaining2);
1557                 if (rebalance > 3 << 3 && itheta != 0)
1558                     sbits += rebalance - (3 << 3);
1559
1560                 /* For a stereo split, the high bits of fill are always zero,
1561                  * so no folding will be done to the side. */
1562                 cm |= celt_decode_band(s, rc, band, Y, NULL, N, sbits, blocks,
1563                                        next_lowband2, duration, NULL,
1564                                        next_level, gain * side, NULL,
1565                                        fill >> blocks) << ((B0 >> 1) & (dualstereo - 1));
1566             } else {
1567                 /* For a stereo split, the high bits of fill are always zero,
1568                  * so no folding will be done to the side. */
1569                 cm = celt_decode_band(s, rc, band, Y, NULL, N, sbits, blocks,
1570                                       next_lowband2, duration, NULL,
1571                                       next_level, gain * side, NULL,
1572                                       fill >> blocks) << ((B0 >> 1) & (dualstereo - 1));
1573
1574                 rebalance = sbits - (rebalance - s->remaining2);
1575                 if (rebalance > 3 << 3 && itheta != 16384)
1576                     mbits += rebalance - (3 << 3);
1577
1578                 /* In stereo mode, we do not apply a scaling to the mid because
1579                  * we need the normalized mid for folding later */
1580                 cm |= celt_decode_band(s, rc, band, X, NULL, N, mbits, blocks,
1581                                        lowband, duration, next_lowband_out1,
1582                                        next_level, dualstereo ? 1.0f : (gain * mid),
1583                                        lowband_scratch, fill);
1584             }
1585         }
1586     } else {
1587         /* This is the basic no-split case */
1588         unsigned int q         = celt_bits2pulses(cache, b);
1589         unsigned int curr_bits = celt_pulses2bits(cache, q);
1590         s->remaining2 -= curr_bits;
1591
1592         /* Ensures we can never bust the budget */
1593         while (s->remaining2 < 0 && q > 0) {
1594             s->remaining2 += curr_bits;
1595             curr_bits      = celt_pulses2bits(cache, --q);
1596             s->remaining2 -= curr_bits;
1597         }
1598
1599         if (q != 0) {
1600             /* Finally do the actual quantization */
1601             cm = celt_alg_unquant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
1602                                   s->spread, blocks, gain);
1603         } else {
1604             /* If there's no pulse, fill the band anyway */
1605             int j;
1606             unsigned int cm_mask = (1 << blocks) - 1;
1607             fill &= cm_mask;
1608             if (!fill) {
1609                 for (j = 0; j < N; j++)
1610                     X[j] = 0.0f;
1611             } else {
1612                 if (!lowband) {
1613                     /* Noise */
1614                     for (j = 0; j < N; j++)
1615                         X[j] = (((int32_t)celt_rng(s)) >> 20);
1616                     cm = cm_mask;
1617                 } else {
1618                     /* Folded spectrum */
1619                     for (j = 0; j < N; j++) {
1620                         /* About 48 dB below the "normal" folding level */
1621                         X[j] = lowband[j] + (((celt_rng(s)) & 0x8000) ? 1.0f / 256 : -1.0f / 256);
1622                     }
1623                     cm = fill;
1624                 }
1625                 celt_renormalize_vector(X, N, gain);
1626             }
1627         }
1628     }
1629
1630     /* This code is used by the decoder and by the resynthesis-enabled encoder */
1631     if (dualstereo) {
1632         int j;
1633         if (N != 2)
1634             celt_stereo_merge(X, Y, mid, N);
1635         if (inv) {
1636             for (j = 0; j < N; j++)
1637                 Y[j] *= -1;
1638         }
1639     } else if (level == 0) {
1640         int k;
1641
1642         /* Undo the sample reorganization going from time order to frequency order */
1643         if (B0 > 1)
1644             celt_interleave_hadamard(s->scratch, X, N_B>>recombine,
1645                                      B0<<recombine, longblocks);
1646
1647         /* Undo time-freq changes that we did earlier */
1648         N_B = N_B0;
1649         blocks = B0;
1650         for (k = 0; k < time_divide; k++) {
1651             blocks >>= 1;
1652             N_B <<= 1;
1653             cm |= cm >> blocks;
1654             celt_haar1(X, N_B, blocks);
1655         }
1656
1657         for (k = 0; k < recombine; k++) {
1658             cm = celt_bit_deinterleave[cm];
1659             celt_haar1(X, N0>>k, 1<<k);
1660         }
1661         blocks <<= recombine;
1662
1663         /* Scale output for later folding */
1664         if (lowband_out) {
1665             int j;
1666             float n = sqrtf(N0);
1667             for (j = 0; j < N0; j++)
1668                 lowband_out[j] = n * X[j];
1669         }
1670         cm = av_mod_uintp2(cm, blocks);
1671     }
1672     return cm;
1673 }
1674
1675 static void celt_denormalize(CeltContext *s, CeltFrame *frame, float *data)
1676 {
1677     int i, j;
1678
1679     for (i = s->startband; i < s->endband; i++) {
1680         float *dst = data + (celt_freq_bands[i] << s->duration);
1681         float norm = exp2(frame->energy[i] + celt_mean_energy[i]);
1682
1683         for (j = 0; j < celt_freq_range[i] << s->duration; j++)
1684             dst[j] *= norm;
1685     }
1686 }
1687
1688 static void celt_postfilter_apply_transition(CeltFrame *frame, float *data)
1689 {
1690     const int T0 = frame->pf_period_old;
1691     const int T1 = frame->pf_period;
1692
1693     float g00, g01, g02;
1694     float g10, g11, g12;
1695
1696     float x0, x1, x2, x3, x4;
1697
1698     int i;
1699
1700     if (frame->pf_gains[0]     == 0.0 &&
1701         frame->pf_gains_old[0] == 0.0)
1702         return;
1703
1704     g00 = frame->pf_gains_old[0];
1705     g01 = frame->pf_gains_old[1];
1706     g02 = frame->pf_gains_old[2];
1707     g10 = frame->pf_gains[0];
1708     g11 = frame->pf_gains[1];
1709     g12 = frame->pf_gains[2];
1710
1711     x1 = data[-T1 + 1];
1712     x2 = data[-T1];
1713     x3 = data[-T1 - 1];
1714     x4 = data[-T1 - 2];
1715
1716     for (i = 0; i < CELT_OVERLAP; i++) {
1717         float w = ff_celt_window2[i];
1718         x0 = data[i - T1 + 2];
1719
1720         data[i] +=  (1.0 - w) * g00 * data[i - T0]                          +
1721                     (1.0 - w) * g01 * (data[i - T0 - 1] + data[i - T0 + 1]) +
1722                     (1.0 - w) * g02 * (data[i - T0 - 2] + data[i - T0 + 2]) +
1723                     w         * g10 * x2                                    +
1724                     w         * g11 * (x1 + x3)                             +
1725                     w         * g12 * (x0 + x4);
1726         x4 = x3;
1727         x3 = x2;
1728         x2 = x1;
1729         x1 = x0;
1730     }
1731 }
1732
1733 static void celt_postfilter_apply(CeltFrame *frame,
1734                                   float *data, int len)
1735 {
1736     const int T = frame->pf_period;
1737     float g0, g1, g2;
1738     float x0, x1, x2, x3, x4;
1739     int i;
1740
1741     if (frame->pf_gains[0] == 0.0 || len <= 0)
1742         return;
1743
1744     g0 = frame->pf_gains[0];
1745     g1 = frame->pf_gains[1];
1746     g2 = frame->pf_gains[2];
1747
1748     x4 = data[-T - 2];
1749     x3 = data[-T - 1];
1750     x2 = data[-T];
1751     x1 = data[-T + 1];
1752
1753     for (i = 0; i < len; i++) {
1754         x0 = data[i - T + 2];
1755         data[i] += g0 * x2        +
1756                    g1 * (x1 + x3) +
1757                    g2 * (x0 + x4);
1758         x4 = x3;
1759         x3 = x2;
1760         x2 = x1;
1761         x1 = x0;
1762     }
1763 }
1764
1765 static void celt_postfilter(CeltContext *s, CeltFrame *frame)
1766 {
1767     int len = s->blocksize * s->blocks;
1768
1769     celt_postfilter_apply_transition(frame, frame->buf + 1024);
1770
1771     frame->pf_period_old = frame->pf_period;
1772     memcpy(frame->pf_gains_old, frame->pf_gains, sizeof(frame->pf_gains));
1773
1774     frame->pf_period = frame->pf_period_new;
1775     memcpy(frame->pf_gains, frame->pf_gains_new, sizeof(frame->pf_gains));
1776
1777     if (len > CELT_OVERLAP) {
1778         celt_postfilter_apply_transition(frame, frame->buf + 1024 + CELT_OVERLAP);
1779         celt_postfilter_apply(frame, frame->buf + 1024 + 2 * CELT_OVERLAP,
1780                               len - 2 * CELT_OVERLAP);
1781
1782         frame->pf_period_old = frame->pf_period;
1783         memcpy(frame->pf_gains_old, frame->pf_gains, sizeof(frame->pf_gains));
1784     }
1785
1786     memmove(frame->buf, frame->buf + len, (1024 + CELT_OVERLAP / 2) * sizeof(float));
1787 }
1788
1789 static int parse_postfilter(CeltContext *s, OpusRangeCoder *rc, int consumed)
1790 {
1791     static const float postfilter_taps[3][3] = {
1792         { 0.3066406250f, 0.2170410156f, 0.1296386719f },
1793         { 0.4638671875f, 0.2680664062f, 0.0           },
1794         { 0.7998046875f, 0.1000976562f, 0.0           }
1795     };
1796     int i;
1797
1798     memset(s->frame[0].pf_gains_new, 0, sizeof(s->frame[0].pf_gains_new));
1799     memset(s->frame[1].pf_gains_new, 0, sizeof(s->frame[1].pf_gains_new));
1800
1801     if (s->startband == 0 && consumed + 16 <= s->framebits) {
1802         int has_postfilter = opus_rc_p2model(rc, 1);
1803         if (has_postfilter) {
1804             float gain;
1805             int tapset, octave, period;
1806
1807             octave = opus_rc_unimodel(rc, 6);
1808             period = (16 << octave) + opus_getrawbits(rc, 4 + octave) - 1;
1809             gain   = 0.09375f * (opus_getrawbits(rc, 3) + 1);
1810             tapset = (opus_rc_tell(rc) + 2 <= s->framebits) ?
1811                      opus_rc_getsymbol(rc, celt_model_tapset) : 0;
1812
1813             for (i = 0; i < 2; i++) {
1814                 CeltFrame *frame = &s->frame[i];
1815
1816                 frame->pf_period_new = FFMAX(period, CELT_POSTFILTER_MINPERIOD);
1817                 frame->pf_gains_new[0] = gain * postfilter_taps[tapset][0];
1818                 frame->pf_gains_new[1] = gain * postfilter_taps[tapset][1];
1819                 frame->pf_gains_new[2] = gain * postfilter_taps[tapset][2];
1820             }
1821         }
1822
1823         consumed = opus_rc_tell(rc);
1824     }
1825
1826     return consumed;
1827 }
1828
1829 static void process_anticollapse(CeltContext *s, CeltFrame *frame, float *X)
1830 {
1831     int i, j, k;
1832
1833     for (i = s->startband; i < s->endband; i++) {
1834         int renormalize = 0;
1835         float *xptr;
1836         float prev[2];
1837         float Ediff, r;
1838         float thresh, sqrt_1;
1839         int depth;
1840
1841         /* depth in 1/8 bits */
1842         depth = (1 + s->pulses[i]) / (celt_freq_range[i] << s->duration);
1843         thresh = exp2f(-1.0 - 0.125f * depth);
1844         sqrt_1 = 1.0f / sqrtf(celt_freq_range[i] << s->duration);
1845
1846         xptr = X + (celt_freq_bands[i] << s->duration);
1847
1848         prev[0] = frame->prev_energy[0][i];
1849         prev[1] = frame->prev_energy[1][i];
1850         if (s->coded_channels == 1) {
1851             CeltFrame *frame1 = &s->frame[1];
1852
1853             prev[0] = FFMAX(prev[0], frame1->prev_energy[0][i]);
1854             prev[1] = FFMAX(prev[1], frame1->prev_energy[1][i]);
1855         }
1856         Ediff = frame->energy[i] - FFMIN(prev[0], prev[1]);
1857         Ediff = FFMAX(0, Ediff);
1858
1859         /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because
1860         short blocks don't have the same energy as long */
1861         r = exp2(1 - Ediff);
1862         if (s->duration == 3)
1863             r *= M_SQRT2;
1864         r = FFMIN(thresh, r) * sqrt_1;
1865         for (k = 0; k < 1 << s->duration; k++) {
1866             /* Detect collapse */
1867             if (!(frame->collapse_masks[i] & 1 << k)) {
1868                 /* Fill with noise */
1869                 for (j = 0; j < celt_freq_range[i]; j++)
1870                     xptr[(j << s->duration) + k] = (celt_rng(s) & 0x8000) ? r : -r;
1871                 renormalize = 1;
1872             }
1873         }
1874
1875         /* We just added some energy, so we need to renormalize */
1876         if (renormalize)
1877             celt_renormalize_vector(xptr, celt_freq_range[i] << s->duration, 1.0f);
1878     }
1879 }
1880
1881 static void celt_decode_bands(CeltContext *s, OpusRangeCoder *rc)
1882 {
1883     float lowband_scratch[8 * 22];
1884     float norm[2 * 8 * 100];
1885
1886     int totalbits = (s->framebits << 3) - s->anticollapse_bit;
1887
1888     int update_lowband = 1;
1889     int lowband_offset = 0;
1890
1891     int i, j;
1892
1893     memset(s->coeffs, 0, sizeof(s->coeffs));
1894
1895     for (i = s->startband; i < s->endband; i++) {
1896         int band_offset = celt_freq_bands[i] << s->duration;
1897         int band_size   = celt_freq_range[i] << s->duration;
1898         float *X = s->coeffs[0] + band_offset;
1899         float *Y = (s->coded_channels == 2) ? s->coeffs[1] + band_offset : NULL;
1900
1901         int consumed = opus_rc_tell_frac(rc);
1902         float *norm2 = norm + 8 * 100;
1903         int effective_lowband = -1;
1904         unsigned int cm[2];
1905         int b;
1906
1907         /* Compute how many bits we want to allocate to this band */
1908         if (i != s->startband)
1909             s->remaining -= consumed;
1910         s->remaining2 = totalbits - consumed - 1;
1911         if (i <= s->codedbands - 1) {
1912             int curr_balance = s->remaining / FFMIN(3, s->codedbands-i);
1913             b = av_clip_uintp2(FFMIN(s->remaining2 + 1, s->pulses[i] + curr_balance), 14);
1914         } else
1915             b = 0;
1916
1917         if (celt_freq_bands[i] - celt_freq_range[i] >= celt_freq_bands[s->startband] &&
1918             (update_lowband || lowband_offset == 0))
1919             lowband_offset = i;
1920
1921         /* Get a conservative estimate of the collapse_mask's for the bands we're
1922         going to be folding from. */
1923         if (lowband_offset != 0 && (s->spread != CELT_SPREAD_AGGRESSIVE ||
1924                                     s->blocks > 1 || s->tf_change[i] < 0)) {
1925             int foldstart, foldend;
1926
1927             /* This ensures we never repeat spectral content within one band */
1928             effective_lowband = FFMAX(celt_freq_bands[s->startband],
1929                                       celt_freq_bands[lowband_offset] - celt_freq_range[i]);
1930             foldstart = lowband_offset;
1931             while (celt_freq_bands[--foldstart] > effective_lowband);
1932             foldend = lowband_offset - 1;
1933             while (celt_freq_bands[++foldend] < effective_lowband + celt_freq_range[i]);
1934
1935             cm[0] = cm[1] = 0;
1936             for (j = foldstart; j < foldend; j++) {
1937                 cm[0] |= s->frame[0].collapse_masks[j];
1938                 cm[1] |= s->frame[s->coded_channels - 1].collapse_masks[j];
1939             }
1940         } else
1941             /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1942             always) be non-zero.*/
1943             cm[0] = cm[1] = (1 << s->blocks) - 1;
1944
1945         if (s->dualstereo && i == s->intensitystereo) {
1946             /* Switch off dual stereo to do intensity */
1947             s->dualstereo = 0;
1948             for (j = celt_freq_bands[s->startband] << s->duration; j < band_offset; j++)
1949                 norm[j] = (norm[j] + norm2[j]) / 2;
1950         }
1951
1952         if (s->dualstereo) {
1953             cm[0] = celt_decode_band(s, rc, i, X, NULL, band_size, b / 2, s->blocks,
1954                                      effective_lowband != -1 ? norm + (effective_lowband << s->duration) : NULL, s->duration,
1955             norm + band_offset, 0, 1.0f, lowband_scratch, cm[0]);
1956
1957             cm[1] = celt_decode_band(s, rc, i, Y, NULL, band_size, b/2, s->blocks,
1958                                      effective_lowband != -1 ? norm2 + (effective_lowband << s->duration) : NULL, s->duration,
1959             norm2 + band_offset, 0, 1.0f, lowband_scratch, cm[1]);
1960         } else {
1961             cm[0] = celt_decode_band(s, rc, i, X, Y, band_size, b, s->blocks,
1962             effective_lowband != -1 ? norm + (effective_lowband << s->duration) : NULL, s->duration,
1963             norm + band_offset, 0, 1.0f, lowband_scratch, cm[0]|cm[1]);
1964
1965             cm[1] = cm[0];
1966         }
1967
1968         s->frame[0].collapse_masks[i]                     = (uint8_t)cm[0];
1969         s->frame[s->coded_channels - 1].collapse_masks[i] = (uint8_t)cm[1];
1970         s->remaining += s->pulses[i] + consumed;
1971
1972         /* Update the folding position only as long as we have 1 bit/sample depth */
1973         update_lowband = (b > band_size << 3);
1974     }
1975 }
1976
1977 int ff_celt_decode_frame(CeltContext *s, OpusRangeCoder *rc,
1978                          float **output, int coded_channels, int frame_size,
1979                          int startband,  int endband)
1980 {
1981     int i, j;
1982
1983     int consumed;           // bits of entropy consumed thus far for this frame
1984     int silence = 0;
1985     int transient = 0;
1986     int anticollapse = 0;
1987     IMDCT15Context *imdct;
1988     float imdct_scale = 1.0;
1989
1990     if (coded_channels != 1 && coded_channels != 2) {
1991         av_log(s->avctx, AV_LOG_ERROR, "Invalid number of coded channels: %d\n",
1992                coded_channels);
1993         return AVERROR_INVALIDDATA;
1994     }
1995     if (startband < 0 || startband > endband || endband > CELT_MAX_BANDS) {
1996         av_log(s->avctx, AV_LOG_ERROR, "Invalid start/end band: %d %d\n",
1997                startband, endband);
1998         return AVERROR_INVALIDDATA;
1999     }
2000
2001     s->flushed        = 0;
2002     s->coded_channels = coded_channels;
2003     s->startband      = startband;
2004     s->endband        = endband;
2005     s->framebits      = rc->rb.bytes * 8;
2006
2007     s->duration = av_log2(frame_size / CELT_SHORT_BLOCKSIZE);
2008     if (s->duration > CELT_MAX_LOG_BLOCKS ||
2009         frame_size != CELT_SHORT_BLOCKSIZE * (1 << s->duration)) {
2010         av_log(s->avctx, AV_LOG_ERROR, "Invalid CELT frame size: %d\n",
2011                frame_size);
2012         return AVERROR_INVALIDDATA;
2013     }
2014
2015     if (!s->output_channels)
2016         s->output_channels = coded_channels;
2017
2018     memset(s->frame[0].collapse_masks, 0, sizeof(s->frame[0].collapse_masks));
2019     memset(s->frame[1].collapse_masks, 0, sizeof(s->frame[1].collapse_masks));
2020
2021     consumed = opus_rc_tell(rc);
2022
2023     /* obtain silence flag */
2024     if (consumed >= s->framebits)
2025         silence = 1;
2026     else if (consumed == 1)
2027         silence = opus_rc_p2model(rc, 15);
2028
2029
2030     if (silence) {
2031         consumed = s->framebits;
2032         rc->total_read_bits += s->framebits - opus_rc_tell(rc);
2033     }
2034
2035     /* obtain post-filter options */
2036     consumed = parse_postfilter(s, rc, consumed);
2037
2038     /* obtain transient flag */
2039     if (s->duration != 0 && consumed+3 <= s->framebits)
2040         transient = opus_rc_p2model(rc, 3);
2041
2042     s->blocks    = transient ? 1 << s->duration : 1;
2043     s->blocksize = frame_size / s->blocks;
2044
2045     imdct = s->imdct[transient ? 0 : s->duration];
2046
2047     if (coded_channels == 1) {
2048         for (i = 0; i < CELT_MAX_BANDS; i++)
2049             s->frame[0].energy[i] = FFMAX(s->frame[0].energy[i], s->frame[1].energy[i]);
2050     }
2051
2052     celt_decode_coarse_energy(s, rc);
2053     celt_decode_tf_changes   (s, rc, transient);
2054     celt_decode_allocation   (s, rc);
2055     celt_decode_fine_energy  (s, rc);
2056     celt_decode_bands        (s, rc);
2057
2058     if (s->anticollapse_bit)
2059         anticollapse = opus_getrawbits(rc, 1);
2060
2061     celt_decode_final_energy(s, rc, s->framebits - opus_rc_tell(rc));
2062
2063     /* apply anti-collapse processing and denormalization to
2064      * each coded channel */
2065     for (i = 0; i < s->coded_channels; i++) {
2066         CeltFrame *frame = &s->frame[i];
2067
2068         if (anticollapse)
2069             process_anticollapse(s, frame, s->coeffs[i]);
2070
2071         celt_denormalize(s, frame, s->coeffs[i]);
2072     }
2073
2074     /* stereo -> mono downmix */
2075     if (s->output_channels < s->coded_channels) {
2076         s->dsp->vector_fmac_scalar(s->coeffs[0], s->coeffs[1], 1.0, FFALIGN(frame_size, 16));
2077         imdct_scale = 0.5;
2078     } else if (s->output_channels > s->coded_channels)
2079         memcpy(s->coeffs[1], s->coeffs[0], frame_size * sizeof(float));
2080
2081     if (silence) {
2082         for (i = 0; i < 2; i++) {
2083             CeltFrame *frame = &s->frame[i];
2084
2085             for (j = 0; j < FF_ARRAY_ELEMS(frame->energy); j++)
2086                 frame->energy[j] = CELT_ENERGY_SILENCE;
2087         }
2088         memset(s->coeffs, 0, sizeof(s->coeffs));
2089     }
2090
2091     /* transform and output for each output channel */
2092     for (i = 0; i < s->output_channels; i++) {
2093         CeltFrame *frame = &s->frame[i];
2094         float m = frame->deemph_coeff;
2095
2096         /* iMDCT and overlap-add */
2097         for (j = 0; j < s->blocks; j++) {
2098             float *dst  = frame->buf + 1024 + j * s->blocksize;
2099
2100             imdct->imdct_half(imdct, dst + CELT_OVERLAP / 2, s->coeffs[i] + j,
2101                               s->blocks, imdct_scale);
2102             s->dsp->vector_fmul_window(dst, dst, dst + CELT_OVERLAP / 2,
2103                                       celt_window, CELT_OVERLAP / 2);
2104         }
2105
2106         /* postfilter */
2107         celt_postfilter(s, frame);
2108
2109         /* deemphasis and output scaling */
2110         for (j = 0; j < frame_size; j++) {
2111             float tmp = frame->buf[1024 - frame_size + j] + m;
2112             m = tmp * CELT_DEEMPH_COEFF;
2113             output[i][j] = tmp / 32768.;
2114         }
2115         frame->deemph_coeff = m;
2116     }
2117
2118     if (coded_channels == 1)
2119         memcpy(s->frame[1].energy, s->frame[0].energy, sizeof(s->frame[0].energy));
2120
2121     for (i = 0; i < 2; i++ ) {
2122         CeltFrame *frame = &s->frame[i];
2123
2124         if (!transient) {
2125             memcpy(frame->prev_energy[1], frame->prev_energy[0], sizeof(frame->prev_energy[0]));
2126             memcpy(frame->prev_energy[0], frame->energy,         sizeof(frame->prev_energy[0]));
2127         } else {
2128             for (j = 0; j < CELT_MAX_BANDS; j++)
2129                 frame->prev_energy[0][j] = FFMIN(frame->prev_energy[0][j], frame->energy[j]);
2130         }
2131
2132         for (j = 0; j < s->startband; j++) {
2133             frame->prev_energy[0][j] = CELT_ENERGY_SILENCE;
2134             frame->energy[j]         = 0.0;
2135         }
2136         for (j = s->endband; j < CELT_MAX_BANDS; j++) {
2137             frame->prev_energy[0][j] = CELT_ENERGY_SILENCE;
2138             frame->energy[j]         = 0.0;
2139         }
2140     }
2141
2142     s->seed = rc->range;
2143
2144     return 0;
2145 }
2146
2147 void ff_celt_flush(CeltContext *s)
2148 {
2149     int i, j;
2150
2151     if (s->flushed)
2152         return;
2153
2154     for (i = 0; i < 2; i++) {
2155         CeltFrame *frame = &s->frame[i];
2156
2157         for (j = 0; j < CELT_MAX_BANDS; j++)
2158             frame->prev_energy[0][j] = frame->prev_energy[1][j] = CELT_ENERGY_SILENCE;
2159
2160         memset(frame->energy, 0, sizeof(frame->energy));
2161         memset(frame->buf,    0, sizeof(frame->buf));
2162
2163         memset(frame->pf_gains,     0, sizeof(frame->pf_gains));
2164         memset(frame->pf_gains_old, 0, sizeof(frame->pf_gains_old));
2165         memset(frame->pf_gains_new, 0, sizeof(frame->pf_gains_new));
2166
2167         frame->deemph_coeff = 0.0;
2168     }
2169     s->seed = 0;
2170
2171     s->flushed = 1;
2172 }
2173
2174 void ff_celt_free(CeltContext **ps)
2175 {
2176     CeltContext *s = *ps;
2177     int i;
2178
2179     if (!s)
2180         return;
2181
2182     for (i = 0; i < FF_ARRAY_ELEMS(s->imdct); i++)
2183         ff_imdct15_uninit(&s->imdct[i]);
2184
2185     av_freep(&s->dsp);
2186     av_freep(ps);
2187 }
2188
2189 int ff_celt_init(AVCodecContext *avctx, CeltContext **ps, int output_channels)
2190 {
2191     CeltContext *s;
2192     int i, ret;
2193
2194     if (output_channels != 1 && output_channels != 2) {
2195         av_log(avctx, AV_LOG_ERROR, "Invalid number of output channels: %d\n",
2196                output_channels);
2197         return AVERROR(EINVAL);
2198     }
2199
2200     s = av_mallocz(sizeof(*s));
2201     if (!s)
2202         return AVERROR(ENOMEM);
2203
2204     s->avctx           = avctx;
2205     s->output_channels = output_channels;
2206
2207     for (i = 0; i < FF_ARRAY_ELEMS(s->imdct); i++) {
2208         ret = ff_imdct15_init(&s->imdct[i], i + 3);
2209         if (ret < 0)
2210             goto fail;
2211     }
2212
2213     s->dsp = avpriv_float_dsp_alloc(avctx->flags & AV_CODEC_FLAG_BITEXACT);
2214     if (!s->dsp) {
2215         ret = AVERROR(ENOMEM);
2216         goto fail;
2217     }
2218
2219     ff_celt_flush(s);
2220
2221     *ps = s;
2222
2223     return 0;
2224 fail:
2225     ff_celt_free(&s);
2226     return ret;
2227 }