Welcome to mirror list, hosted at ThFree Co, Russian Federation.

github.com/mumble-voip/speexdsp.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--libspeex/lpc.c4
-rw-r--r--libspeex/lsp.c291
-rw-r--r--libspeex/lsp.h17
-rw-r--r--libspeex/speex.c10
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);
}