New post-filter based on a new way of moving poles in the LPC polynomial
authorjmvalin <jmvalin@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Wed, 17 Jul 2002 06:35:05 +0000 (06:35 +0000)
committerjmvalin <jmvalin@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Wed, 17 Jul 2002 06:35:05 +0000 (06:35 +0000)
git-svn-id: http://svn.xiph.org/trunk/speex@3644 0101bb08-14d6-0310-b084-bc0e0c8e3800

Speex.spec
configure.in
libspeex/Makefile.am
libspeex/filters.c
libspeex/filters.h
libspeex/modes.c
libspeex/post_filter.c
libspeex/roots.c [new file with mode: 0644]
libspeex/roots.h [new file with mode: 0644]

index 2051649..55e9e81 100644 (file)
@@ -1,5 +1,5 @@
 %define name     Speex
-%define ver      0.5.0
+%define ver      0.5.1
 %define rel      1
 
 Summary: An open-source, patent-free speech codec
index 3acfa46..4362580 100644 (file)
@@ -4,7 +4,7 @@ AC_INIT(libspeex/speex.h)
 
 SPEEX_MAJOR_VERSION=0
 SPEEX_MINOR_VERSION=5
-SPEEX_MICRO_VERSION=0
+SPEEX_MICRO_VERSION=1
 SPEEX_VERSION=$SPEEX_MAJOR_VERSION.$SPEEX_MINOR_VERSION.$SPEEX_MICRO_VERSION
 SPEEX_BINARY_AGE=0
 SPEEX_INTERFACE_AGE=0
index e223c1c..d439716 100644 (file)
@@ -1,6 +1,6 @@
 ## Process this file with automake to produce Makefile.in. -*-Makefile-*-
 
-# $Id: Makefile.am,v 1.36 2002/07/02 05:14:28 jmvalin Exp $
+# $Id: Makefile.am,v 1.37 2002/07/17 06:35:05 jmvalin Exp $
 
 # Disable automatic dependency tracking if using other tools than gcc and gmake
 #AUTOMAKE_OPTIONS = no-dependencies
@@ -32,7 +32,8 @@ libspeex_la_SOURCES = nb_celp.c \
        hexc_10_32_table.c \
        misc.c \
        speex_header.c \
-       post_filter.c
+       post_filter.c \
+       roots.c
 
 
 include_HEADERS =  speex.h \
@@ -52,7 +53,8 @@ noinst_HEADERS = lsp.h \
        sb_celp.h \
        vbr.h \
        post_filter.h \
-       misc.h
+       misc.h \
+       roots.h
        
        
 libspeex_la_LDFLAGS = -release $(LT_RELEASE)
index ba0e465..1c9143e 100644 (file)
 
 #include "filters.h"
 #include <stdio.h>
+#include "roots.h"
+#include <math.h>
 
 #define min(a,b) ((a) < (b) ? (a) : (b))
 
+#define MAX_ORD 20
+
 void print_vec(float *vec, int len, char *name)
 {
    int i;
@@ -42,6 +46,51 @@ void bw_lpc(float gamma, float *lpc_in, float *lpc_out, int order)
    }
 }
 
+void poly(float *re, float *im, float *p, int ord)
+{
+   int i,j;
+
+   float p_re[MAX_ORD], p_im[MAX_ORD];
+   for(i=0;i<ord+1;i++)
+      p_re[i]=p_im[i]=0;
+   p_re[0]=1;
+   for (i=0;i<ord;i++)
+   {
+      for (j=i;j>=0;j--)
+      {
+         /* complex version of: p[j+1] -= p[j]*root[i] */
+         p_re[j+1] -= p_re[j]*re[i] - p_im[j]*im[i];
+         p_im[j+1] -= p_re[j]*im[i] + p_im[j]*re[i];
+      }
+   }
+   for (i=0;i<ord+1;i++)
+      p[i]=p_re[i];
+}
+
+/*LPC polynomial "flatifier"*/
+void lpc_flat(float gamma, float *lpc_in, float *lpc_out, int order)
+{
+   int i;
+   float re[10], im[10], conv[10];
+   float alpha;
+   alpha = 1/(4-4*gamma);
+   poly_roots(lpc_in, re, im, conv, 10, 20, 7);
+   for (i=0;i<order;i++)
+   {
+      float fact,tmp;
+      float radius = sqrt(re[i]*re[i]+im[i]*im[i]);
+      tmp=1-radius;
+      if (tmp>2-2*gamma)
+         fact = tmp;
+      else
+         fact = alpha*tmp*tmp-gamma+1;
+      fact = (1-fact)/(radius+.001);
+      re[i]*=fact;
+      im[i]*=fact;
+   }
+   poly(re, im, lpc_out, order);
+}
+
 void syn_filt(float *x, float *a, float *y, int N, int ord)
 {
    int i,j;
index 4fffec6..af73b3d 100644 (file)
@@ -25,6 +25,11 @@ void print_vec(float *vec, int len, char *name);
 /* Apply bandwidth expansion on LPC coef */
 void bw_lpc(float gamma, float *lpc_in, float *lpc_out, int order);
 
+void poly(float *re, float *im, float *p, int ord);
+
+/*LPC polynomial "flatifier"*/
+void lpc_flat(float gamma, float *lpc_in, float *lpc_out, int order);
+
 /* Synthesis filter using the past of y[n] (negative indices) as memory */
 void syn_filt(float *x, float *a, float *y, int N, int ord);
 
index 8db10f4..baa7d83 100644 (file)
@@ -45,14 +45,21 @@ extern float hexc_10_32_table[];
 static pf_params pf_params_nb = {
    0.65,      /* formant enhancement numerator */
    0.7,      /* formant enhancement denominator */
+   0.2       /* pitch enhancement factor */
+};
+
+/* Post-filter parameters for narrowband */
+static pf_params pf_params_med = {
+   0.65,      /* formant enhancement numerator */
+   0.72,      /* formant enhancement denominator */
    0.4       /* pitch enhancement factor */
 };
 
 /* Post-filter parameters for low bit-rate narrowband */
 static pf_params pf_params_lbr = {
    0.65,      /* formant enhancement numerator */
-   0.72,      /* formant enhancement denominator */
-   0.4       /* pitch enhancement factor */
+   0.75,      /* formant enhancement denominator */
+   0.6       /* pitch enhancement factor */
 };
 
 /* Post-filter parameters for wideband */
@@ -213,7 +220,7 @@ static SpeexSubmode nb_submode4 = {
    split_cb_nogain_unquant,
    &split_cb_nb_med,
    nb_post_filter,
-   &pf_params_nb
+   &pf_params_med
 };
 
 static SpeexSubmode nb_submode5 = {
index 66170a0..3581edd 100644 (file)
@@ -75,10 +75,10 @@ float *stack)
    formant_den =   voiced_fact    * params->formant_enh_den 
                  + (1-voiced_fact)* params->formant_enh_num;
 
-   /*Short-term post-filter: A(z/g1)/A(z/.g1)*/
-   bw_lpc (formant_num, ak, awk, p);
+   /*Short-term post-filter using "flatified" versions of ak*/
+   lpc_flat (formant_num, ak, awk, p);
    residue_mem(new_exc, awk, tmp_exc, nsf, p, mem);
-   bw_lpc (formant_den, ak, awk, p);
+   lpc_flat (formant_den, ak, awk, p);
    syn_filt_mem(tmp_exc, awk, new_exc, nsf, p, mem2);
 
    /*Gain after enhancement*/
diff --git a/libspeex/roots.c b/libspeex/roots.c
new file mode 100644 (file)
index 0000000..ed683d1
--- /dev/null
@@ -0,0 +1,259 @@
+/* Copyright (C) 1981-1999 Ken Turkowski. <turk@computer.org>
+ *
+ * All rights reserved.
+ *
+ * Warranty Information
+ *  Even though I have reviewed this software, I make no warranty
+ *  or representation, either express or implied, with respect to this
+ *  software, its quality, accuracy, merchantability, or fitness for a
+ *  particular purpose.  As a result, this software is provided "as is,"
+ *  and you, its user, are assuming the entire risk as to its quality
+ *  and accuracy.
+ *
+ * This code may be used and freely distributed as long as it includes
+ * this copyright notice and the above warranty information.
+
+
+   Code slightly modified by Jean-Marc Valin (2002)
+
+   Speex License:
+
+   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
+
+ */
+
+#include <math.h>
+
+/*******************************************************************************
+ * FindPolynomialRoots
+ *
+ * The Bairstow and Newton correction formulae are used for a simultaneous
+ * linear and quadratic iterated synthetic division.  The coefficients of
+ * a polynomial of degree n are given as a[i] (i=0,i,..., n) where a[0] is
+ * the constant term.  The coefficients are scaled by dividing them by
+ * their geometric mean.  The Bairstow or Newton iteration method will
+ * nearly always converge to the number of figures carried, fig, either to
+ * root values or to their reciprocals.  If the simultaneous Newton and
+ * Bairstow iteration fails to converge on root values or their
+ * reciprocals in maxiter iterations, the convergence requirement will be
+ * successively reduced by one decimal figure.  This program anticipates
+ * and protects against loss of significance in the quadratic synthetic
+ * division.  (Refer to "On Programming the Numerical Solution of
+ * Polynomial Equations," by K. W. Ellenberger, Commun. ACM 3 (Dec. 1960),
+ * 644-647.)  The real and imaginary part of each root is stated as u[i]
+ * and v[i], respectively, together with the corresponding constant,
+ * conv[i], used in the convergence test.  This program has been used
+ * successfully for over a year on the Bendix G15-D (Intercard System) and
+ * has recently been coded for the IBM 709 (Fortran system).
+ *
+ * ACM algorithm #30 - Numerical Solution of the Polynomial Equation
+ * K. W. Ellenberger
+ * Missle Division, North American Aviation, Downey, California
+ * Converted to C, modified, optimized, and structured by
+ * Ken Turkowski
+ * CADLINC, Inc., Palo Alto, California
+ *******************************************************************************/
+
+#define MAXN 20
+
+void
+poly_roots(
+       const float             *a,                     /* Coefficients */
+       float                   *u,                     /* Real component of each root */
+       float                   *v,                     /* Imaginary component of each root */
+       float                   *conv,          /* Convergence constant associated with each root */
+       register long   n,                      /* Degree of polynomial (order-1) */
+       long                    maxiter,        /* Maximum number of iterations */
+       long                    fig                     /* The number of decimal figures to be computed */
+)
+{
+       int i;
+       register int j;
+       float h[MAXN + 3], b[MAXN + 3], c[MAXN + 3], d[MAXN + 3], e[MAXN + 3];
+       /* [-2 : n] */
+       float K, ps, qs, pt, qt, s, rev, r=0;
+       int t;
+       float p=0, q=0;
+
+       /* Zero elements with negative indices */
+       b[2 + -1] = b[2 + -2] =
+       c[2 + -1] = c[2 + -2] =
+       d[2 + -1] = d[2 + -2] =
+       e[2 + -1] = e[2 + -2] =
+       h[2 + -1] = h[2 + -2] = 0.0;
+
+       /* Copy polynomial coefficients to working storage */
+       for (j = 0; j <= n; j++)
+               h[2 + j] = *a++;                                                /* Note reversal of coefficients */
+
+       t = 1;
+       K = pow(10.0, (double)(fig));                           /* Relative accuracy */
+
+       for (; h[2 + n] == 0.0; n--) {                          /* Look for zero high-order coeff. */
+               *u++ = 0.0;
+               *v++ = 0.0;
+               *conv++ = K;
+       }
+
+INIT:
+       if (n == 0)
+               return;
+
+       ps = qs = pt = qt = s = 0.0;
+       rev = 1.0;
+       K = pow(10.0, (double)(fig));
+
+       if (n == 1) {
+               r = -h[2 + 1] / h[2 + 0];
+               goto LINEAR;
+       }
+
+       for (j = n; j >= 0; j--)                                        /* Find geometric mean of coeff's */
+               if (h[2 + j] != 0.0)
+                       s += log(fabs(h[2 + j]));
+       s = exp(s / (n + 1));
+
+       for (j = n; j >= 0; j--)                                        /* Normalize coeff's by mean */
+               h[2 + j] /= s;
+
+       if (fabs(h[2 + 1] / h[2 + 0]) < fabs(h[2 + n - 1] / h[2 + n])) {
+REVERSE:
+               t = -t;
+               for (j = (n - 1) / 2; j >= 0; j--) {
+                       s = h[2 + j];
+                       h[2 + j] = h[2 + n - j];
+                       h[2 + n - j] = s;
+               }
+       }
+       if (qs != 0.0) {
+               p = ps;
+               q = qs;
+       } else {
+               if (h[2 + n - 2] == 0.0) {
+                       q = 1.0;
+                       p = -2.0;
+               } else {
+                       q = h[2 + n] / h[2 + n - 2];
+                       p = (h[2 + n - 1] - q * h[2 + n - 3]) / h[2 + n - 2];
+               }
+               if (n == 2)
+                       goto QADRTIC;
+               r = 0.0;
+       }
+ITERATE:
+       for (i = maxiter; i > 0; i--) {
+
+               for (j = 0; j <= n; j++) {                              /* BAIRSTOW */
+                       b[2 + j] = h[2 + j] - p * b[2 + j - 1] - q * b[2 + j - 2];
+                       c[2 + j] = b[2 + j] - p * c[2 + j - 1] - q * c[2 + j - 2];
+               }
+               if ((h[2 + n - 1] != 0.0) && (b[2 + n - 1] != 0.0)) {
+                       if (fabs(h[2 + n - 1] / b[2 + n - 1]) >= K) {
+                               b[2 + n] = h[2 + n] - q * b[2 + n - 2];
+                       }
+                       if (b[2 + n] == 0.0)
+                               goto QADRTIC;
+                       if (K < fabs(h[2 + n] / b[2 + n]))
+                               goto QADRTIC;
+               }
+
+               for (j = 0; j <= n; j++) {                              /* NEWTON */
+                       d[2 + j] = h[2 + j] + r * d[2 + j - 1];/* Calculate polynomial at r */
+                       e[2 + j] = d[2 + j] + r * e[2 + j - 1];/* Calculate derivative at r */
+               }
+               if (d[2 + n] == 0.0)
+                       goto LINEAR;
+               if (K < fabs(h[2 + n] / d[2 + n]))
+                       goto LINEAR;
+
+               c[2 + n - 1] = -p * c[2 + n - 2] - q * c[2 + n - 3];
+               s = c[2 + n - 2] * c[2 + n - 2] - c[2 + n - 1] * c[2 + n - 3];
+               if (s == 0.0) {
+                       p -= 2.0;
+                       q *= (q + 1.0);
+               } else {
+                       p += (b[2 + n - 1] * c[2 + n - 2] - b[2 + n] * c[2 + n - 3]) / s;
+                       q += (-b[2 + n - 1] * c[2 + n - 1] + b[2 + n] * c[2 + n - 2]) / s;
+               }
+               if (e[2 + n - 1] == 0.0)
+                       r -= 1.0;                                                       /* Minimum step */
+               else
+                       r -= d[2 + n] / e[2 + n - 1];           /* Newton's iteration */
+       }
+       ps = pt;
+       qs = qt;
+       pt = p;
+       qt = q;
+       if (rev < 0.0)
+               K /= 10.0;
+       rev = -rev;
+       goto REVERSE;
+
+LINEAR:
+       if (t < 0)
+               r = 1.0 / r;
+       n--;
+       *u++ = r;
+       *v++ = 0.0;
+       *conv++ = K;
+
+       for (j = n; j >= 0; j--) {                                      /* Polynomial deflation by lin-nomial */
+               if ((d[2 + j] != 0.0) && (fabs(h[2 + j] / d[2 + j]) < K))
+                       h[2 + j] = d[2 + j];
+               else
+                       h[2 + j] = 0.0;
+       }
+
+       if (n == 0)
+               return;
+       goto ITERATE;
+
+QADRTIC:
+       if (t < 0) {
+               p /= q;
+               q = 1.0 / q;
+       }
+       n -= 2;
+
+       if (0.0 < (q - (p * p / 4.0))) {                        /* Two complex roots */
+               *(u + 1) = *u = -p / 2.0;
+               u += 2;
+               s = sqrt(q - (p * p / 4.0));
+               *v++ = s;
+               *v++ = -s;
+       } else {                                                                        /* Two real roots */
+               s = sqrt(((p * p / 4.0)) - q);
+               if (p < 0.0)
+                       *u++ = -p / 2.0 + s;
+               else
+                       *u++ = -p / 2.0 - s;
+               *u++ = q / u[-1];
+               *v++ = 0.0;
+               *v++ = 0.0;
+       }
+       *conv++ = K;
+       *conv++ = K;
+
+       for (j = n; j >= 0; j--) {                                      /* Polynomial deflation by quadratic */
+               if ((b[2 + j] != 0.0) && (fabs(h[2 + j] / b[2 + j]) < K))
+                       h[2 + j] = b[2 + j];
+               else
+                       h[2 + j] = 0.0;
+       }
+       goto INIT;
+}
+
+
+#undef MAXN
diff --git a/libspeex/roots.h b/libspeex/roots.h
new file mode 100644 (file)
index 0000000..5cc6616
--- /dev/null
@@ -0,0 +1,46 @@
+/* Copyright (C) 1981-1999 Ken Turkowski. <turk@computer.org>
+ *
+ * All rights reserved.
+ *
+ * Warranty Information
+ *  Even though I have reviewed this software, I make no warranty
+ *  or representation, either express or implied, with respect to this
+ *  software, its quality, accuracy, merchantability, or fitness for a
+ *  particular purpose.  As a result, this software is provided "as is,"
+ *  and you, its user, are assuming the entire risk as to its quality
+ *  and accuracy.
+ *
+ * This code may be used and freely distributed as long as it includes
+ * this copyright notice and the above warranty information.
+
+
+   Code slightly modified by Jean-Marc Valin (2002)
+
+   Speex License:
+
+   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
+
+ */
+
+void
+poly_roots(
+       const float             *a,                     /* Coefficients */
+       float                   *u,                     /* Real component of each root */
+       float                   *v,                     /* Imaginary component of each root */
+       float                   *conv,          /* Convergence constant associated with each root */
+       register long   n,                      /* Degree of polynomial (order-1) */
+       long                    maxiter,        /* Maximum number of iterations */
+       long                    fig                     /* The number of decimal figures to be computed */
+        );