[xiph-commits] r11522 - trunk/speex/libspeex
jm at svn.xiph.org
jm at svn.xiph.org
Mon Jun 5 04:23:37 PDT 2006
Author: jm
Date: 2006-06-05 04:23:33 -0700 (Mon, 05 Jun 2006)
New Revision: 11522
Added:
trunk/speex/libspeex/lsp_bfin.h
Modified:
trunk/speex/libspeex/Makefile.am
trunk/speex/libspeex/lsp.c
Log:
Rewrite of the lsp_to_lpc() algorithm by David Rowe. Removes a bunch of
multiplications by zero. Also, Blackfin implementation of cheb_poly_eva().
Modified: trunk/speex/libspeex/Makefile.am
===================================================================
--- trunk/speex/libspeex/Makefile.am 2006-06-05 02:05:44 UTC (rev 11521)
+++ trunk/speex/libspeex/Makefile.am 2006-06-05 11:23:33 UTC (rev 11522)
@@ -24,7 +24,7 @@
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@
Modified: trunk/speex/libspeex/lsp.c
===================================================================
--- trunk/speex/libspeex/lsp.c 2006-06-05 02:05:44 UTC (rev 11521)
+++ trunk/speex/libspeex/lsp.c 2006-06-05 11:23:33 UTC (rev 11522)
@@ -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 @@
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 @@
#ifdef FIXED_POINT
-
-
#define FREQ_SCALE 16384
/*#define ANGLE2X(a) (32768*cos(((a)/8192.)))*/
@@ -72,6 +106,10 @@
/*#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 @@
/*---------------------------------------------------------------------------*\
- 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 @@
return sum;
}
+#endif
#else
@@ -159,10 +199,10 @@
/*---------------------------------------------------------------------------*\
- 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 @@
return(roots);
}
-
/*---------------------------------------------------------------------------*\
FUNCTION....: lsp_to_lpc()
@@ -355,8 +394,7 @@
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 @@
/* 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]);
+
+ xout1 = xout2 = 0;
+
+ /* now generate remaining rows */
+
+ for(i=1;i<m;i++) {
+
+ 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 */
+
+ 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);
}
- /* Set pointers up */
+ /* process last row to extra a{k} */
- pw = Wp;
- xin1 = 1048576;
- xin2 = 1048576;
+ for(j=1;j<=lpcrdr;j++) {
+ int shift = QIMP-13;
- /* reconstruct P(z) and Q(z) by cascading second order
- polynomials in form 1 - 2xz(-1) +z(-2), where x is the
- LSP coefficient */
+ /* 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 */
- 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;
+ if (a < -32767) a = 32767;
+ if (a > 32767) a = 32767;
+ ak[j-1] = (short)a;
+
+ }
- xin1 = 0;
- xin2 = 0;
- }
}
+
#else
void lsp_to_lpc(spx_lsp_t *freq,spx_coef_t *ak,int lpcrdr, char *stack)
Added: trunk/speex/libspeex/lsp_bfin.h
===================================================================
--- trunk/speex/libspeex/lsp_bfin.h 2006-06-05 02:05:44 UTC (rev 11521)
+++ trunk/speex/libspeex/lsp_bfin.h 2006-06-05 11:23:33 UTC (rev 11522)
@@ -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
+
+
+
More information about the commits
mailing list