Working on some stability issues (appears to be solved by making the pitch
authorJean-Marc Valin <jean-marc.valin@usherbrooke.ca>
Tue, 11 Dec 2007 13:45:15 +0000 (00:45 +1100)
committerJean-Marc Valin <jean-marc.valin@usherbrooke.ca>
Tue, 11 Dec 2007 13:45:15 +0000 (00:45 +1100)
projection less aggressive). Also, fixed a 64-bit overflow in the stereo mode
and added a "band rotation" function.

libcelt/bands.c
libcelt/bands.h
libcelt/celt.c
libcelt/modes.c
libcelt/vq.c
libcelt/vq.h

index 75fd66e..b1a7580 100644 (file)
 #include "vq.h"
 #include "cwrs.h"
 
+/* Applies a series of rotations so that pulses are spread like a two-sided expo
+nential */
+static void exp_rotation(float *X, int len, float theta, int dir)
+{
+   int i;
+   float c, s;
+   c = cos(theta);
+   s = sin(theta);
+   if (dir > 0)
+   {
+      for (i=0;i<(len/2)-1;i++)
+      {
+         float x1, x2;
+         x1 = X[2*i];
+         x2 = X[2*i+2];
+         X[2*i] = c*x1 - s*x2;
+         X[2*i+2] = c*x2 + s*x1;
+         
+         x1 = X[2*i+1];
+         x2 = X[2*i+3];
+         X[2*i+1] = c*x1 - s*x2;
+         X[2*i+3] = c*x2 + s*x1;
+      }
+      for (i=(len/2)-3;i>=0;i--)
+      {
+         float x1, x2;
+         x1 = X[2*i];
+         x2 = X[2*i+2];
+         X[2*i] = c*x1 - s*x2;
+         X[2*i+2] = c*x2 + s*x1;
+         
+         x1 = X[2*i+1];
+         x2 = X[2*i+3];
+         X[2*i+1] = c*x1 - s*x2;
+         X[2*i+3] = c*x2 + s*x1;
+      }
+
+   } else {
+      for (i=0;i<(len/2)-2;i++)
+      {
+         float x1, x2;
+         x1 = X[2*i];
+         x2 = X[2*i+2];
+         X[2*i] = c*x1 + s*x2;
+         X[2*i+2] = c*x2 - s*x1;
+         
+         x1 = X[2*i+1];
+         x2 = X[2*i+3];
+         X[2*i+1] = c*x1 + s*x2;
+         X[2*i+3] = c*x2 - s*x1;
+      }
+      
+      for (i=(len/2)-2;i>=0;i--)
+      {
+         float x1, x2;
+         x1 = X[2*i];
+         x2 = X[2*i+2];
+         X[2*i] = c*x1 + s*x2;
+         X[2*i+2] = c*x2 - s*x1;
+         
+         x1 = X[2*i+1];
+         x2 = X[2*i+3];
+         X[2*i+1] = c*x1 + s*x2;
+         X[2*i+3] = c*x2 - s*x1;
+      }
+   }
+}
+
 
 /* Compute the energy in each of the bands */
 void compute_band_energies(const CELTMode *m, float *X, float *bank)
@@ -159,16 +227,18 @@ void quant_bands(const CELTMode *m, float *X, float *P, ec_enc *enc)
       q = m->nbPulses[i];
       if (q>0) {
          float n = sqrt(B*(eBands[i+1]-eBands[i]));
-         alg_quant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], enc);
+         alg_quant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, enc);
          for (j=B*eBands[i];j<B*eBands[i+1];j++)
             norm[j] = X[j] * n;
-         //printf ("%f ", log2(ncwrs(B*(eBands[i+1]-eBands[i]), q))/(B*(eBands[i+1]-eBands[i])));
+         //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q))/(B*(eBands[i+1]-eBands[i])));
+         //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q)));
       } else {
          float n = sqrt(B*(eBands[i+1]-eBands[i]));
          copy_quant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), -q, norm, B, eBands[i], enc);
          for (j=B*eBands[i];j<B*eBands[i+1];j++)
             norm[j] = X[j] * n;
-         //printf ("%f ", (1+log2(eBands[i]-(eBands[i+1]-eBands[i]))+log2(ncwrs(B*(eBands[i+1]-eBands[i]), -q)))/(B*(eBands[i+1]-eBands[i])));
+         //printf ("%f ", (1+log2(eBands[i]-(eBands[i+1]-eBands[i]))+log2(ncwrs64(B*(eBands[i+1]-eBands[i]), -q)))/(B*(eBands[i+1]-eBands[i])));
+         //printf ("%f ", (1+log2(eBands[i]-(eBands[i+1]-eBands[i]))+log2(ncwrs64(B*(eBands[i+1]-eBands[i]), -q))));
       }
    }
    //printf ("\n");
@@ -189,7 +259,7 @@ void unquant_bands(const CELTMode *m, float *X, float *P, ec_dec *dec)
       q = m->nbPulses[i];
       if (q>0) {
          float n = sqrt(B*(eBands[i+1]-eBands[i]));
-         alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], dec);
+         alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, dec);
          for (j=B*eBands[i];j<B*eBands[i+1];j++)
             norm[j] = X[j] * n;
       } else {
@@ -204,3 +274,17 @@ void unquant_bands(const CELTMode *m, float *X, float *P, ec_dec *dec)
    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
       X[i] = 0;
 }
+
+void band_rotation(const CELTMode *m, float *X, int dir)
+{
+   int i, B;
+   const int *eBands = m->eBands;
+   B = m->nbMdctBlocks*m->nbChannels;
+   for (i=0;i<m->nbEBands;i++)
+   {
+      float theta;
+      theta = pow(.1f,1.f*abs(m->nbPulses[i])/(B*(eBands[i+1]-eBands[i])));
+      exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, dir);
+   }
+   //printf ("\n");
+}
index 922bea9..01d69fc 100644 (file)
@@ -51,4 +51,6 @@ void quant_bands(const CELTMode *m, float *X, float *P, ec_enc *enc);
 
 void unquant_bands(const CELTMode *m, float *X, float *P, ec_dec *dec);
 
+void band_rotation(const CELTMode *m, float *X, int dir);
+
 #endif /* BANDS_H */
index 19ef1f8..0bcf53b 100644 (file)
@@ -266,6 +266,9 @@ int celt_encode(CELTEncoder *st, short *pcm)
       normalise_bands(st->mode, P, bandEp);
    }
    
+   band_rotation(st->mode, X, -1);
+   band_rotation(st->mode, P, -1);
+   
    quant_energy(st->mode, bandE, st->oldBandE, &st->enc);
    
    /* Pitch prediction */
@@ -295,6 +298,8 @@ int celt_encode(CELTEncoder *st, short *pcm)
       //printf ("\n");
    }
    
+   band_rotation(st->mode, X, 1);
+
    /* Synthesis */
    denormalise_bands(st->mode, X, bandE);
 
@@ -495,6 +500,7 @@ int celt_decode(CELTDecoder *st, char *data, int len, short *pcm)
       compute_band_energies(st->mode, P, bandEp);
       normalise_bands(st->mode, P, bandEp);
    }
+   band_rotation(st->mode, P, -1);
 
    /* Get the pitch gains */
    unquant_pitch(gains, st->mode->nbPBands, &dec);
@@ -505,6 +511,8 @@ int celt_decode(CELTDecoder *st, char *data, int len, short *pcm)
    /* Decode fixed codebook and merge with pitch */
    unquant_bands(st->mode, X, P, &dec);
 
+   band_rotation(st->mode, X, 1);
+   
    /* Synthesis */
    denormalise_bands(st->mode, X, bandE);
 
index 36c48de..44b1d60 100644 (file)
@@ -39,7 +39,7 @@ const int qbank1[NBANDS128+2] =   {0, 2, 4, 6, 8, 12, 16, 20, 24, 28, 36, 44, 52
 
 const int qpulses1[NBANDS128] =   {7, 5, 5, 5, 4,  5,  4,  5,  5,  4, -2, 0, 0, 0,  0};
 const int qpulses2[NBANDS128] =   {28,24,20,16,24,20, 18, 12, 10,  10,-7, -4, 0, 0,  0};
-const int qpulses2b[NBANDS128] =   {32,28,24,20,28,24, 22, 18, 16,  15,-12, -12, 12, 12,  0};
+const int qpulses2s[NBANDS128] =  {38,30,24,20,24,20, 18, 16, 14, 20,-20,-14, -8, -8,  -5};
 
 const int pbank1[PBANDS128+2] =   {0, 4, 8, 12, 20, PITCH_END128, 128};
 
@@ -94,7 +94,7 @@ const CELTMode mode3 = {
    qpulses2     /**< nbPulses */
 };
 
-/* Stereo mode (doesn't work yet) */
+/* Stereo mode around 120 kbps */
 const CELTMode mode4 = {
    256,         /**< frameSize */
    128,         /**< mdctSize */
@@ -107,7 +107,7 @@ const CELTMode mode4 = {
    
    qbank1,      /**< eBands */
    pbank1,      /**< pBands*/
-   qpulses1     /**< nbPulses */
+   qpulses2s     /**< nbPulses */
 };
 
 const CELTMode const *celt_mode1 = &mode1;
index 8469e29..4073587 100644 (file)
@@ -38,9 +38,9 @@
 /* Improved algebraic pulse-base quantiser. The signal x is replaced by the sum of the pitch 
    a combination of pulses such that its norm is still equal to 1. The only difference with 
    the quantiser above is that the search is more complete. */
-void alg_quant(float *x, int N, int K, float *p, ec_enc *enc)
+void alg_quant(float *x, int N, int K, float *p, float alpha, ec_enc *enc)
 {
-   int L = 5;
+   int L = 3;
    //float tata[200];
    float y[L][N];
    int iy[L][N];
@@ -55,7 +55,6 @@ void alg_quant(float *x, int N, int K, float *p, ec_enc *enc)
    float Rpp=0, Rxp=0;
    float gain[L];
    int maxL = 1;
-   float alpha = .9;
    
    for (j=0;j<N;j++)
       Rpp += p[j]*p[j];
@@ -188,9 +187,40 @@ void alg_quant(float *x, int N, int K, float *p, ec_enc *enc)
    //   printf ("%d ", iy[0][i]);
    pulse2comb(N, K, comb, signs, iy[0]); 
    ec_enc_uint64(enc,icwrs64(N, K, comb, signs),ncwrs64(N, K));
+   
+   /* Recompute the gain in one pass (to reduce errors) */
+   if (0) {
+      float Ryp=0;
+      float Rpp=0;
+      float Ryy=0;
+      float g=0;
+      for (i=0;i<N;i++)
+         Rpp += p[i]*p[i];
+      
+      for (i=0;i<N;i++)
+         Ryp += iy[0][i]*p[i];
+      
+      for (i=0;i<N;i++)
+         y[0][i] = iy[0][i] - alpha*Ryp*p[i];
+      
+      /* Recompute after the projection (I think it's right) */
+      Ryp = 0;
+      for (i=0;i<N;i++)
+         Ryp += y[0][i]*p[i];
+      
+      for (i=0;i<N;i++)
+         Ryy += y[0][i]*y[0][i];
+      
+      g = (sqrt(Ryp*Ryp + Ryy - Ryy*Rpp) - Ryp)/Ryy;
+        
+      for (i=0;i<N;i++)
+         x[i] = p[i] + g*y[0][i];
+      
+   }
+
 }
 
-static const float pg[5] = {1.f, .82f, .75f, 0.7f, 0.6f};
+static const float pg[5] = {1.f, .6f, .45f, 0.35f, 0.25f};
 
 /* Finds the right offset into Y and copy it */
 void copy_quant(float *x, int N, int K, float *Y, int B, int N0, ec_enc *enc)
@@ -257,11 +287,11 @@ void copy_quant(float *x, int N, int K, float *Y, int B, int N0, ec_enc *enc)
       E = .8/sqrt(E);
       for (j=0;j<N;j++)
          P[j] *= E;
-      alg_quant(x, N, K, P, enc);
+      alg_quant(x, N, K, P, 0, enc);
    }
 }
 
-void alg_unquant(float *x, int N, int K, float *p, ec_dec *dec)
+void alg_unquant(float *x, int N, int K, float *p, float alpha, ec_dec *dec)
 {
    int i;
    celt_uint64_t id;
@@ -269,7 +299,6 @@ void alg_unquant(float *x, int N, int K, float *p, ec_dec *dec)
    int signs[K];
    int iy[N];
    float y[N];
-   float alpha = .9;
    float Rpp=0, Ryp=0, Ryy=0;
    float g;
    
@@ -344,6 +373,6 @@ void copy_unquant(float *x, int N, int K, float *Y, int B, int N0, ec_dec *dec)
       E = .8/sqrt(E);
       for (j=0;j<N;j++)
          P[j] *= E;
-      alg_unquant(x, N, K, P, dec);
+      alg_unquant(x, N, K, P, 0, dec);
    }
 }
index 600c946..ae75f7e 100644 (file)
@@ -39,9 +39,9 @@
 /* Algebraic pulse-base quantiser. The signal x is replaced by the sum of the pitch 
    a combination of pulses such that its norm is still equal to 1. The only difference with 
    the quantiser above is that the search is more complete. */
-void alg_quant(float *x, int N, int K, float *p, ec_enc *enc);
+void alg_quant(float *x, int N, int K, float *p, float alpha, ec_enc *enc);
 
-void alg_unquant(float *x, int N, int K, float *p, ec_dec *dec);
+void alg_unquant(float *x, int N, int K, float *p, float alpha, ec_dec *dec);
 
 /* Finds the right offset into Y and copy it */
 void copy_quant(float *x, int N, int K, float *Y, int B, int N0, ec_enc *enc);