/********************************************************************
* *
* THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
- * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
- * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH *
- * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
* *
- * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2000 *
- * by Monty <monty@xiph.org> and the XIPHOPHORUS Company *
- * http://www.xiph.org/ *
+ * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 *
+ * by the Xiph.Org Foundation http://www.xiph.org/ *
* *
********************************************************************
function: LPC low level routines
- last mod: $Id: lpc.c,v 1.28 2000/12/21 21:04:39 xiphmont Exp $
+ last mod: $Id$
********************************************************************/
/* Input : n elements of time doamin data
Output: m lpc coefficients, excitation energy */
-float vorbis_lpc_from_data(float *data,float *lpc,int n,int m){
- float *aut=alloca(sizeof(float)*(m+1));
- float error;
+float vorbis_lpc_from_data(float *data,float *lpci,int n,int m){
+ double *aut=alloca(sizeof(*aut)*(m+1));
+ double *lpc=alloca(sizeof(*lpc)*(m));
+ double error;
+ double epsilon;
int i,j;
/* autocorrelation, p+1 lag coefficients */
-
j=m+1;
while(j--){
- float d=0;
- for(i=j;i<n;i++)d+=data[i]*data[i-j];
+ double d=0; /* double needed for accumulator depth */
+ for(i=j;i<n;i++)d+=(double)data[i]*data[i-j];
aut[j]=d;
}
-
+
/* Generate lpc coefficients from autocorr values */
- error=aut[0];
- if(error==0){
- memset(lpc,0,m*sizeof(float));
- return 0;
- }
-
+ /* set our noise floor to about -100dB */
+ error=aut[0] * (1. + 1e-10);
+ epsilon=1e-9*aut[0]+1e-10;
+
for(i=0;i<m;i++){
- float r=-aut[i+1];
+ double r= -aut[i+1];
+
+ if(error<epsilon){
+ memset(lpc+i,0,(m-i)*sizeof(*lpc));
+ goto done;
+ }
/* Sum up this iteration's reflection coefficient; note that in
Vorbis we don't save it. If anyone wants to recycle this code
each iteration. */
for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
- r/=error;
+ r/=error;
/* Update LPC coefficients and total error */
-
+
lpc[i]=r;
for(j=0;j<i/2;j++){
- float tmp=lpc[j];
+ double tmp=lpc[j];
+
lpc[j]+=r*lpc[i-1-j];
lpc[i-1-j]+=r*tmp;
}
- if(i%2)lpc[j]+=lpc[j]*r;
-
- error*=1.0f-r*r;
- }
-
- /* we need the error value to know how big an impulse to hit the
- filter with later */
-
- return error;
-}
+ if(i&1)lpc[j]+=lpc[j]*r;
-/* Input : n element envelope spectral curve
- Output: m lpc coefficients, excitation energy */
+ error*=1.-r*r;
-float vorbis_lpc_from_curve(float *curve,float *lpc,lpc_lookup *l){
- int n=l->ln;
- int m=l->m;
- float *work=alloca(sizeof(float)*(n+n));
- float fscale=.5f/n;
- int i,j;
-
- /* input is a real curve. make it complex-real */
- /* This mixes phase, but the LPC generation doesn't care. */
- for(i=0;i<n;i++){
- work[i*2]=curve[i]*fscale;
- work[i*2+1]=0;
}
- work[n*2-1]=curve[n-1]*fscale;
-
- n*=2;
- drft_backward(&l->fft,work);
-
- /* The autocorrelation will not be circular. Shift, else we lose
- most of the power in the edges. */
-
- for(i=0,j=n/2;i<n/2;){
- float temp=work[i];
- work[i++]=work[j];
- work[j++]=temp;
- }
-
- /* we *could* shave speed here by skimping on the edges (thus
- speeding up the autocorrelation in vorbis_lpc_from_data) but we
- don't right now. */
-
- return(vorbis_lpc_from_data(work,lpc,n,m));
-}
-void lpc_init(lpc_lookup *l,long mapped, int m){
- memset(l,0,sizeof(lpc_lookup));
+ done:
- l->ln=mapped;
- l->m=m;
+ /* slightly damp the filter */
+ {
+ double g = .99;
+ double damp = g;
+ for(j=0;j<m;j++){
+ lpc[j]*=damp;
+ damp*=g;
+ }
+ }
- /* we cheat decoding the LPC spectrum via FFTs */
- drft_init(&l->fft,mapped*2);
+ for(j=0;j<m;j++)lpci[j]=(float)lpc[j];
-}
+ /* we need the error value to know how big an impulse to hit the
+ filter with later */
-void lpc_clear(lpc_lookup *l){
- if(l){
- drft_clear(&l->fft);
- }
+ return error;
}
void vorbis_lpc_predict(float *coeff,float *prime,int m,
float *data,long n){
- /* in: coeff[0...m-1] LPC coefficients
+ /* in: coeff[0...m-1] LPC coefficients
prime[0...m-1] initial values (allocated size of n+m-1)
out: data[0...n-1] data samples */
long i,j,o,p;
float y;
- float *work=alloca(sizeof(float)*(m+n));
+ float *work=alloca(sizeof(*work)*(m+n));
if(!prime)
for(i=0;i<m;i++)
p=m;
for(j=0;j<m;j++)
y-=work[o++]*coeff[--p];
-
+
data[i]=work[o]=y;
}
}
-
-
-
-
-