libcelt/celt.c

   1 /* (C) 2007 Jean-Marc Valin, CSIRO
   2 */
   3 /*
   4    Redistribution and use in source and binary forms, with or without
   5    modification, are permitted provided that the following conditions
   6    are met:
   7
   8    - Redistributions of source code must retain the above copyright
   9    notice, this list of conditions and the following disclaimer.
  10
  11    - Redistributions in binary form must reproduce the above copyright
  12    notice, this list of conditions and the following disclaimer in the
  13    documentation and/or other materials provided with the distribution.
  14
  15    - Neither the name of the Xiph.org Foundation nor the names of its
  16    contributors may be used to endorse or promote products derived from
  17    this software without specific prior written permission.
  18
  19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
  23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  30 */
  31
  32 #ifdef HAVE_CONFIG_H
  33 #include "config.h"
  34 #endif
  35
  36 #include "os_support.h"
  37 #include "mdct.h"
  38 #include <math.h>
  39 #include "celt.h"
  40 #include "pitch.h"
  41 #include "kiss_fftr.h"
  42 #include "bands.h"
  43 #include "modes.h"
  44 #include "entcode.h"
  45 #include "quant_pitch.h"
  46 #include "quant_bands.h"
  47 #include "psy.h"
  48 #include "rate.h"
  49
  50 #define MAX_PERIOD 1024
  51
  52 #ifndef M_PI
  53 #define M_PI 3.14159263
  54 #endif
  55
  56 /** Encoder state
  57  @brief Encoder state
  58  */
  59 struct CELTEncoder {
  60    const CELTMode *mode;     /**< Mode used by the encoder */
  61    int frame_size;
  62    int block_size;
  63    int nb_blocks;
  64    int overlap;
  65    int channels;
  66
  67    ec_byte_buffer buf;
  68    ec_enc         enc;
  69
  70    float preemph;
  71    float *preemph_memE;
  72    float *preemph_memD;
  73
  74    mdct_lookup mdct_lookup;
  75    kiss_fftr_cfg fft;
  76    struct PsyDecay psy;
  77
  78    float *window;
  79    float *in_mem;
  80    float *mdct_overlap;
  81    float *out_mem;
  82
  83    float *oldBandE;
  84 };
  85
  86
  87
  88 CELTEncoder *celt_encoder_create(const CELTMode *mode)
  89 {
  90    int i, N, B, C, N4;
  91    CELTEncoder *st;
  92
  93    if (check_mode(mode) != CELT_OK)
  94       return NULL;
  95
  96    N = mode->mdctSize;
  97    B = mode->nbMdctBlocks;
  98    C = mode->nbChannels;
  99    st = celt_alloc(sizeof(CELTEncoder));
 100
 101    st->mode = mode;
 102    st->frame_size = B*N;
 103    st->block_size = N;
 104    st->nb_blocks  = B;
 105    st->overlap = mode->overlap;
 106
 107    N4 = (N-st->overlap)/2;
 108    ec_byte_writeinit(&st->buf);
 109    ec_enc_init(&st->enc,&st->buf);
 110
 111    mdct_init(&st->mdct_lookup, 2*N);
 112    st->fft = kiss_fftr_alloc(MAX_PERIOD*C, 0, 0);
 113    psydecay_init(&st->psy, MAX_PERIOD*C/2, st->mode->Fs);
 114
 115    st->window = celt_alloc(2*N*sizeof(float));
 116    st->in_mem = celt_alloc(N*C*sizeof(float));
 117    st->mdct_overlap = celt_alloc(N*C*sizeof(float));
 118    st->out_mem = celt_alloc(MAX_PERIOD*C*sizeof(float));
 119    for (i=0;i<2*N;i++)
 120       st->window[i] = 0;
 121    for (i=0;i<st->overlap;i++)
 122       st->window[N4+i] = st->window[2*N-N4-i-1]
 123             = sin(.5*M_PI* sin(.5*M_PI*(i+.5)/st->overlap) * sin(.5*M_PI*(i+.5)/st->overlap));
 124    for (i=0;i<2*N4;i++)
 125       st->window[N-N4+i] = 1;
 126    st->oldBandE = celt_alloc(C*mode->nbEBands*sizeof(float));
 127
 128    st->preemph = 0.8;
 129    st->preemph_memE = celt_alloc(C*sizeof(float));;
 130    st->preemph_memD = celt_alloc(C*sizeof(float));;
 131
 132    return st;
 133 }
 134
 135 void celt_encoder_destroy(CELTEncoder *st)
 136 {
 137    if (st == NULL)
 138    {
 139       celt_warning("NULL passed to celt_encoder_destroy");
 140       return;
 141    }
 142    if (check_mode(st->mode) != CELT_OK)
 143       return;
 144
 145    ec_byte_writeclear(&st->buf);
 146
 147    mdct_clear(&st->mdct_lookup);
 148    kiss_fft_free(st->fft);
 149    psydecay_clear(&st->psy);
 150
 151    celt_free(st->window);
 152    celt_free(st->in_mem);
 153    celt_free(st->mdct_overlap);
 154    celt_free(st->out_mem);
 155
 156    celt_free(st->oldBandE);
 157
 158    celt_free(st->preemph_memE);
 159    celt_free(st->preemph_memD);
 160
 161    celt_free(st);
 162 }
 163
 164 /** Apply window and compute the MDCT for all sub-frames and all channels in a frame */
 165 static float compute_mdcts(mdct_lookup *mdct_lookup, float *window, float *in, float *out, int N, int B, int C)
 166 {
 167    int i, c;
 168    float E = 1e-15;
 169    VARDECL(float *x);
 170    VARDECL(float *tmp);
 171    ALLOC(x, 2*N, float);
 172    ALLOC(tmp, N, float);
 173    for (c=0;c<C;c++)
 174    {
 175       for (i=0;i<B;i++)
 176       {
 177          int j;
 178          for (j=0;j<2*N;j++)
 179          {
 180             x[j] = window[j]*in[C*i*N+C*j+c];
 181             E += x[j]*x[j];
 182          }
 183          mdct_forward(mdct_lookup, x, tmp);
 184          /* Interleaving the sub-frames */
 185          for (j=0;j<N;j++)
 186             out[C*B*j+C*i+c] = tmp[j];
 187       }
 188    }
 189    return E;
 190 }
 191
 192 /** Compute the IMDCT and apply window for all sub-frames and all channels in a frame */
 193 static void compute_inv_mdcts(mdct_lookup *mdct_lookup, float *window, float *X, float *out_mem, float *mdct_overlap, int N, int overlap, int B, int C)
 194 {
 195    int i, c, N4;
 196    VARDECL(float *x);
 197    VARDECL(float *tmp);
 198    ALLOC(x, 2*N, float);
 199    ALLOC(tmp, N, float);
 200    N4 = (N-overlap)/2;
 201    for (c=0;c<C;c++)
 202    {
 203       for (i=0;i<B;i++)
 204       {
 205          int j;
 206          /* De-interleaving the sub-frames */
 207          for (j=0;j<N;j++)
 208             tmp[j] = X[C*B*j+C*i+c];
 209          mdct_backward(mdct_lookup, tmp, x);
 210          for (j=0;j<2*N;j++)
 211             x[j] = window[j]*x[j];
 212          for (j=0;j<overlap;j++)
 213             out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] = x[N4+j]+mdct_overlap[C*j+c];
 214          for (j=0;j<2*N4;j++)
 215             out_mem[C*(MAX_PERIOD+(i-B)*N)+C*(j+overlap)+c] = x[j+N4+overlap];
 216          for (j=0;j<overlap;j++)
 217             mdct_overlap[C*j+c] = x[N+N4+j];
 218       }
 219    }
 220 }
 221
 222 int celt_encode(CELTEncoder *st, celt_int16_t *pcm, unsigned char *compressed, int nbCompressedBytes)
 223 {
 224    int i, c, N, B, C, N4;
 225    int has_pitch;
 226    int pitch_index;
 227    float curr_power, pitch_power;
 228    VARDECL(float *in);
 229    VARDECL(float *X);
 230    VARDECL(float *P);
 231    VARDECL(float *mask);
 232    VARDECL(float *bandE);
 233    VARDECL(float *gains);
 234
 235    if (check_mode(st->mode) != CELT_OK)
 236       return CELT_INVALID_MODE;
 237
 238    N = st->block_size;
 239    B = st->nb_blocks;
 240    C = st->mode->nbChannels;
 241    ALLOC(in, (B+1)*C*N, float);
 242    ALLOC(X, B*C*N, float);         /**< Interleaved signal MDCTs */
 243    ALLOC(P, B*C*N, float);         /**< Interleaved pitch MDCTs*/
 244    ALLOC(mask, B*C*N, float);      /**< Masking curve */
 245    ALLOC(bandE,st->mode->nbEBands*C, float);
 246    ALLOC(gains,st->mode->nbPBands, float);
 247
 248    N4 = (N-st->overlap)/2;
 249
 250    for (c=0;c<C;c++)
 251    {
 252       for (i=0;i<N4;i++)
 253          in[C*i+c] = 0;
 254       for (i=0;i<st->overlap;i++)
 255          in[C*(i+N4)+c] = st->in_mem[C*i+c];
 256       for (i=0;i<B*N;i++)
 257       {
 258          float tmp = pcm[C*i+c];
 259          in[C*(i+st->overlap+N4)+c] = tmp - st->preemph*st->preemph_memE[c];
 260          st->preemph_memE[c] = tmp;
 261       }
 262       for (i=N*(B+1)-N4;i<N*(B+1);i++)
 263          in[C*i+c] = 0;
 264       for (i=0;i<st->overlap;i++)
 265          st->in_mem[C*i+c] = in[C*(N*(B+1)-N4-st->overlap+i)+c];
 266    }
 267    /*for (i=0;i<(B+1)*C*N;i++) printf ("%f(%d) ", in[i], i); printf ("\n");*/
 268    /* Compute MDCTs */
 269    curr_power = compute_mdcts(&st->mdct_lookup, st->window, in, X, N, B, C);
 270
 271 #if 0 /* Mask disabled until it can be made to do something useful */
 272    compute_mdct_masking(X, mask, B*C*N, st->Fs);
 273
 274    /* Invert and stretch the mask to length of X
 275       For some reason, I get better results by using the sqrt instead,
 276       although there's no valid reason to. Must investigate further */
 277    for (i=0;i<B*C*N;i++)
 278       mask[i] = 1/(.1+mask[i]);
 279 #else
 280    for (i=0;i<B*C*N;i++)
 281       mask[i] = 1;
 282 #endif
 283    /* Pitch analysis */
 284    for (c=0;c<C;c++)
 285    {
 286       for (i=0;i<N;i++)
 287       {
 288          in[C*i+c] *= st->window[i];
 289          in[C*(B*N+i)+c] *= st->window[N+i];
 290       }
 291    }
 292    find_spectral_pitch(st->fft, &st->psy, in, st->out_mem, MAX_PERIOD, (B+1)*N, C, &pitch_index);
 293
 294    /* Compute MDCTs of the pitch part */
 295    pitch_power = compute_mdcts(&st->mdct_lookup, st->window, st->out_mem+pitch_index*C, P, N, B, C);
 296
 297    /*printf ("%f %f\n", curr_power, pitch_power);*/
 298    /*int j;
 299    for (j=0;j<B*N;j++)
 300       printf ("%f ", X[j]);
 301    for (j=0;j<B*N;j++)
 302       printf ("%f ", P[j]);
 303    printf ("\n");*/
 304
 305    /* Band normalisation */
 306    compute_band_energies(st->mode, X, bandE);
 307    normalise_bands(st->mode, X, bandE);
 308    /*for (i=0;i<st->mode->nbEBands;i++)printf("%f ", bandE[i]);printf("\n");*/
 309    /*for (i=0;i<N*B*C;i++)printf("%f ", X[i]);printf("\n");*/
 310
 311    quant_energy(st->mode, bandE, st->oldBandE, nbCompressedBytes*8/3, &st->enc);
 312
 313    if (C==2)
 314    {
 315       stereo_mix(st->mode, X, bandE, 1);
 316    }
 317
 318    /* Check if we can safely use the pitch (i.e. effective gain isn't too high) */
 319    if (curr_power + 1e5f < 10.f*pitch_power)
 320    {
 321       /* Normalise the pitch vector as well (discard the energies) */
 322       VARDECL(float *bandEp);
 323       ALLOC(bandEp, st->mode->nbEBands*st->mode->nbChannels, float);
 324       compute_band_energies(st->mode, P, bandEp);
 325       normalise_bands(st->mode, P, bandEp);
 326
 327       if (C==2)
 328          stereo_mix(st->mode, P, bandE, 1);
 329       /* Simulates intensity stereo */
 330       /*for (i=30;i<N*B;i++)
 331          X[i*C+1] = P[i*C+1] = 0;*/
 332
 333       /* Pitch prediction */
 334       compute_pitch_gain(st->mode, X, P, gains, bandE);
 335       has_pitch = quant_pitch(gains, st->mode->nbPBands, &st->enc);
 336       if (has_pitch)
 337          ec_enc_uint(&st->enc, pitch_index, MAX_PERIOD-(B+1)*N);
 338    } else {
 339       /* No pitch, so we just pretend we found a gain of zero */
 340       for (i=0;i<st->mode->nbPBands;i++)
 341          gains[i] = 0;
 342       ec_enc_uint(&st->enc, 0, 128);
 343       for (i=0;i<B*C*N;i++)
 344          P[i] = 0;
 345    }
 346
 347
 348    pitch_quant_bands(st->mode, X, P, gains);
 349
 350    /*for (i=0;i<B*N;i++) printf("%f ",P[i]);printf("\n");*/
 351    /* Compute residual that we're going to encode */
 352    for (i=0;i<B*C*N;i++)
 353       X[i] -= P[i];
 354
 355    /*float sum=0;
 356    for (i=0;i<B*N;i++)
 357       sum += X[i]*X[i];
 358    printf ("%f\n", sum);*/
 359    /* Residual quantisation */
 360    quant_bands(st->mode, X, P, mask, nbCompressedBytes*8, &st->enc);
 361
 362    if (C==2)
 363       stereo_mix(st->mode, X, bandE, -1);
 364
 365    renormalise_bands(st->mode, X);
 366    /* Synthesis */
 367    denormalise_bands(st->mode, X, bandE);
 368
 369
 370    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
 371
 372    compute_inv_mdcts(&st->mdct_lookup, st->window, X, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
 373    /* De-emphasis and put everything back at the right place in the synthesis history */
 374    for (c=0;c<C;c++)
 375    {
 376       for (i=0;i<B;i++)
 377       {
 378          int j;
 379          for (j=0;j<N;j++)
 380          {
 381             float tmp = st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] + st->preemph*st->preemph_memD[c];
 382             st->preemph_memD[c] = tmp;
 383             if (tmp > 32767) tmp = 32767;
 384             if (tmp < -32767) tmp = -32767;
 385             pcm[C*i*N+C*j+c] = (short)floor(.5+tmp);
 386          }
 387       }
 388    }
 389
 390    if (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8 - 7)
 391       celt_warning_int ("many unused bits: ", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));
 392    /*printf ("%d\n", ec_enc_tell(&st->enc, 0)-8*nbCompressedBytes);*/
 393    /* Finishing the stream with a 0101... pattern so that the decoder can check is everything's right */
 394    {
 395       int val = 0;
 396       while (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8)
 397       {
 398          ec_enc_uint(&st->enc, val, 2);
 399          val = 1-val;
 400       }
 401    }
 402    ec_enc_done(&st->enc);
 403    {
 404       unsigned char *data;
 405       int nbBytes = ec_byte_bytes(&st->buf);
 406       if (nbBytes > nbCompressedBytes)
 407       {
 408          celt_warning_int ("got too many bytes:", nbBytes);
 409          return CELT_INTERNAL_ERROR;
 410       }
 411       /*printf ("%d\n", *nbBytes);*/
 412       data = ec_byte_get_buffer(&st->buf);
 413       for (i=0;i<nbBytes;i++)
 414          compressed[i] = data[i];
 415       for (;i<nbCompressedBytes;i++)
 416          compressed[i] = 0;
 417    }
 418    /* Reset the packing for the next encoding */
 419    ec_byte_reset(&st->buf);
 420    ec_enc_init(&st->enc,&st->buf);
 421
 422    return nbCompressedBytes;
 423 }
 424
 425
 426 /****************************************************************************/
 427 /*                                                                          */
 428 /*                                DECODER                                   */
 429 /*                                                                          */
 430 /****************************************************************************/
 431
 432
 433 /** Decoder state
 434  @brief Decoder state
 435  */
 436 struct CELTDecoder {
 437    const CELTMode *mode;
 438    int frame_size;
 439    int block_size;
 440    int nb_blocks;
 441    int overlap;
 442
 443    ec_byte_buffer buf;
 444    ec_enc         enc;
 445
 446    float preemph;
 447    float *preemph_memD;
 448
 449    mdct_lookup mdct_lookup;
 450
 451    float *window;
 452    float *mdct_overlap;
 453    float *out_mem;
 454
 455    float *oldBandE;
 456
 457    int last_pitch_index;
 458 };
 459
 460 CELTDecoder *celt_decoder_create(const CELTMode *mode)
 461 {
 462    int i, N, B, C, N4;
 463    CELTDecoder *st;
 464
 465    if (check_mode(mode) != CELT_OK)
 466       return NULL;
 467
 468    N = mode->mdctSize;
 469    B = mode->nbMdctBlocks;
 470    C = mode->nbChannels;
 471    st = celt_alloc(sizeof(CELTDecoder));
 472
 473    st->mode = mode;
 474    st->frame_size = B*N;
 475    st->block_size = N;
 476    st->nb_blocks  = B;
 477    st->overlap = mode->overlap;
 478
 479    N4 = (N-st->overlap)/2;
 480
 481    mdct_init(&st->mdct_lookup, 2*N);
 482
 483    st->window = celt_alloc(2*N*sizeof(float));
 484    st->mdct_overlap = celt_alloc(N*C*sizeof(float));
 485    st->out_mem = celt_alloc(MAX_PERIOD*C*sizeof(float));
 486
 487    for (i=0;i<2*N;i++)
 488       st->window[i] = 0;
 489    for (i=0;i<st->overlap;i++)
 490       st->window[N4+i] = st->window[2*N-N4-i-1]
 491             = sin(.5*M_PI* sin(.5*M_PI*(i+.5)/st->overlap) * sin(.5*M_PI*(i+.5)/st->overlap));
 492    for (i=0;i<2*N4;i++)
 493       st->window[N-N4+i] = 1;
 494
 495    st->oldBandE = celt_alloc(C*mode->nbEBands*sizeof(float));
 496
 497    st->preemph = 0.8;
 498    st->preemph_memD = celt_alloc(C*sizeof(float));;
 499
 500    st->last_pitch_index = 0;
 501    return st;
 502 }
 503
 504 void celt_decoder_destroy(CELTDecoder *st)
 505 {
 506    if (st == NULL)
 507    {
 508       celt_warning("NULL passed to celt_encoder_destroy");
 509       return;
 510    }
 511    if (check_mode(st->mode) != CELT_OK)
 512       return;
 513
 514    mdct_clear(&st->mdct_lookup);
 515
 516    celt_free(st->window);
 517    celt_free(st->mdct_overlap);
 518    celt_free(st->out_mem);
 519
 520    celt_free(st->oldBandE);
 521
 522    celt_free(st->preemph_memD);
 523
 524    celt_free(st);
 525 }
 526
 527 /** Handles lost packets by just copying past data with the same offset as the last
 528     pitch period */
 529 static void celt_decode_lost(CELTDecoder *st, short *pcm)
 530 {
 531    int i, c, N, B, C;
 532    int pitch_index;
 533    VARDECL(float *X);
 534    N = st->block_size;
 535    B = st->nb_blocks;
 536    C = st->mode->nbChannels;
 537    ALLOC(X,C*B*N, float);         /**< Interleaved signal MDCTs */
 538
 539    pitch_index = st->last_pitch_index;
 540
 541    /* Use the pitch MDCT as the "guessed" signal */
 542    compute_mdcts(&st->mdct_lookup, st->window, st->out_mem+pitch_index*C, X, N, B, C);
 543
 544    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
 545    /* Compute inverse MDCTs */
 546    compute_inv_mdcts(&st->mdct_lookup, st->window, X, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
 547
 548    for (c=0;c<C;c++)
 549    {
 550       for (i=0;i<B;i++)
 551       {
 552          int j;
 553          for (j=0;j<N;j++)
 554          {
 555             float tmp = st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] + st->preemph*st->preemph_memD[c];
 556             st->preemph_memD[c] = tmp;
 557             if (tmp > 32767) tmp = 32767;
 558             if (tmp < -32767) tmp = -32767;
 559             pcm[C*i*N+C*j+c] = (short)floor(.5+tmp);
 560          }
 561       }
 562    }
 563 }
 564
 565 int celt_decode(CELTDecoder *st, unsigned char *data, int len, celt_int16_t *pcm)
 566 {
 567    int i, c, N, B, C;
 568    int has_pitch;
 569    int pitch_index;
 570    ec_dec dec;
 571    ec_byte_buffer buf;
 572    VARDECL(float *X);
 573    VARDECL(float *P);
 574    VARDECL(float *bandE);
 575    VARDECL(float *gains);
 576
 577    if (check_mode(st->mode) != CELT_OK)
 578       return CELT_INVALID_MODE;
 579
 580    N = st->block_size;
 581    B = st->nb_blocks;
 582    C = st->mode->nbChannels;
 583
 584    ALLOC(X, C*B*N, float);         /**< Interleaved signal MDCTs */
 585    ALLOC(P, C*B*N, float);         /**< Interleaved pitch MDCTs*/
 586    ALLOC(bandE, st->mode->nbEBands*C, float);
 587    ALLOC(gains, st->mode->nbPBands, float);
 588
 589    if (check_mode(st->mode) != CELT_OK)
 590       return CELT_INVALID_MODE;
 591    if (data == NULL)
 592    {
 593       celt_decode_lost(st, pcm);
 594       return 0;
 595    }
 596
 597    ec_byte_readinit(&buf,data,len);
 598    ec_dec_init(&dec,&buf);
 599
 600    /* Get band energies */
 601    unquant_energy(st->mode, bandE, st->oldBandE, len*8/3, &dec);
 602
 603    /* Get the pitch gains */
 604    has_pitch = unquant_pitch(gains, st->mode->nbPBands, &dec);
 605
 606    /* Get the pitch index */
 607    if (has_pitch)
 608    {
 609       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-(B+1)*N);
 610       st->last_pitch_index = pitch_index;
 611    } else {
 612       /* FIXME: We could be more intelligent here and just not compute the MDCT */
 613       pitch_index = 0;
 614    }
 615
 616    /* Pitch MDCT */
 617    compute_mdcts(&st->mdct_lookup, st->window, st->out_mem+pitch_index*C, P, N, B, C);
 618
 619    {
 620       VARDECL(float *bandEp);
 621       ALLOC(bandEp, st->mode->nbEBands*C, float);
 622       compute_band_energies(st->mode, P, bandEp);
 623       normalise_bands(st->mode, P, bandEp);
 624    }
 625
 626    if (C==2)
 627       stereo_mix(st->mode, P, bandE, 1);
 628
 629    /* Apply pitch gains */
 630    pitch_quant_bands(st->mode, X, P, gains);
 631
 632    /* Decode fixed codebook and merge with pitch */
 633    unquant_bands(st->mode, X, P, len*8, &dec);
 634
 635    if (C==2)
 636       stereo_mix(st->mode, X, bandE, -1);
 637
 638    renormalise_bands(st->mode, X);
 639
 640    /* Synthesis */
 641    denormalise_bands(st->mode, X, bandE);
 642
 643
 644    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
 645    /* Compute inverse MDCTs */
 646    compute_inv_mdcts(&st->mdct_lookup, st->window, X, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
 647
 648    for (c=0;c<C;c++)
 649    {
 650       for (i=0;i<B;i++)
 651       {
 652          int j;
 653          for (j=0;j<N;j++)
 654          {
 655             float tmp = st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] + st->preemph*st->preemph_memD[c];
 656             st->preemph_memD[c] = tmp;
 657             if (tmp > 32767) tmp = 32767;
 658             if (tmp < -32767) tmp = -32767;
 659             pcm[C*i*N+C*j+c] = (short)floor(.5+tmp);
 660          }
 661       }
 662    }
 663
 664    {
 665       int val = 0;
 666       while (ec_dec_tell(&dec, 0) < len*8)
 667       {
 668          if (ec_dec_uint(&dec, 2) != val)
 669          {
 670             celt_warning("decode error");
 671             return CELT_CORRUPTED_DATA;
 672          }
 673          val = 1-val;
 674       }
 675    }
 676
 677    return 0;
 678    /*printf ("\n");*/
 679 }
 680