libcelt/bands.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 <math.h>
  37 #include "bands.h"
  38 #include "modes.h"
  39 #include "vq.h"
  40 #include "cwrs.h"
  41
  42
  43 void exp_rotation(celt_norm_t *X, int len, float theta, int dir, int stride, int iter)
  44 {
  45    int i, k;
  46    celt_word16_t c, s;
  47    c = Q15ONE*cos(theta);
  48    s = dir*Q15ONE*sin(theta);
  49    for (k=0;k<iter;k++)
  50    {
  51       /* We could use MULT16_16_P15 instead of MULT16_16_Q15 for more accuracy,
  52          but at this point, I really don't think it's necessary */
  53       for (i=0;i<len-stride;i++)
  54       {
  55          celt_norm_t x1, x2;
  56          x1 = X[i];
  57          x2 = X[i+stride];
  58          X[i] = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
  59          X[i+stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
  60       }
  61       for (i=len-2*stride-1;i>=0;i--)
  62       {
  63          celt_norm_t x1, x2;
  64          x1 = X[i];
  65          x2 = X[i+stride];
  66          X[i] = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
  67          X[i+stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
  68       }
  69    }
  70 }
  71
  72 /* Compute the amplitude (sqrt energy) in each of the bands */
  73 void compute_band_energies(const CELTMode *m, celt_sig_t *X, celt_ener_t *bank)
  74 {
  75    int i, c, B, C;
  76    const int *eBands = m->eBands;
  77    B = m->nbMdctBlocks;
  78    C = m->nbChannels;
  79    for (c=0;c<C;c++)
  80    {
  81       for (i=0;i<m->nbEBands;i++)
  82       {
  83          int j;
  84          float sum = 1e-10;
  85          for (j=B*eBands[i];j<B*eBands[i+1];j++)
  86             sum += SIG_SCALING_1*SIG_SCALING_1*X[j*C+c]*X[j*C+c];
  87          bank[i*C+c] = ENER_SCALING*sqrt(sum);
  88          /*printf ("%f ", bank[i*C+c]);*/
  89       }
  90    }
  91    /*printf ("\n");*/
  92 }
  93
  94 /* Normalise each band such that the energy is one. */
  95 void normalise_bands(const CELTMode *m, celt_sig_t *freq, celt_norm_t *X, celt_ener_t *bank)
  96 {
  97    int i, c, B, C;
  98    const int *eBands = m->eBands;
  99    B = m->nbMdctBlocks;
 100    C = m->nbChannels;
 101    for (c=0;c<C;c++)
 102    {
 103       for (i=0;i<m->nbEBands;i++)
 104       {
 105          int j;
 106          float g = 1.f/(1e-10+ENER_SCALING_1*bank[i*C+c]*sqrt(C));
 107          for (j=B*eBands[i];j<B*eBands[i+1];j++)
 108             X[j*C+c] = NORM_SCALING*SIG_SCALING_1*freq[j*C+c]*g;
 109       }
 110    }
 111    for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
 112       X[i] = 0;
 113 }
 114
 115 #ifdef FIXED_POINT
 116 void renormalise_bands(const CELTMode *m, celt_norm_t *X)
 117 {
 118    int i;
 119    VARDECL(celt_ener_t *tmpE);
 120    VARDECL(celt_sig_t *freq);
 121    ALLOC(tmpE, m->nbEBands*m->nbChannels, celt_ener_t);
 122    ALLOC(freq, m->nbMdctBlocks*m->nbChannels*m->eBands[m->nbEBands+1], celt_sig_t);
 123    for (i=0;i<m->nbMdctBlocks*m->nbChannels*m->eBands[m->nbEBands+1];i++)
 124       freq[i] = SHL32(EXTEND32(X[i]), 10);
 125    compute_band_energies(m, freq, tmpE);
 126    normalise_bands(m, freq, X, tmpE);
 127 }
 128 #else
 129 void renormalise_bands(const CELTMode *m, celt_norm_t *X)
 130 {
 131    VARDECL(celt_ener_t *tmpE);
 132    ALLOC(tmpE, m->nbEBands*m->nbChannels, celt_ener_t);
 133    compute_band_energies(m, X, tmpE);
 134    normalise_bands(m, X, X, tmpE);
 135 }
 136 #endif
 137
 138 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
 139 void denormalise_bands(const CELTMode *m, celt_norm_t *X, celt_sig_t *freq, celt_ener_t *bank)
 140 {
 141    int i, c, B, C;
 142    const int *eBands = m->eBands;
 143    B = m->nbMdctBlocks;
 144    C = m->nbChannels;
 145    for (c=0;c<C;c++)
 146    {
 147       for (i=0;i<m->nbEBands;i++)
 148       {
 149          int j;
 150          float g = ENER_SCALING_1*sqrt(C)*bank[i*C+c];
 151          for (j=B*eBands[i];j<B*eBands[i+1];j++)
 152             freq[j*C+c] = NORM_SCALING_1*SIG_SCALING*X[j*C+c] * g;
 153       }
 154    }
 155    for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
 156       freq[i] = 0;
 157 }
 158
 159
 160 /* Compute the best gain for each "pitch band" */
 161 void compute_pitch_gain(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, celt_pgain_t *gains, celt_ener_t *bank)
 162 {
 163    int i, B;
 164    const int *eBands = m->eBands;
 165    const int *pBands = m->pBands;
 166    VARDECL(float *w);
 167    B = m->nbMdctBlocks*m->nbChannels;
 168    ALLOC(w, B*eBands[m->nbEBands], float);
 169    for (i=0;i<m->nbEBands;i++)
 170    {
 171       int j;
 172       for (j=B*eBands[i];j<B*eBands[i+1];j++)
 173          w[j] = bank[i]*ENER_SCALING_1;
 174    }
 175
 176
 177    for (i=0;i<m->nbPBands;i++)
 178    {
 179       float Sxy=0;
 180       float Sxx = 0;
 181       int j;
 182       float gain;
 183       for (j=B*pBands[i];j<B*pBands[i+1];j++)
 184       {
 185          Sxy += 1.f*X[j]*P[j]*w[j];
 186          Sxx += 1.f*X[j]*X[j]*w[j];
 187       }
 188       gain = Sxy/(1e-10*NORM_SCALING*NORM_SCALING+Sxx);
 189       if (gain > 1.f)
 190          gain = 1.f;
 191       if (gain < 0.0f)
 192          gain = 0.0f;
 193       /* We need to be a bit conservative, otherwise residual doesn't quantise well */
 194       gain *= .9f;
 195       gains[i] = PGAIN_SCALING*gain;
 196       /*printf ("%f ", 1-sqrt(1-gain*gain));*/
 197    }
 198    /*if(rand()%10==0)
 199    {
 200       for (i=0;i<m->nbPBands;i++)
 201          printf ("%f ", 1-sqrt(1-gains[i]*gains[i]));
 202       printf ("\n");
 203    }*/
 204    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
 205       P[i] = 0;
 206 }
 207
 208 /* Apply the (quantised) gain to each "pitch band" */
 209 void pitch_quant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, celt_pgain_t *gains)
 210 {
 211    int i, B;
 212    const int *pBands = m->pBands;
 213    B = m->nbMdctBlocks*m->nbChannels;
 214    for (i=0;i<m->nbPBands;i++)
 215    {
 216       int j;
 217       for (j=B*pBands[i];j<B*pBands[i+1];j++)
 218          P[j] *= PGAIN_SCALING_1*gains[i];
 219       /*printf ("%f ", gain);*/
 220    }
 221    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
 222       P[i] = 0;
 223 }
 224
 225
 226 /* Quantisation of the residual */
 227 void quant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, celt_mask_t *W, int total_bits, ec_enc *enc)
 228 {
 229    int i, j, B, bits;
 230    const int *eBands = m->eBands;
 231    float alpha = .7;
 232    VARDECL(celt_norm_t *norm);
 233    VARDECL(int *pulses);
 234    VARDECL(int *offsets);
 235
 236    B = m->nbMdctBlocks*m->nbChannels;
 237
 238    ALLOC(norm, B*eBands[m->nbEBands+1], celt_norm_t);
 239    ALLOC(pulses, m->nbEBands, int);
 240    ALLOC(offsets, m->nbEBands, int);
 241
 242    for (i=0;i<m->nbEBands;i++)
 243       offsets[i] = 0;
 244    /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
 245    bits = total_bits - ec_enc_tell(enc, 0) - 1;
 246    compute_allocation(m, offsets, bits, pulses);
 247
 248    /*printf("bits left: %d\n", bits);
 249    for (i=0;i<m->nbEBands;i++)
 250       printf ("%d ", pulses[i]);
 251    printf ("\n");*/
 252    /*printf ("%d %d\n", ec_enc_tell(enc, 0), compute_allocation(m, m->nbPulses));*/
 253    for (i=0;i<m->nbEBands;i++)
 254    {
 255       int q;
 256       float theta, n;
 257       q = pulses[i];
 258       /*Scale factor of .0625f is just there to prevent overflows in fixed-point
 259        (has no effect on float)*/
 260       n = .0625f*sqrt(B*(eBands[i+1]-eBands[i]));
 261       theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+q);
 262
 263       /* If pitch isn't available, use intra-frame prediction */
 264       if (eBands[i] >= m->pitchEnd || q<=0)
 265       {
 266          q -= 1;
 267          alpha = 0;
 268          if (q<0)
 269             intra_fold(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), norm, P+B*eBands[i], B, eBands[i], eBands[m->nbEBands+1]);
 270          else
 271             intra_prediction(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], enc);
 272       } else {
 273          alpha = .7;
 274       }
 275
 276       if (q > 0)
 277       {
 278          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
 279          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
 280          alg_quant(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], alpha, enc);
 281          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
 282       }
 283       for (j=B*eBands[i];j<B*eBands[i+1];j++)
 284          norm[j] = X[j] * n;
 285    }
 286    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
 287       X[i] = 0;
 288 }
 289
 290 /* Decoding of the residual */
 291 void unquant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, int total_bits, ec_dec *dec)
 292 {
 293    int i, j, B, bits;
 294    const int *eBands = m->eBands;
 295    float alpha = .7;
 296    VARDECL(celt_norm_t *norm);
 297    VARDECL(int *pulses);
 298    VARDECL(int *offsets);
 299
 300    B = m->nbMdctBlocks*m->nbChannels;
 301
 302    ALLOC(norm, B*eBands[m->nbEBands+1], celt_norm_t);
 303    ALLOC(pulses, m->nbEBands, int);
 304    ALLOC(offsets, m->nbEBands, int);
 305
 306    for (i=0;i<m->nbEBands;i++)
 307       offsets[i] = 0;
 308    /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
 309    bits = total_bits - ec_dec_tell(dec, 0) - 1;
 310    compute_allocation(m, offsets, bits, pulses);
 311
 312    for (i=0;i<m->nbEBands;i++)
 313    {
 314       int q;
 315       float theta, n;
 316       q = pulses[i];
 317       /*Scale factor of .0625f is just there to prevent overflows in fixed-point
 318       (has no effect on float)*/
 319       n = .0625f*sqrt(B*(eBands[i+1]-eBands[i]));
 320       theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+q);
 321
 322       /* If pitch isn't available, use intra-frame prediction */
 323       if (eBands[i] >= m->pitchEnd || q<=0)
 324       {
 325          q -= 1;
 326          alpha = 0;
 327          if (q<0)
 328             intra_fold(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), norm, P+B*eBands[i], B, eBands[i], eBands[m->nbEBands+1]);
 329          else
 330             intra_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], dec);
 331       } else {
 332          alpha = .7;
 333       }
 334
 335       if (q > 0)
 336       {
 337          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
 338          alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], alpha, dec);
 339          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
 340       }
 341       for (j=B*eBands[i];j<B*eBands[i+1];j++)
 342          norm[j] = X[j] * n;
 343    }
 344    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
 345       X[i] = 0;
 346 }
 347
 348 void stereo_mix(const CELTMode *m, celt_norm_t *X, celt_ener_t *bank, int dir)
 349 {
 350    int i, B, C;
 351    const int *eBands = m->eBands;
 352    B = m->nbMdctBlocks;
 353    C = m->nbChannels;
 354    for (i=0;i<m->nbEBands;i++)
 355    {
 356       int j;
 357       float left, right;
 358       float a1, a2;
 359       left = bank[i*C];
 360       right = bank[i*C+1];
 361       a1 = left/sqrt(.01+left*left+right*right);
 362       a2 = dir*right/sqrt(.01+left*left+right*right);
 363       for (j=B*eBands[i];j<B*eBands[i+1];j++)
 364       {
 365          float r, l;
 366          l = X[j*C];
 367          r = X[j*C+1];
 368          X[j*C] = a1*l + a2*r;
 369          X[j*C+1] = a1*r - a2*l;
 370       }
 371    }
 372    for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
 373       X[i] = 0;
 374
 375 }