Pitch prediction stuff...
authorjmvalin <jmvalin@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Thu, 28 Feb 2002 21:12:08 +0000 (21:12 +0000)
committerjmvalin <jmvalin@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Thu, 28 Feb 2002 21:12:08 +0000 (21:12 +0000)
git-svn-id: http://svn.xiph.org/trunk/speex@3118 0101bb08-14d6-0310-b084-bc0e0c8e3800

libspeex/cb_search.c
libspeex/ltp.c
libspeex/ltp.h
libspeex/speex.c

index 85360ef..8ffc64d 100644 (file)
@@ -81,7 +81,7 @@ int   nsf                       /* number of samples in subframe */
     d = 0.0; e = 0.0;
     for(i=0; i<nsf; i++) {
       d += target[i]*resp[i];
-      e += resp[i]*resp[i]+1;
+      e += resp[i]*resp[i];
     }
     g = d/(e+.1);
     score = g*d;
index 7f1a754..9bcefec 100644 (file)
 #define abs(x) ((x)<0 ? -(x) : (x))
 
 
-/** Computes a 3-tap pitch predictor */
-int three_tap_ltp(float *x, int len, int start, int end, float *gain)
 
-   /*  x:     time-domain signal (note, x[-end] must be valid)
-       len:   length of the x signal
-       start: smallest pitch period possible
-       end:   largest pitch period possible
-       gain:  return value for the pitch predictor gain
-    */
+/** Finds the best quantized 3-tap pitch predictor by analysis by synthesis */
+float pitch_search_3tap_unquant(
+float target[],                 /* Target vector */
+float ak[],                     /* LPCs for this subframe */
+float awk1[],                   /* Weighted LPCs #1 for this subframe */
+float awk2[],                   /* Weighted LPCs #2 for this subframe */
+float exc[],                    /* Overlapping codebook */
+int   start,                    /* Smallest pitch value allowed */
+int   end,                      /* Largest pitch value allowed */
+float *gain,                    /* 3-tab gains of optimum entry */
+int   *pitch,                   /* Index of optimum entry */
+int   p,                        /* Number of LPC coeffs */
+int   nsf                       /* Number of samples in subframe */
+)
 {
-   int i, period, best_period=0;
-   float total, score[3]={0,0,0}, best_score=-1, corr[3]={0,0,0}, energy[3]={0,0,0};
-   float best_energy[3], best_corr[3], best_gain=0;
+   int i,j;
+   float tmp[3*nsf];
+   float *x[3];
+   float corr[3];
    float A[3][3];
-   /*Check all periods */
-   for (period = start; period <= end; period++)
-   {
-      corr[0]=energy[0]=0;
-      for (i=0;i<len;i++)
-      {
-         corr[0] += x[i]*x[i-period];
-         energy[0] += x[i-period]*x[i-period];
-      }
-      score[0] = corr[0]*corr[0]/(energy[0]+1);
-      /* Smooth the "correlation score" */
-      total = score[0]+2*score[1]+score[2];
-      if (total > best_score && period >= start+3 && period <= end-3)
-      {
-         best_score=total;
-         best_period=period;
-         for (i=0;i<3;i++)
-         {
-            best_corr[i]=corr[i];
-            best_energy[i]=energy[i];
-         }
-         best_gain = corr[1]/(energy[1]+1);
-      }
-      score[2]=score[1];
-      score[1]=score[0];
-      corr[2]=corr[1];
-      corr[1]=corr[0];
-      energy[2]=energy[1];
-      energy[1]=energy[0];
-   }
-   /* build the correlation matrix */
-   A[0][0]=best_energy[0]+1;
-   A[1][1]=best_energy[1]+1;
-   A[2][2]=best_energy[2]+1;
-   A[0][1]=0;
-   A[0][2]=0;
-   for (i=0;i<len;i++)
+
+   x[0]=tmp;
+   x[1]=tmp+nsf;
+   x[2]=tmp+2*nsf;
+
+   /* Perform closed-loop 1-tap search*/
+   overlap_cb_search(target, ak, awk1, awk2,
+                     &exc[-end], end-start+1, gain, pitch, p,
+                     nsf);
+   /* Real pitch value */
+   *pitch=end-*pitch;
+   
+   
+   for (i=0;i<3;i++)
    {
-      A[0][1] += x[i-best_period+1]*x[i-best_period];
-      A[0][2] += x[i-best_period+2]*x[i-best_period];
+      residue_zero(&exc[-*pitch-1+i],awk1,x[i],nsf,p);
+      syn_filt_zero(x[i],ak,x[i],nsf,p);
+      syn_filt_zero(x[i],awk2,x[i],nsf,p);
    }
-   A[1][0]=A[0][1];
-   A[2][0]=A[0][2];
-   A[1][2]=A[0][1] - x[-best_period+1]*x[-best_period] 
-                   + x[len-best_period]*x[len-best_period+1];
-   A[2][1]=A[1][2];
-   /*for (i=0;i<3;i++)
-      printf ("%f %f %f\n", A[i][0], A[i][1], A[i][2]);
-      printf ("\n%f %f %f\n", best_corr[0], best_corr[1], best_corr[2]);*/
-   /*Solve the linear system to find gains*/
+
+   for (i=0;i<3;i++)
+      corr[i]=xcorr(x[i],target,nsf);
+   
+   for (i=0;i<3;i++)
+      for (j=0;j<=i;j++)
+         A[i][j]=A[j][i]=xcorr(x[i],x[j],nsf);
+   A[0][0]+=1;
+   A[1][1]+=1;
+   A[2][2]+=1;
    {
       float tmp=A[1][0]/A[0][0];
       for (i=0;i<3;i++)
          A[1][i] -= tmp*A[0][i];
-      best_corr[1] -= tmp*best_corr[0];
+      corr[1] -= tmp*corr[0];
 
       tmp=A[2][0]/A[0][0];
       for (i=0;i<3;i++)
          A[2][i] -= tmp*A[0][i];
-      best_corr[2] -= tmp*best_corr[0];
+      corr[2] -= tmp*corr[0];
       
       tmp=A[2][1]/A[1][1];
       A[2][2] -= tmp*A[1][2];
-      best_corr[2] -= tmp*best_corr[1];
+      corr[2] -= tmp*corr[1];
 
-      best_corr[2] /= A[2][2];
-      best_corr[1] = (best_corr[1] - A[1][2]*best_corr[2])/A[1][1];
-      best_corr[0] = (best_corr[0] - A[0][2]*best_corr[2] - A[0][1]*best_corr[1])/A[0][0];
+      corr[2] /= A[2][2];
+      corr[1] = (corr[1] - A[1][2]*corr[2])/A[1][1];
+      corr[0] = (corr[0] - A[0][2]*corr[2] - A[0][1]*corr[1])/A[0][0];
       /*printf ("\n%f %f %f\n", best_corr[0], best_corr[1], best_corr[2]);*/
 
    }
    /* Put gains in right order */
-   gain[0]=best_corr[2];gain[1]=best_corr[1];gain[2]=best_corr[0];
-   
-   /* Check that 3-tap predictor "makes sense" */
-   /*if (!((abs(gain[0]) + abs(gain[1]) + abs(gain[2])<2) 
-         && abs(gain[0]+gain[1]+gain[2]) < 1.2 
-         && abs(gain[0]+gain[1]+gain[2]) > -.2 ))
-   {
-      gain[0]=0;gain[1]=best_gain;gain[2]=0;
-      if (best_gain > 1.2)
-         gain[1]=1.2;
-      else if (best_gain<-.2)
-         gain[1]=-.2;
-      else 
-         gain[1]=best_gain;
-         }*/
-   return best_period-2;
-}
-
+   gain[0]=corr[2];gain[1]=corr[1];gain[2]=corr[0];
 
-/** Finds the best 3-tap pitch predictor from a codebook*/
-int ltp_closed_loop(float *x, int len, int start, int end, float *gain)
-
-   /*  x:     time-domain signal (note, x[-end] must be valid)
-       len:   length of the x signal
-       start: smallest pitch period possible
-       end:   largest pitch period possible
-       gain:  return value for the pitch predictor gain
-    */
-{
-   int i, period, best_period=0;
-   float total, score[3]={0,0,0}, best_score=-1, corr[3]={0,0,0}, energy[3]={0,0,0};
-   float best_energy[3], best_corr[3], best_gain=0;
-   float A[3][3];
-   /*Check all periods */
-   for (period = start; period <= end; period++)
-   {
-      corr[0]=energy[0]=0;
-      for (i=0;i<len;i++)
-      {
-         corr[0] += x[i]*x[i-period];
-         energy[0] += x[i-period]*x[i-period];
-      }
-      score[0] = corr[0]*corr[0]/(energy[0]+1);
-      /* Smooth the "correlation score" */
-      total = score[0]+2*score[1]+score[2];
-      if (total > best_score && period >= start+3 && period <= end-3)
-      {
-         best_score=total;
-         best_period=period;
-         for (i=0;i<3;i++)
-         {
-            best_corr[i]=corr[i];
-            best_energy[i]=energy[i];
-         }
-         best_gain = corr[1]/(energy[1]+1);
-      }
-      score[2]=score[1];
-      score[1]=score[0];
-      corr[2]=corr[1];
-      corr[1]=corr[0];
-      energy[2]=energy[1];
-      energy[1]=energy[0];
-   }
-   /* build the correlation matrix */
-   A[0][0]=best_energy[0]+1;
-   A[1][1]=best_energy[1]+1;
-   A[2][2]=best_energy[2]+1;
-   A[0][1]=0;
-   A[0][2]=0;
-   for (i=0;i<len;i++)
-   {
-      A[0][1] += x[i-best_period+1]*x[i-best_period];
-      A[0][2] += x[i-best_period+2]*x[i-best_period];
-   }
-   A[1][0]=A[0][1];
-   A[2][0]=A[0][2];
-   A[1][2]=A[0][1] - x[-best_period+1]*x[-best_period] 
-                   + x[len-best_period]*x[len-best_period+1];
-   A[2][1]=A[1][2];
+   --*pitch;
 
    {
-      int j;
-      float C[9];
-      float *ptr=gain_cdbk_nb;
-      int best_cdbk=0;
-      float best_sum=0;
-      C[0]=best_corr[2];
-      C[1]=best_corr[1];
-      C[2]=best_corr[0];
-      C[3]=A[1][2];
-      C[4]=A[0][1];
-      C[5]=A[0][2];
-      C[6]=A[2][2];
-      C[7]=A[1][1];
-      C[8]=A[0][0];
-      
-      for (i=0;i<127;i++)
-      {
-         float sum=0;
-         ptr = gain_cdbk_nb+12*i;
-         for (j=0;j<9;j++)
-            sum+=C[j]*ptr[j+3];
-         if (sum>best_sum || i==0)
-         {
-            best_sum=sum;
-            best_cdbk=i;
-         }
-      }
-      gain[0] = gain_cdbk_nb[best_cdbk*12];
-      gain[1] = gain_cdbk_nb[best_cdbk*12+1];
-      gain[2] = gain_cdbk_nb[best_cdbk*12+2];
+      float tmp1=0,tmp2=0;
+      for (i=0;i<nsf;i++)
+         tmp1+=target[i]*target[i];
+      for (i=0;i<nsf;i++)
+         tmp2+=(target[i]-gain[2]*x[0][i]-gain[1]*x[1][i]-gain[0]*x[2][i])
+         * (target[i]-gain[2]*x[0][i]-gain[1]*x[1][i]-gain[0]*x[2][i]);
+      printf ("prediction gain = %f\n",tmp1/(tmp2+1));
+      return tmp1/(tmp2+1);
    }
-   return best_period-2;
 }
 
 
 /** Finds the best quantized 3-tap pitch predictor by analysis by synthesis */
-void pitch_search_3tap(
+float pitch_search_3tap(
 float target[],                 /* Target vector */
 float ak[],                     /* LPCs for this subframe */
 float awk1[],                   /* Weighted LPCs #1 for this subframe */
@@ -238,6 +125,7 @@ int   start,                    /* Smallest pitch value allowed */
 int   end,                      /* Largest pitch value allowed */
 float *gain,                    /* 3-tab gains of optimum entry */
 int   *pitch,                   /* Index of optimum entry */
+int   *gain_index,              /* Index of optimum gain */
 int   p,                        /* Number of LPC coeffs */
 int   nsf                       /* Number of samples in subframe */
 )
@@ -305,6 +193,8 @@ int   nsf                       /* Number of samples in subframe */
       gain[0] = gain_cdbk_nb[best_cdbk*12];
       gain[1] = gain_cdbk_nb[best_cdbk*12+1];
       gain[2] = gain_cdbk_nb[best_cdbk*12+2];
+      *gain_index=best_cdbk;
+
    }
    --*pitch;
 
index 863e9a4..8681134 100644 (file)
 extern float gain_cdbk_nb[];
 
 
-/** Computes a 3-tap pitch predictor */
-int three_tap_ltp(float *x, int len, int start, int end, float *gain);
-
-/** Finds the best 3-tap pitch predictor from a codebook*/
-int ltp_closed_loop(float *x, int len, int start, int end, float *gain);
+/** Finds the best quantized 3-tap pitch predictor by analysis by synthesis */
+float pitch_search_3tap(
+float target[],                 /* Target vector */
+float ak[],                     /* LPCs for this subframe */
+float awk1[],                   /* Weighted LPCs #1 for this subframe */
+float awk2[],                   /* Weighted LPCs #2 for this subframe */
+float exc[],                    /* Overlapping codebook */
+int   start,                    /* Smallest pitch value allowed */
+int   end,                      /* Largest pitch value allowed */
+float *gain,                    /* 3-tab gains of optimum entry */
+int   *pitch,                   /* Best pitch delay */
+int   *gain_index,              /* Index of optimum gain */
+int   p,                        /* Number of LPC coeffs */
+int   nsf                       /* Number of samples in subframe */
+);
 
 /** Finds the best quantized 3-tap pitch predictor by analysis by synthesis */
-void pitch_search_3tap(
+float pitch_search_3tap_unquant(
 float target[],                 /* Target vector */
 float ak[],                     /* LPCs for this subframe */
 float awk1[],                   /* Weighted LPCs #1 for this subframe */
index acee9ac..54bf800 100644 (file)
@@ -203,7 +203,7 @@ void encode(EncState *st, float *in, int *outSize, void *bits)
    for (sub=0;sub<st->nbSubframes;sub++)
    {
       float tmp, tmp1,tmp2,gain[3];
-      int pitch, offset;
+      int pitch, offset, pitch_gain_index;
       float *sp, *sw, *res, *exc, *target;
       
       /* Offset relative to start of frame */
@@ -278,7 +278,7 @@ void encode(EncState *st, float *in, int *outSize, void *bits)
 
       /* Perform adaptive codebook search (3-tap pitch predictor) */
       pitch_search_3tap(target, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2,
-                        exc, 20, 147, &gain[0], &pitch, st->lpcSize,
+                        exc, 20, 147, &gain[0], &pitch, &pitch_gain_index, st->lpcSize,
                         st->subframeSize);
       for (i=0;i<st->subframeSize;i++)
         exc[i]=gain[0]*exc[i-pitch]+gain[1]*exc[i-pitch-1]+gain[2]*exc[i-pitch-2];
@@ -308,6 +308,34 @@ void encode(EncState *st, float *in, int *outSize, void *bits)
 
 
 #else
+
+#if 0 /* Code to calculate the exact excitation after pitch prediction  */
+      for (i=0;i<st->subframeSize;i++)
+         st->buf2[i]=target[i];
+      pitch_search_3tap(target, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2,
+                        exc, 20, 147, &gain[0], &pitch, &pitch_gain_index, st->lpcSize,
+                        st->subframeSize);
+      for (i=0;i<st->subframeSize;i++)
+        exc[i]=gain[0]*exc[i-pitch]+gain[1]*exc[i-pitch-1]+gain[2]*exc[i-pitch-2];
+      printf ("3-tap pitch = %d, gains = [%f %f %f]\n",pitch, gain[0], gain[1], gain[2]);
+
+      /* Update target for adaptive codebook contribution */
+      residue_zero(exc, st->bw_lpc1, res, st->subframeSize, st->lpcSize);
+      syn_filt_zero(res, st->interp_qlpc, res, st->subframeSize, st->lpcSize);
+      syn_filt_zero(res, st->bw_lpc2, res, st->subframeSize, st->lpcSize);
+      for (i=0;i<st->subframeSize;i++)
+        target[i]-=res[i];
+      syn_filt_zero(target, st->bw_lpc1, res, st->subframeSize, st->lpcSize);
+      residue_zero(res, st->interp_qlpc, exc, st->subframeSize, st->lpcSize);
+      residue_zero(exc, st->bw_lpc2, exc, st->subframeSize, st->lpcSize);
+      for (i=0;i<st->subframeSize;i++)
+         printf ("%f ", exc[i]);
+      printf ("\n");
+      
+      for (i=0;i<st->subframeSize;i++)
+         target[i]=st->buf2[i];
+#endif
+
       /* We're cheating to get perfect reconstruction */
       syn_filt_zero(target, st->bw_lpc1, res, st->subframeSize, st->lpcSize);
       residue_zero(res, st->interp_qlpc, exc, st->subframeSize, st->lpcSize);