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