Rewrite of the lsp_to_lpc() algorithm by David Rowe. Removes a bunch of
authorjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Mon, 5 Jun 2006 11:23:33 +0000 (11:23 +0000)
committerjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Mon, 5 Jun 2006 11:23:33 +0000 (11:23 +0000)
multiplications by zero. Also, Blackfin implementation of cheb_poly_eva().

git-svn-id: http://svn.xiph.org/trunk/speex@11522 0101bb08-14d6-0310-b084-bc0e0c8e3800

libspeex/Makefile.am
libspeex/lsp.c
libspeex/lsp_bfin.h [new file with mode: 0644]

index 4d9a227..f178eba 100644 (file)
@@ -24,7 +24,7 @@ noinst_HEADERS = lsp.h        nb_celp.h       lpc.h   lpc_bfin.h      ltp.h   quant_lsp.h \
                                ltp_bfin.h      filters_sse.h   filters_arm4.h  filters_bfin.h  math_approx.h \
                                smallft.h       arch.h  fixed_arm4.h    fixed_arm5e.h   fixed_bfin.h    fixed_debug.h \
                                fixed_generic.h         cb_search_sse.h         cb_search_arm4.h        cb_search_bfin.h vorbis_psy.h \
-               fftwrap.h pseudofloat.h
+               fftwrap.h pseudofloat.h lsp_bfin.h
 
 
 libspeex_la_LDFLAGS = -no-undefined -version-info @SPEEX_LT_CURRENT@:@SPEEX_LT_REVISION@:@SPEEX_LT_AGE@
index 081f659..1023dcb 100644 (file)
@@ -1,7 +1,6 @@
 /*---------------------------------------------------------------------------*\
 Original copyright
-       FILE........: AKSLSPD.C
-       TYPE........: Turbo C
+       FILE........: lsp.c
        AUTHOR......: David Rowe
        DATE CREATED: 24/2/93
 
@@ -43,6 +42,43 @@ Heavily modified by Jean-Marc Valin (fixed-point, optimizations,
    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
+/*---------------------------------------------------------------------------*\
+
+  Introduction to Line Spectrum Pairs (LSPs)
+  ------------------------------------------
+
+  LSPs are used to encode the LPC filter coefficients {ak} for
+  transmission over the channel.  LSPs have several properties (like
+  less sensitivity to quantisation noise) that make them superior to
+  direct quantisation of {ak}.
+
+  A(z) is a polynomial of order lpcrdr with {ak} as the coefficients.
+
+  A(z) is transformed to P(z) and Q(z) (using a substitution and some
+  algebra), to obtain something like:
+
+    A(z) = 0.5[P(z)(z+z^-1) + Q(z)(z-z^-1)]  (1)
+
+  As you can imagine A(z) has complex zeros all over the z-plane. P(z)
+  and Q(z) have the very neat property of only having zeros _on_ the
+  unit circle.  So to find them we take a test point z=exp(jw) and
+  evaluate P (exp(jw)) and Q(exp(jw)) using a grid of points between 0
+  and pi.
+
+  The zeros (roots) of P(z) also happen to alternate, which is why we
+  swap coefficients as we find roots.  So the process of finding the
+  LSP frequencies is basically finding the roots of 5th order
+  polynomials.
+
+  The root so P(z) and Q(z) occur in symmetrical pairs at +/-w, hence
+  the name Line Spectrum Pairs (LSPs).
+
+  To convert back to ak we just evaluate (1), "clocking" an impulse
+  thru it lpcrdr times gives us the impulse response of A(z) which is
+  {ak}.
+
+\*---------------------------------------------------------------------------*/
+
 #ifdef HAVE_CONFIG_H
 #include "config.h"
 #endif
@@ -62,8 +98,6 @@ Heavily modified by Jean-Marc Valin (fixed-point, optimizations,
 
 #ifdef FIXED_POINT
 
-
-
 #define FREQ_SCALE 16384
 
 /*#define ANGLE2X(a) (32768*cos(((a)/8192.)))*/
@@ -72,6 +106,10 @@ Heavily modified by Jean-Marc Valin (fixed-point, optimizations,
 /*#define X2ANGLE(x) (acos(.00006103515625*(x))*LSP_SCALING)*/
 #define X2ANGLE(x) (spx_acos(x))
 
+#ifdef BFIN_ASM
+#include "lsp_bfin.h"
+#endif
+
 #else
 
 /*#define C1 0.99940307
@@ -87,17 +125,18 @@ Heavily modified by Jean-Marc Valin (fixed-point, optimizations,
 
 /*---------------------------------------------------------------------------*\
 
-       FUNCTION....: cheb_poly_eva()
+   FUNCTION....: cheb_poly_eva()
 
-       AUTHOR......: David Rowe
-       DATE CREATED: 24/2/93
+   AUTHOR......: David Rowe
+   DATE CREATED: 24/2/93
 
-    This function evaluates a series of Chebyshev polynomials
+   This function evaluates a series of Chebyshev polynomials
 
 \*---------------------------------------------------------------------------*/
 
 #ifdef FIXED_POINT
 
+#ifndef OVERRIDE_CHEB_POLY_EVA
 static inline spx_word32_t cheb_poly_eva(
   spx_word16_t *coef, /* P or Q coefs in Q13 format               */
   spx_word16_t     x, /* cos of freq (-1.0 to 1.0) in Q14 format  */
@@ -131,6 +170,7 @@ static inline spx_word32_t cheb_poly_eva(
     
     return sum;
 }
+#endif
 
 #else
 
@@ -159,10 +199,10 @@ static float cheb_poly_eva(spx_word32_t *coef, spx_word16_t x, int m, char *stac
 
 /*---------------------------------------------------------------------------*\
 
-       FUNCTION....: lpc_to_lsp()
+    FUNCTION....: lpc_to_lsp()
 
-       AUTHOR......: David Rowe
-       DATE CREATED: 24/2/93
+    AUTHOR......: David Rowe
+    DATE CREATED: 24/2/93
 
     This function converts LPC coefficients to LSP
     coefficients.
@@ -347,7 +387,6 @@ int lpc_to_lsp (spx_coef_t *a,int lpcrdr,spx_lsp_t *freq,int nb,spx_word16_t del
     return(roots);
 }
 
-
 /*---------------------------------------------------------------------------*\
 
        FUNCTION....: lsp_to_lpc()
@@ -355,8 +394,7 @@ int lpc_to_lsp (spx_coef_t *a,int lpcrdr,spx_lsp_t *freq,int nb,spx_word16_t del
        AUTHOR......: David Rowe
        DATE CREATED: 24/2/93
 
-    lsp_to_lpc: This function converts LSP coefficients to LPC
-    coefficients.
+        Converts LSP coefficients to LPC coefficients.
 
 \*---------------------------------------------------------------------------*/
 
@@ -366,77 +404,119 @@ void lsp_to_lpc(spx_lsp_t *freq,spx_coef_t *ak,int lpcrdr, char *stack)
 /*  float *freq        array of LSP frequencies in the x domain        */
 /*  float *ak          array of LPC coefficients                       */
 /*  int lpcrdr         order of LPC coefficients                       */
-
-
 {
     int i,j;
-    spx_word32_t xout1,xout2,xin1,xin2;
-    VARDECL(spx_word32_t *Wp);
-    spx_word32_t *pw,*n1,*n2,*n3,*n4=NULL;
+    spx_word32_t xout1,xout2,xin;
+    spx_word32_t mult, a;
     VARDECL(spx_word16_t *freqn);
+    VARDECL(spx_word32_t **xp);
+    VARDECL(spx_word32_t *xpmem);
+    VARDECL(spx_word32_t **xq);
+    VARDECL(spx_word32_t *xqmem);
     int m = lpcrdr>>1;
+
+    /* 
+    
+       Reconstruct P(z) and Q(z) by cascading second order polynomials
+       in form 1 - 2cos(w)z(-1) + z(-2), where w is the LSP frequency.
+       In the time domain this is:
+
+       y(n) = x(n) - 2cos(w)x(n-1) + x(n-2)
     
+       This is what the ALLOCS below are trying to do:
+
+         int xp[m+1][lpcrdr+1+2]; // P matrix in QIMP
+         int xq[m+1][lpcrdr+1+2]; // Q matrix in QIMP
+
+       These matrices store the output of each stage on each row.  The
+       final (m-th) row has the output of the final (m-th) cascaded
+       2nd order filter.  The first row is the impulse input to the
+       system (not written as it is known).
+
+       The version below takes advantage of the fact that a lot of the
+       outputs are zero or known, for example if we put an inpulse
+       into the first section the "clock" it 10 times only the first 3
+       outputs samples are non-zero (it's an FIR filter).
+    */
+
+    ALLOC(xp, (m+1), spx_word32_t*);
+    ALLOC(xpmem, (m+1)*(lpcrdr+1+2), spx_word32_t);
+
+    ALLOC(xq, (m+1), spx_word32_t*);
+    ALLOC(xqmem, (m+1)*(lpcrdr+1+2), spx_word32_t);
+    
+    for(i=0; i<=m; i++) {
+      xp[i] = xpmem + i*(lpcrdr+1+2);
+      xq[i] = xqmem + i*(lpcrdr+1+2);
+    }
+
+    /* work out 2cos terms in Q14 */
+
     ALLOC(freqn, lpcrdr, spx_word16_t);
-    for (i=0;i<lpcrdr;i++)
+    for (i=0;i<lpcrdr;i++) 
        freqn[i] = ANGLE2X(freq[i]);
 
-    ALLOC(Wp, 4*m+2, spx_word32_t);
-    pw = Wp;
+    #define QIMP  21   /* scaling for impulse */
 
+    xin = 1<<(QIMP-1); /* 0.5 in QIMP format */
+   
+    /* first col and last non-zero values of each row are trivial */
+    
+    for(i=0;i<=m;i++) {
+     xp[i][1] = 0;
+     xp[i][2] = xin;
+     xp[i][2+2*i] = xin;
+     xq[i][1] = 0;
+     xq[i][2] = xin;
+     xq[i][2+2*i] = xin;
+    }
 
-    /* initialise contents of array */
+    /* 2nd row (first output row) is trivial */
 
-    for(i=0;i<=4*m+1;i++){             /* set contents of buffer to 0 */
-       *pw++ = 0;
-    }
+    xp[1][3] = -MULT16_32_Q14(freqn[0],xp[0][2]);
+    xq[1][3] = -MULT16_32_Q14(freqn[1],xq[0][2]);
 
-    /* Set pointers up */
+    xout1 = xout2 = 0;
 
-    pw = Wp;
-    xin1 = 1048576;
-    xin2 = 1048576;
+    /* now generate remaining rows */
 
-    /* reconstruct P(z) and Q(z) by  cascading second order
-      polynomials in form 1 - 2xz(-1) +z(-2), where x is the
-      LSP coefficient */
+    for(i=1;i<m;i++) {
 
-    for(j=0;j<=lpcrdr;j++){
-       spx_word16_t *fr=freqn;
-       for(i=0;i<m;i++){
-           n1 = pw+(i<<2);
-           n2 = n1 + 1;
-           n3 = n2 + 1;
-           n4 = n3 + 1;
-           xout1 = ADD32(SUB32(xin1, MULT16_32_Q14(*fr,*n1)), *n2);
-            fr++;
-            xout2 = ADD32(SUB32(xin2, MULT16_32_Q14(*fr,*n3)), *n4);
-            fr++;
-           *n2 = *n1;
-           *n4 = *n3;
-           *n1 = xin1;
-           *n3 = xin2;
-           xin1 = xout1;
-           xin2 = xout2;
-       }
-       xout1 = xin1 + *(n4+1);
-       xout2 = xin2 - *(n4+2);
-        /* FIXME: perhaps apply bandwidth expansion in case of overflow? */
-       if (j>0)
-       {
-        if (xout1 + xout2>SHL32(EXTEND32(32766),8))
-           ak[j-1] = 32767;
-        else if (xout1 + xout2 < -SHL32(EXTEND32(32766),8))
-           ak[j-1] = -32767;
-        else
-           ak[j-1] = EXTRACT16(PSHR32(ADD32(xout1,xout2),8));
-       } else {/*speex_warning_int("ak[0] = ", EXTRACT16(PSHR32(ADD32(xout1,xout2),8)));*/}
-       *(n4+1) = xin1;
-       *(n4+2) = xin2;
+      for(j=1;j<2*(i+1)-1;j++) {
+       mult = MULT16_32_Q14(freqn[2*i],xp[i][j+1]);
+       xp[i+1][j+2] = ADD32(SUB32(xp[i][j+2], mult), xp[i][j]);
+       mult = MULT16_32_Q14(freqn[2*i+1],xq[i][j+1]);
+       xq[i+1][j+2] = ADD32(SUB32(xq[i][j+2], mult), xq[i][j]);
+      }
+
+      /* for last col xp[i][j+2] = xq[i][j+2] = 0 */
 
-       xin1 = 0;
-       xin2 = 0;
+      mult = MULT16_32_Q14(freqn[2*i],xp[i][j+1]);
+      xp[i+1][j+2] = SUB32(xp[i][j], mult);
+      mult = MULT16_32_Q14(freqn[2*i+1],xq[i][j+1]);
+      xq[i+1][j+2] = SUB32(xq[i][j], mult);
     }
+
+    /* process last row to extra a{k} */
+
+    for(j=1;j<=lpcrdr;j++) {
+      int shift = QIMP-13;
+
+      /* final filter sections */
+      a = PSHR32(xp[m][j+2] + xout1 + xq[m][j+2] - xout2, shift); 
+      xout1 = xp[m][j+2];
+      xout2 = xq[m][j+2];
+      
+      /* hard limit ak's to +/- 32767 */
+
+      if (a < -32767) a = 32767;
+      if (a > 32767) a = 32767;
+      ak[j-1] = (short)a;
+     
+    }
+
 }
+
 #else
 
 void lsp_to_lpc(spx_lsp_t *freq,spx_coef_t *ak,int lpcrdr, char *stack)
diff --git a/libspeex/lsp_bfin.h b/libspeex/lsp_bfin.h
new file mode 100644 (file)
index 0000000..20e5052
--- /dev/null
@@ -0,0 +1,89 @@
+/* Copyright (C) 2006 David Rowe */
+/**
+   @file lsp_bfin.h
+   @author David Rowe
+   @brief LSP routines optimised for the Blackfin
+*/
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+   
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+   
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+   
+   - Neither the name of the Xiph.org Foundation nor the names of its
+   contributors may be used to endorse or promote products derived from
+   this software without specific prior written permission.
+   
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#define OVERRIDE_CHEB_POLY_EVA
+#ifdef OVERRIDE_CHEB_POLY_EVA
+static inline spx_word32_t cheb_poly_eva(
+  spx_word16_t *coef, /* P or Q coefs in Q13 format               */
+  spx_word16_t     x, /* cos of freq (-1.0 to 1.0) in Q14 format  */
+  int              m, /* LPC order/2                              */
+  char         *stack
+)
+{
+    spx_word32_t sum;
+
+   __asm__ __volatile__
+     (
+      "P0 = %2;\n\t"           /* P0: coef[m], coef[m-1],..., coef[0] */
+      "R4 = 8192;\n\t"         /* R4: rounding constant               */
+      "R2 = %1;\n\t"           /* R2: x  */
+
+      "R5 = -16383;\n\t"
+      "R2 = MAX(R2,R5);\n\t"
+      "R5 = 16383;\n\t"
+      "R2 = MIN(R2,R5);\n\t"
+
+      "R3 = W[P0--] (X);\n\t"  /* R3: sum */
+      "R5 = W[P0--] (X);\n\t"
+      "R5 = R5.L * R2.L (IS);\n\t"
+      "R5 = R5 + R4;\n\t"
+      "R5 >>>= 14;\n\t"
+      "R3 = R3 + R5;\n\t" 
+      
+      "R0 = R2;\n\t"           /* R0: b0 */
+      "R1 = 16384;\n\t"        /* R1: b1 */
+      "LOOP cpe%= LC0 = %3;\n\t"
+      "LOOP_BEGIN cpe%=;\n\t"
+        "P1 = R0;\n\t" 
+        "R0 = R2.L * R0.L (IS) || R5 = W[P0--] (X);\n\t"
+        "R0 >>>= 13;\n\t"
+        "R0 = R0 - R1;\n\t"
+        "R1 = P1;\n\t"
+        "R5 = R5.L * R0.L (IS);\n\t"
+        "R5 = R5 + R4;\n\t"
+        "R5 >>>= 14;\n\t"
+        "R3 = R3 + R5;\n\t"
+      "LOOP_END cpe%=;\n\t"
+      "%0 = R3;\n\t"
+      : "=&d" (sum)
+      : "a" (x), "a" (&coef[m]), "a" (m-1)
+      : "R0", "R1", "R3", "R2", "R4", "R5", "P0", "P1"
+      );
+    return sum;
+}
+#endif
+
+
+