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, B;
  95    const int *eBands = m->eBands;
  96    B = m->nbMdctBlocks*m->nbChannels;
  97    for (i=0;i<m->nbEBands;i++)
  98    {
  99       int j;
 100       bank[i] = 1e-10;
 101       for (j=B*eBands[i];j<B*eBands[i+1];j++)
 102          bank[i] += X[j]*X[j];
 103       bank[i] = sqrt(bank[i]);
 104    }
 105 }
 106
 107 /* Normalise each band such that the energy is one. */
 108 void normalise_bands(const CELTMode *m, float *X, float *bank)
 109 {
 110    int i, B;
 111    const int *eBands = m->eBands;
 112    B = m->nbMdctBlocks*m->nbChannels;
 113    for (i=0;i<m->nbEBands;i++)
 114    {
 115       int j;
 116       float x = 1.f/(1e-10+bank[i]);
 117       for (j=B*eBands[i];j<B*eBands[i+1];j++)
 118          X[j] *= x;
 119    }
 120    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
 121       X[i] = 0;
 122 }
 123
 124 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
 125 void denormalise_bands(const CELTMode *m, float *X, float *bank)
 126 {
 127    int i, B;
 128    const int *eBands = m->eBands;
 129    B = m->nbMdctBlocks*m->nbChannels;
 130    for (i=0;i<m->nbEBands;i++)
 131    {
 132       int j;
 133       float x = bank[i];
 134       for (j=B*eBands[i];j<B*eBands[i+1];j++)
 135          X[j] *= x;
 136    }
 137    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
 138       X[i] = 0;
 139 }
 140
 141
 142 /* Compute the best gain for each "pitch band" */
 143 void compute_pitch_gain(const CELTMode *m, float *X, float *P, float *gains, float *bank)
 144 {
 145    int i, B;
 146    const int *eBands = m->eBands;
 147    const int *pBands = m->pBands;
 148    B = m->nbMdctBlocks*m->nbChannels;
 149    float w[B*eBands[m->nbEBands]];
 150    for (i=0;i<m->nbEBands;i++)
 151    {
 152       int j;
 153       for (j=B*eBands[i];j<B*eBands[i+1];j++)
 154          w[j] = bank[i];
 155    }
 156
 157
 158    for (i=0;i<m->nbPBands;i++)
 159    {
 160       float Sxy=0;
 161       float Sxx = 0;
 162       int j;
 163       float gain;
 164       for (j=B*pBands[i];j<B*pBands[i+1];j++)
 165       {
 166          Sxy += X[j]*P[j]*w[j];
 167          Sxx += X[j]*X[j]*w[j];
 168       }
 169       gain = Sxy/(1e-10+Sxx);
 170       //gain = Sxy/(2*(pbank[i+1]-pbank[i]));
 171       //if (i<3)
 172       //gain *= 1+.02*gain;
 173       if (gain > .90)
 174          gain = .90;
 175       if (gain < 0.0)
 176          gain = 0.0;
 177
 178       gains[i] = gain;
 179    }
 180    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
 181       P[i] = 0;
 182 }
 183
 184 /* Apply the (quantised) gain to each "pitch band" */
 185 void pitch_quant_bands(const CELTMode *m, float *X, float *P, float *gains)
 186 {
 187    int i, B;
 188    const int *pBands = m->pBands;
 189    B = m->nbMdctBlocks*m->nbChannels;
 190    for (i=0;i<m->nbPBands;i++)
 191    {
 192       int j;
 193       for (j=B*pBands[i];j<B*pBands[i+1];j++)
 194          P[j] *= gains[i];
 195       //printf ("%f ", gain);
 196    }
 197    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
 198       P[i] = 0;
 199 }
 200
 201 /* Quantisation of the residual */
 202 void quant_bands(const CELTMode *m, float *X, float *P, float *W, ec_enc *enc)
 203 {
 204    int i, j, B;
 205    const int *eBands = m->eBands;
 206    B = m->nbMdctBlocks*m->nbChannels;
 207    float norm[B*eBands[m->nbEBands+1]];
 208
 209    for (i=0;i<m->nbEBands;i++)
 210    {
 211       int q;
 212       float theta, n;
 213       q = m->nbPulses[i];
 214       n = sqrt(B*(eBands[i+1]-eBands[i]));
 215       theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+abs(m->nbPulses[i]));
 216
 217       if (q<=0) {
 218          q = -q;
 219          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);
 220       }
 221
 222       if (q != 0)
 223       {
 224          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
 225          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
 226          alg_quant(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, enc);
 227          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
 228       }
 229       for (j=B*eBands[i];j<B*eBands[i+1];j++)
 230          norm[j] = X[j] * n;
 231       //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q))/(B*(eBands[i+1]-eBands[i])));
 232       //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q)));
 233    }
 234    //printf ("\n");
 235    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
 236       X[i] = 0;
 237 }
 238
 239 /* Decoding of the residual */
 240 void unquant_bands(const CELTMode *m, float *X, float *P, ec_dec *dec)
 241 {
 242    int i, j, B;
 243    const int *eBands = m->eBands;
 244    B = m->nbMdctBlocks*m->nbChannels;
 245    float norm[B*eBands[m->nbEBands+1]];
 246
 247    for (i=0;i<m->nbEBands;i++)
 248    {
 249       int q;
 250       float theta, n;
 251       q = m->nbPulses[i];
 252       n = sqrt(B*(eBands[i+1]-eBands[i]));
 253       theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+abs(m->nbPulses[i]));
 254
 255       if (q<=0) {
 256          q = -q;
 257          intra_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], dec);
 258       }
 259
 260       if (q != 0)
 261       {
 262          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
 263          alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, dec);
 264          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
 265       }
 266       for (j=B*eBands[i];j<B*eBands[i+1];j++)
 267          norm[j] = X[j] * n;
 268    }
 269    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
 270       X[i] = 0;
 271 }