Some pitch prediction work.

libcelt/bands.c

@@ -33,7 +33,11 @@
 #include "bands.h"
 
 const int qbank[NBANDS+2] =   {0, 2, 4, 6, 8, 12, 16, 20, 24, 28, 36, 44, 52, 68, 84, 116, 128};
+int qpulses[] = {4, 5, 4, 4, 3,  3,  3,  3,  3,  4,  4,  4,  0,  0,  0}; //c: 134 bits
+#define WAVEFORM_END 52
 
+/* Start frequency of each band */
+int pbank[] = {0, 4, 8, 12, 20, WAVEFORM_END, 128};
 
 /* Compute the energy in each of the bands */
 void compute_bands(float *X, int B, float *bank)
@@ -56,7 +60,7 @@ void normalise_bands(float *X, int B, float *bank)
    for (i=0;i<NBANDS;i++)
    {
       int j;
-      float x = 1.f/bank[i];
+      float x = 1.f/(1e-10+bank[i]);
       for (j=B*qbank[i];j<B*qbank[i+1];j++)
          X[j] *= x;
    }
@@ -80,6 +84,61 @@ void denormalise_bands(float *X, int B, float *bank)
 }
 
 
+/* Compute the best gain for each "pitch band" */
+void compute_pitch_gain(float *X, int B, float *P, float *gains, float *bank)
+{
+   int i;
+   float w[B*qbank[NBANDS]];
+   for (i=0;i<NBANDS;i++)
+   {
+      int j;
+      for (j=B*qbank[i];j<B*qbank[i+1];j++)
+         w[j] = bank[i];
+   }
+
+   
+   for (i=0;i<PBANDS;i++)
+   {
+      float Sxy=0;
+      float Sxx = 0;
+      int j;
+      float gain;
+      for (j=B*pbank[i];j<B*pbank[i+1];j++)
+      {
+         Sxy += X[j]*P[j]*w[j];
+         Sxx += X[j]*X[j]*w[j];
+      }
+      gain = Sxy/(1e-10+Sxx);
+      //gain = Sxy/(2*(pbank[i+1]-pbank[i]));
+      //if (i<3)
+      //gain *= 1+.02*gain;
+      if (gain > .90)
+         gain = .90;
+      if (gain < 0.0)
+         gain = 0.0;
+      
+      gains[i] = gain;
+   }
+   for (i=B*pbank[PBANDS];i<B*pbank[PBANDS+1];i++)
+      P[i] = 0;
+}
+
+/* Apply the (quantised) gain to each "pitch band" */
+void pitch_quant_bands(float *X, int B, float *P, float *gains)
+{
+   int i;
+   for (i=0;i<PBANDS;i++)
+   {
+      int j;
+      for (j=B*pbank[i];j<B*pbank[i+1];j++)
+         P[j] *= gains[i];
+      //printf ("%f ", gain);
+   }
+   for (i=B*pbank[PBANDS];i<B*pbank[PBANDS+1];i++)
+      P[i] = 0;
+}
+
+
 /* Scales the pulse-codebook entry in each band such that unit-energy is conserved when 
    adding the pitch */
 void pitch_renormalise_bands(float *X, int B, float *P)
@@ -105,7 +164,7 @@ void pitch_renormalise_bands(float *X, int B, float *P)
       Rxx = 0;
       for (j=B*qbank[i];j<B*qbank[i+1];j++)
       {
-         X[j*2-1] = P[j]+gain1*X[j];
+         X[j] = P[j]+gain1*X[j];
          Rxx += X[j]*X[j];
       }
    }

libcelt/bands.h

@@ -32,7 +32,10 @@
 #ifndef BANDS_H
 #define BANDS_H
 
+/* Number of constant-energy bands */
 #define NBANDS 15
+/* Number of bands only for the pitch prediction */
+#define PBANDS 5
 
 void compute_bands(float *X, int B, float *bands);
 
@@ -40,6 +43,10 @@ void normalise_bands(float *X, int B, float *bands);
 
 void denormalise_bands(float *X, int B, float *bands);
 
+void compute_pitch_gain(float *X, int B, float *P, float *gains, float *bank);
+
+void pitch_quant_bands(float *X, int B, float *P, float *gains);
+
 void pitch_renormalise_bands(float *X, int B, float *P);
 
 #endif /* BANDS_H */

libcelt/celt.c

@@ -55,8 +55,6 @@ struct CELTState_ {
    float *mdct_overlap;
    float *out_mem;
 
-   float *bandE;
-
 };
 
 
@@ -78,7 +76,6 @@ CELTState *celt_encoder_new(int blockSize, int blocksPerFrame)
    st->in_mem = celt_alloc(N*sizeof(float));
    st->mdct_overlap = celt_alloc(N*sizeof(float));
    st->out_mem = celt_alloc(MAX_PERIOD*sizeof(float));
-   st->bandE = celt_alloc(NBANDS*sizeof(float));
    for (i=0;i<N;i++)
       st->window[i] = st->window[2*N-i-1] = sin(.5*M_PI* sin(.5*M_PI*(i+.5)/N) * sin(.5*M_PI*(i+.5)/N));
    return st;
@@ -99,8 +96,6 @@ void celt_encoder_destroy(CELTState *st)
    celt_free(st->mdct_overlap);
    celt_free(st->out_mem);
    
-   celt_free(st->bandE);
-
    celt_free(st);
 }
 
@@ -131,6 +126,9 @@ int celt_encode(CELTState *st, short *pcm)
    
    float X[B*N]; /**< Interleaved signal MDCTs */
    float P[B*N]; /**< Interleaved pitch MDCTs*/
+   float bandEp[NBANDS];
+   float bandE[NBANDS];
+   float gains[PBANDS];
    int pitch_index;
    
    /* FIXME: Add preemphasis */
@@ -164,18 +162,34 @@ int celt_encode(CELTState *st, short *pcm)
    printf ("\n");*/
    
    /* Band normalisation */
-   compute_bands(X, B, st->bandE);
+   compute_bands(X, B, bandE);
-   normalise_bands(X, B, st->bandE);
+   //for (i=0;i<NBANDS;i++) printf("%f ",bandE[i]);printf("\n"); 
+   normalise_bands(X, B, bandE);
    
-   compute_bands(P, B, st->bandE);
+   compute_bands(P, B, bandEp);
-   normalise_bands(P, B, st->bandE);
+   normalise_bands(P, B, bandEp);
 
    /* Pitch prediction */
+   compute_pitch_gain(X, B, P, gains, bandE);
+   //quantise_pitch(gains, PBANDS);
+   pitch_quant_bands(X, B, P, gains);
+   
+   for (i=0;i<B*N;i++) printf("%f ",P[i]);printf("\n");
+   /* Subtract the pitch prediction from the signal to encode */
+   for (i=0;i<B*N;i++)
+      X[i] -= P[i];
 
    /* Residual quantisation */
+   if (1) {
+      float tmpE[NBANDS];
+      compute_bands(X, B, tmpE);
+      normalise_bands(X, B, tmpE);
+      pitch_renormalise_bands(X, B, P);
+   }
+   //quant_bands(X, P);
    
    /* Synthesis */
-   denormalise_bands(X, B, st->bandE);
+   denormalise_bands(X, B, bandE);
 
    CELT_MOVE(st->out_mem, st->out_mem+B*N, MAX_PERIOD-B*N);
    /* Compute inverse MDCTs */