Added 3-tap pitch predictor by analysis by synthesis
authorjmvalin <jmvalin@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Wed, 27 Feb 2002 19:32:58 +0000 (19:32 +0000)
committerjmvalin <jmvalin@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Wed, 27 Feb 2002 19:32:58 +0000 (19:32 +0000)
git-svn-id: http://svn.xiph.org/trunk/speex@3111 0101bb08-14d6-0310-b084-bc0e0c8e3800

configure.in
libspeex/Makefile.am
libspeex/cb_search.c
libspeex/filters.c
libspeex/filters.h
libspeex/ltp.c
libspeex/ltp.h
libspeex/speex.c
libspeex/stack_alloc.h [new file with mode: 0644]

index ce38e2f..e5006ca 100644 (file)
@@ -32,6 +32,7 @@ dnl AC_DISABLE_STATIC
 AM_PROG_LIBTOOL
 
 AC_CHECK_LIB(m, sin)
+AC_FUNC_ALLOCA
 dnl AC_DEFINE(VERSION, ${VERSION})
 dnl AC_DEFINE_UNQUOTED(VERSION, ${VERSION})
 
index 5400829..3cf81f9 100644 (file)
@@ -1,6 +1,6 @@
 ## Process this file with automake to produce Makefile.in. -*-Makefile-*-
 
-# $Id: Makefile.am,v 1.4 2002/02/25 08:46:21 jmvalin Exp $
+# $Id: Makefile.am,v 1.5 2002/02/27 19:32:58 jmvalin Exp $
 
 # Disable automatic dependency tracking if using other tools than gcc and gmake
 #AUTOMAKE_OPTIONS = no-dependencies
@@ -26,7 +26,8 @@ noinst_HEADERS = lsp.h \
        ltp.h \
        quant_lsp.h \
        cb_search.h \
-       filters.h
+       filters.h \
+       stack_alloc.h
        
        
 libspeex_la_LDFLAGS = -release $(LT_RELEASE)
index ad6ef4c..9f9e6dd 100644 (file)
@@ -6,7 +6,7 @@
     COMPANY.....: Voicetronix
     DATE CREATED: 19/2/02
 
-    General gain-shape codebbok search.
+    General gain-shape codebook search.
 
 \*-----------------------------------------------------------------------*/
 
index f6baf53..fbe8790 100644 (file)
@@ -95,3 +95,12 @@ void residue_mem(float *x, float *a, float *y, int N, int ord, float *mem)
    for (i=0;i<ord;i++)
       mem[i]=x[N-i-1];
 }
+
+float xcorr(float *x, float *y, int len)
+{
+   int i;
+   float sum=0;
+   for (i=0;i<len;i++)
+      sum += x[i]*y[i];
+   return sum;
+}
index dbbb017..14337b8 100644 (file)
@@ -38,4 +38,7 @@ void residue_zero(float *x, float *a, float *y, int N, int ord);
 /* Analysis (FIR) filter using memory */
 void residue_mem(float *x, float *a, float *y, int N, int ord, float *mem);
 
+/* Cross correlation */
+float xcorr(float *x, float *y, int len);
+
 #endif
index aa65941..7f1a754 100644 (file)
 */
 
 #include "ltp.h"
-
+#include "cb_search.h"
+#include "stack_alloc.h"
+#include "filters.h"
 
 #define abs(x) ((x)<0 ? -(x) : (x))
 
-/** Computes the open-loop pitch prediction. Returns pitch period and pitch gain */
-int open_loop_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
-    */
-{
-   int i, period, best_period=0;
-   float score, best_score=-1, corr, energy;
-   
-   for (period = start; period <= end; period++)
-   {
-      corr=energy=0;
-      for (i=0;i<len;i++)
-      {
-         corr += x[i]*x[i-period];
-         energy += x[i-period]*x[i-period];
-      }
-      score = corr*corr/(energy+1);
-      if (score > best_score)
-      {
-         best_score=score;
-         best_period=period;
-         *gain = corr/(energy+1);
-      }
-      
-   }
-   return best_period;
-}
 
 /** Computes a 3-tap pitch predictor */
 int three_tap_ltp(float *x, int len, int start, int end, float *gain)
@@ -256,24 +227,94 @@ int ltp_closed_loop(float *x, int len, int start, int end, float *gain)
 }
 
 
-
-/** In place 3-tap pitch predictor (FIR)*/
-void predictor_three_tap(float *x, int len, int period, float *gain)
+/** Finds the best quantized 3-tap pitch predictor by analysis by synthesis */
+void 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,                   /* Index of optimum entry */
+int   p,                        /* Number of LPC coeffs */
+int   nsf                       /* Number of samples in subframe */
+)
 {
-   int i;
-   for (i=len-1;i>=0;i--)
-   /*for (i=0;i<len;i++)*/
+   int i,j;
+   float tmp[3*nsf];
+   float *x[3];
+   float corr[3];
+   float A[3][3];
+
+   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++)
    {
-      x[i] -= gain[0]*x[i-period] + gain[1]*x[i-period-1] + gain[2]*x[i-period-2];
+      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);
    }
-}
 
-/** In place 3-tap inverse pitch predictor (IIR)*/
-void inverse_three_tap(float *x, int len, int period, float *gain)
-{
-   int i;
-   for (i=0;i<len;i++)
+   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);
+   
+   {
+      int j;
+      float C[9];
+      float *ptr=gain_cdbk_nb;
+      int best_cdbk=0;
+      float best_sum=0;
+      C[0]=corr[2];
+      C[1]=corr[1];
+      C[2]=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];
+   }
+   --*pitch;
+
    {
-      x[i] += gain[0]*x[i-period] + gain[1]*x[i-period-1] + gain[2]*x[i-period-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));
    }
 }
index d2bb558..863e9a4 100644 (file)
@@ -20,9 +20,6 @@
 
 extern float gain_cdbk_nb[];
 
-/** Computes the open-loop pitch prediction. Returns pitch period and pitch gain */
-int open_loop_ltp(float *x, int len, int start, int end, float *gain);
-
 
 /** Computes a 3-tap pitch predictor */
 int three_tap_ltp(float *x, int len, int start, int end, float *gain);
@@ -30,9 +27,17 @@ 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);
 
-/** In place 3-tap pitch predictor (FIR)*/
-void predictor_three_tap(float *x, int len, int period, float *gain);
-
-
-/** In place 3-tap inverse pitch predictor (IIR)*/
-void inverse_three_tap(float *x, int len, int period, float *gain);
+/** Finds the best quantized 3-tap pitch predictor by analysis by synthesis */
+void 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,                   /* Index of optimum entry */
+int   p,                        /* Number of LPC coeffs */
+int   nsf                       /* Number of samples in subframe */
+);
index ffdd5c1..b2849a0 100644 (file)
@@ -268,14 +268,16 @@ void encode(EncState *st, float *in, int *outSize, void *bits)
       /* Compute target signal */
       for (i=0;i<st->subframeSize;i++)
          target[i]=sw[i]-res[i];
-#if 1
-      /* Perform adaptive codebook search (pitch prediction) */
-      overlap_cb_search(target, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2,
-                        &exc[-120], 100, &gain[0], &pitch, st->lpcSize,
+
+#if 1 /*If set to 0, we compute the excitation directly from the target, i.e. we're cheating */
+
+      /* 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,
                         st->subframeSize);
-      printf ("gain = %f, pitch = %d\n",gain[0], 120-pitch);
       for (i=0;i<st->subframeSize;i++)
-        exc[i]=gain[0]*exc[i-120+pitch];
+        exc[i]=gain[0]*exc[i-pitch]+gain[1]*exc[i-pitch-1]+gain[2]*exc[i-pitch-2];
+      printf ("3tap 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);
@@ -298,7 +300,10 @@ void encode(EncState *st, float *in, int *outSize, void *bits)
       syn_filt_zero(res, st->bw_lpc2, res, st->subframeSize, st->lpcSize);
       for (i=0;i<st->subframeSize;i++)
          target[i]-=gain[0]*res[i];
+
+
 #else
+      /* 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);
       residue_zero(exc, st->bw_lpc2, exc, st->subframeSize, st->lpcSize);
@@ -308,11 +313,10 @@ void encode(EncState *st, float *in, int *outSize, void *bits)
       syn_filt_mem(exc, st->interp_qlpc, sp, st->subframeSize, st->lpcSize, st->mem5);
 
       /* Compute weighted signal again, from synthesized speech (not sure it's the right thing) */
-#if 0
-      residue_mem(sp, st->bw_lpc1, sw, st->subframeSize, st->lpcSize, st->mem6);
+      residue_mem(sp, st->bw_lpc1, sw, st->subframeSize, st->lpcSize, st->mem1);
+      syn_filt_mem(sw, st->bw_lpc2, sw, st->subframeSize, st->lpcSize, st->mem2);
       for (i=0;i<st->lpcSize;i++)
         st->mem3[i]=sw[st->subframeSize-i-1];
-#endif
    }
 
    /* Store the LSPs for interpolation in the next frame */
diff --git a/libspeex/stack_alloc.h b/libspeex/stack_alloc.h
new file mode 100644 (file)
index 0000000..d58015c
--- /dev/null
@@ -0,0 +1,53 @@
+/* Copyright (C) 2002 Jean-Marc Valin 
+   File: stack_alloc.h
+   
+   Temporary memory allocation on stack
+
+   This library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+   
+   This library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+   
+   You should have received a copy of the GNU Lesser General Public
+   License along with this library; if not, write to the Free Software
+   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*/
+
+#ifndef STACK_ALLOC_H
+#define STACK_ALLOC_H
+
+
+#ifdef __GNUC__
+#define VAR_ARRAY
+#endif
+
+
+#if defined (VAR_ARRAY)  /* Prefered method is variable-size arrays is supported */
+
+#define DYN_VEC(type, num, var) type var[num];
+
+#elif defined (HAVE_ALLOCA_H)  /* Second best: alloca */
+
+#include <alloca.h>
+
+#define DYN_VEC(type, num, var) type *var=(type*)alloca((num)*sizeof(type));
+
+#elif defined WIN32  /* On Win32, it's _alloca */
+
+#include <malloc.h>
+#define DYN_VEC(type, num, var) type *var=(type*)_alloca((num)*sizeof(type));
+
+#else  /* When all else fails, allocate on the heap (but it's going to be slow) */
+
+#error Cannot allocate memory on stack
+
+#endif
+
+
+
+#endif