Oops... now the LSPs are there. Also, lpcSize now represents the order
authorjmvalin <jmvalin@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Thu, 14 Feb 2002 01:57:22 +0000 (01:57 +0000)
committerjmvalin <jmvalin@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Thu, 14 Feb 2002 01:57:22 +0000 (01:57 +0000)
instead of the filter length (including a[0]=1). Cleaner that way and
more like what everybody else is doing.

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

libspeex/lpc.c
libspeex/lsp.c [new file with mode: 0644]
libspeex/lsp.h [new file with mode: 0644]
libspeex/speex.c

index 21b3b0b..bff5f4b 100644 (file)
@@ -40,9 +40,6 @@ wld(
 {
    int i, j;  float r, error = ac[0];
 
-   lpc[0]=1; /*for compatibility*/
-   lpc++; /*for compatibility*/
-
    if (ac[0] == 0) {
       for (i = 0; i < p; i++) ref[i] = 0; return 0; }
 
@@ -82,7 +79,6 @@ void autocorr(
               int lag, int   n)
 {
    float d; int i;
-   lag++;   /*Compensate for the old routine*/
    while (lag--) {
       for (i = lag, d = 0; i < n; i++) d += x[i] * x[i-lag];
       ac[lag] = d;
diff --git a/libspeex/lsp.c b/libspeex/lsp.c
new file mode 100644 (file)
index 0000000..a4b884b
--- /dev/null
@@ -0,0 +1,291 @@
+/*---------------------------------------------------------------------------*\
+
+       FILE........: AKSLSPD.C
+       TYPE........: Turbo C
+       COMPANY.....: Voicetronix
+       AUTHOR......: David Rowe
+       DATE CREATED: 24/2/93
+
+
+    This file contains functions for LPC to LSP conversion and
+    LSP to LPC conversion. Note that the LSP coefficients are not in
+    radians format but in the x domain of the unit circle.
+
+\*---------------------------------------------------------------------------*/
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "lsp.h"
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: cheb_poly_eva()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 24/2/93
+
+    This function evalutes a series of chebyshev polynomials
+
+\*---------------------------------------------------------------------------*/
+
+
+
+float cheb_poly_eva(float *coef,float x,int m)
+/*  float coef[]       coefficients of the polynomial to be evaluated  */
+/*  float x            the point where polynomial is to be evaluated   */
+/*  int m              order of the polynomial                         */
+
+
+{
+    int i;
+    float *T,*t,*u,*v,sum;
+
+    /* Allocate memory for chebyshev series formulation */
+
+    if((T = (float *)malloc((m/2+1)*sizeof(float))) == NULL){
+       printf("not enough memory to allocate buffer\n");
+       exit(1);
+    }
+
+    /* Initialise pointers */
+
+    t = T;                             /* T[i-2]                       */
+    *t++ = 1.0;
+    u = t--;                           /* T[i-1]                       */
+    *u++ = x;
+    v = u--;                           /* T[i]                         */
+
+    /* Evaluate chebyshev series formulation using iterative approach  */
+
+    for(i=2;i<=m/2;i++)
+       *v++ = (2*x)*(*u++) - *t++;     /* T[i] = 2*x*T[i-1] - T[i-2]   */
+
+    sum=0.0;                           /* initialise sum to zero       */
+    t = T;                             /* reset pointer                */
+
+    /* Evaluate polynomial and return value also free memory space */
+
+    for(i=0;i<=m/2;i++)
+       sum+=coef[(m/2)-i]**t++;
+
+    free(T);
+    return sum;
+}
+
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: lpc_to_lsp()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 24/2/93
+
+    This function converts LPC coefficients to LSP
+    coefficients.
+
+\*---------------------------------------------------------------------------*/
+
+
+int lpc_to_lsp (float *a,int lpcrdr,float *freq,int nb,float delta)
+/*  float *a                   lpc coefficients                        */
+/*  int lpcrdr                 order of LPC coefficients (10)          */
+/*  float *freq                LSP frequencies in the x domain         */
+/*  int nb                     number of sub-intervals (4)             */
+/*  float delta                        grid spacing interval (0.02)            */
+
+
+{
+
+    float psuml,psumr,psumm,temp_xr,xl,xr,xm;
+    float temp_psumr,temp_qsumr;
+    int i,j,m,flag,k;
+    float *Q;                  /* ptrs for memory allocation           */
+    float *P;
+    float *px;                 /* ptrs of respective P'(z) & Q'(z)     */
+    float *qx;
+    float *p;
+    float *q;
+    float *pt;                 /* ptr used for cheb_poly_eval()
+                               whether P' or Q'                        */
+    int roots=0;               /* DR 8/2/94: number of roots found     */
+    flag = 1;                  /*  program is searching for a root when,
+                               1 else has found one                    */
+    m = lpcrdr/2;              /* order of P'(z) & Q'(z) polynimials   */
+
+
+    /* Allocate memory space for polynomials */
+
+    if((Q = (float *) malloc((m+1)*sizeof(float))) == NULL){
+       printf("not enough memory to allocate buffer\n");
+       exit(1);
+    }
+
+    if((P = (float *) malloc((m+1)*sizeof(float))) == NULL){
+       printf("not enough memory to allocate buffer\n");
+       exit(1);
+    }
+
+    /* determine P'(z)'s and Q'(z)'s coefficients where
+      P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */
+
+    px = P;                      /* initilaise ptrs                    */
+    qx = Q;
+    p = px;
+    q = qx;
+    *px++ = 1.0;
+    *qx++ = 1.0;
+    for(i=1;i<=m;i++){
+       *px++ = a[i]+a[lpcrdr+1-i]-*p++;
+       *qx++ = a[i]-a[lpcrdr+1-i]+*q++;
+    }
+    px = P;
+    qx = Q;
+    for(i=0;i<m;i++){
+       *px = 2**px;
+       *qx = 2**qx;
+        px++;
+        qx++;
+    }
+    px = P;                    /* re-initialise ptrs                   */
+    qx = Q;
+
+    /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z).
+    Keep alternating between the two polynomials as each zero is found         */
+
+    xr = 0;                    /* initialise xr to zero                */
+    xl = 1.0;                  /* start at point xl = 1                */
+
+
+    for(j=0;j<lpcrdr;j++){
+       if(j%2)                 /* determines whether P' or Q' is eval. */
+           pt = qx;
+       else
+           pt = px;
+
+       psuml = cheb_poly_eva(pt,xl,lpcrdr);    /* evals poly. at xl    */
+       flag = 1;
+       while(flag && (xr >= -1.0)){
+           xr = xl - delta ;                   /* interval spacing     */
+           psumr = cheb_poly_eva(pt,xr,lpcrdr);/* poly(xl-delta_x)     */
+           temp_psumr = psumr;
+           temp_xr = xr;
+
+    /* if no sign change increment xr and re-evaluate poly(xr). Repeat til
+    sign change.
+    if a sign change has occurred the interval is bisected and then
+    checked again for a sign change which determines in which
+    interval the zero lies in.
+    If there is no sign change between poly(xm) and poly(xl) set interval
+    between xm and xr else set interval between xl and xr and repeat till
+    root is located within the specified limits                        */
+
+           if((psumr*psuml)<0.0){
+               roots++;
+
+               psumm=psuml;
+               for(k=0;k<=nb;k++){
+                   xm = (xl+xr)/2;             /* bisect the interval  */
+                   psumm=cheb_poly_eva(pt,xm,lpcrdr);
+                   if(psumm*psuml>0.){
+                       psuml=psumm;
+                       xl=xm;
+                   }
+                   else{
+                       psumr=psumm;
+                       xr=xm;
+                   }
+               }
+
+              /* once zero is found, reset initial interval to xr      */
+              freq[j] = (xm);
+              xl = xm;
+              flag = 0;                /* reset flag for next search   */
+           }
+           else{
+               psuml=temp_psumr;
+               xl=temp_xr;
+           }
+       }
+    }
+    free(P);                           /* free memory space            */
+    free(Q);
+
+    return(roots);
+}
+
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: lsp_to_lpc()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 24/2/93
+
+    lsp_to_lpc: This function converts LSP coefficients to LPC
+    coefficients.
+
+\*---------------------------------------------------------------------------*/
+
+
+void lsp_to_lpc(float *freq,float *ak,int lpcrdr)
+/*  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;
+    float xout1,xout2,xin1,xin2;
+    float *Wp;
+    float *pw,*n1,*n2,*n3,*n4;
+    int m = lpcrdr/2;
+
+    if((Wp = (float *) malloc((4*m+2)*sizeof(float))) == NULL){
+       printf("not enough memory to allocate buffer\n");
+       exit(1);
+    }
+    pw = Wp;
+
+    /* initialise contents of array */
+
+    for(i=0;i<=4*m+1;i++){             /* set contents of buffer to 0 */
+       *pw++ = 0.0;
+    }
+
+    /* Set pointers up */
+
+    pw = Wp;
+    xin1 = 1.0;
+    xin2 = 1.0;
+
+    /* 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(j=0;j<=lpcrdr;j++){
+       for(i=0;i<m;i++){
+           n1 = pw+(i*4);
+           n2 = n1 + 1;
+           n3 = n2 + 1;
+           n4 = n3 + 1;
+           xout1 = xin1 - 2*(freq[2*i]) * *n1 + *n2;
+           xout2 = xin2 - 2*(freq[2*i+1]) * *n3 + *n4;
+           *n2 = *n1;
+           *n4 = *n3;
+           *n1 = xin1;
+           *n3 = xin2;
+           xin1 = xout1;
+           xin2 = xout2;
+       }
+       xout1 = xin1 + *(n4+1);
+       xout2 = xin2 - *(n4+2);
+       ak[j] = (xout1 + xout2)*0.5;
+       *(n4+1) = xin1;
+       *(n4+2) = xin2;
+
+       xin1 = 0.0;
+       xin2 = 0.0;
+    }
+    free(Wp);
+}
diff --git a/libspeex/lsp.h b/libspeex/lsp.h
new file mode 100644 (file)
index 0000000..a800446
--- /dev/null
@@ -0,0 +1,17 @@
+/*---------------------------------------------------------------------------*\\r
+\r
+       FILE........: AK2LSPD.H\r
+       TYPE........: Turbo C header file\r
+       COMPANY.....: Voicetronix\r
+       AUTHOR......: James Whitehall\r
+       DATE CREATED: 21/11/95\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+#ifndef __AK2LSPD__\r
+#define __AK2LSPD__\r
+\r
+int lpc_to_lsp (float *a, int lpcrdr, float *freq, int nb, float delta);\r
+void lsp_to_lpc(float *freq, float *ak, int lpcrdr);\r
+\r
+#endif /* __AK2LSPD__ */\r
index 5c2edf4..0a4b8f9 100644 (file)
@@ -26,7 +26,7 @@ void encoder_init(EncState *st)
       st->window[i]=.5*(1-cos(2*M_PI*i/st->windowSize));
    st->buf2 = malloc(st->windowSize*sizeof(float));
    st->lpc = malloc(st->lpcSize*sizeof(float));
-   st->autocorr = malloc(st->lpcSize*sizeof(float));
+   st->autocorr = malloc((st->lpcSize+1)*sizeof(float));
    st->lsp = malloc(st->lpcSize*sizeof(float));
    st->rc = malloc(st->lpcSize*sizeof(float));
 }
@@ -56,15 +56,15 @@ void encode(EncState *st, float *in, int *outSize, void *bits)
    for (i=0;i<st->windowSize;i++)
       st->buf2[i] = st->inBuf[i] * st->window[i];
    
-   autocorr(st->buf2, st->autocorr, st->lpcSize-1, st->windowSize);
+   autocorr(st->buf2, st->autocorr, st->lpcSize+1, st->windowSize);
    st->autocorr[0] += 1;        /* prevents nan */
    st->autocorr[0] *= 1.0001;   /* 40 dB noise floor */
    /*Should do lag windowing here */
-   error = wld(st->lpc, st->autocorr, st->rc, st->lpcSize-1);
+   error = wld(st->lpc, st->autocorr, st->rc, st->lpcSize);
    printf ("prediction error = %f, R[0] = %f, gain = %f\n", error, st->autocorr[0], 
            st->autocorr[0]/error);
-   roots=lpc_to_lsp (st->lpc, st->lpcSize-1, st->lsp, 4, 0.02);
-   for (i=0;i<st->lpcSize-1;i++)
+   roots=lpc_to_lsp (st->lpc, st->lpcSize, st->lsp, 4, 0.02);
+   for (i=0;i<roots;i++)
       printf("%f ", st->lsp[i]);
    printf ("\nfound %d roots\n", roots);
 }