diff options
author | jmvalin <jmvalin@0101bb08-14d6-0310-b084-bc0e0c8e3800> | 2002-02-14 04:57:22 +0300 |
---|---|---|
committer | jmvalin <jmvalin@0101bb08-14d6-0310-b084-bc0e0c8e3800> | 2002-02-14 04:57:22 +0300 |
commit | ca57d81f1698fa8d688c2053289ddc15b2b72b87 (patch) | |
tree | b295228eda8154e8ff7491905ea97452e2a3d517 /libspeex | |
parent | 782fadcab36cbbf886c0a73a7ac0199800962267 (diff) |
Oops... now the LSPs are there. Also, lpcSize now represents the order
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
Diffstat (limited to 'libspeex')
-rw-r--r-- | libspeex/lpc.c | 4 | ||||
-rw-r--r-- | libspeex/lsp.c | 291 | ||||
-rw-r--r-- | libspeex/lsp.h | 17 | ||||
-rw-r--r-- | libspeex/speex.c | 10 |
4 files changed, 313 insertions, 9 deletions
diff --git a/libspeex/lpc.c b/libspeex/lpc.c index 21b3b0b..bff5f4b 100644 --- a/libspeex/lpc.c +++ b/libspeex/lpc.c @@ -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 index 0000000..a4b884b --- /dev/null +++ b/libspeex/lsp.c @@ -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 index 0000000..a800446 --- /dev/null +++ b/libspeex/lsp.h @@ -0,0 +1,17 @@ +/*---------------------------------------------------------------------------*\
+
+ FILE........: AK2LSPD.H
+ TYPE........: Turbo C header file
+ COMPANY.....: Voicetronix
+ AUTHOR......: James Whitehall
+ DATE CREATED: 21/11/95
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef __AK2LSPD__
+#define __AK2LSPD__
+
+int lpc_to_lsp (float *a, int lpcrdr, float *freq, int nb, float delta);
+void lsp_to_lpc(float *freq, float *ak, int lpcrdr);
+
+#endif /* __AK2LSPD__ */
diff --git a/libspeex/speex.c b/libspeex/speex.c index 5c2edf4..0a4b8f9 100644 --- a/libspeex/speex.c +++ b/libspeex/speex.c @@ -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); } |