/********************************************************************
* *
- * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
- * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
- * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
- * PLEASE READ THESE TERMS DISTRIBUTING. *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
+ * 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 OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999 *
- * by 1999 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 https://xiph.org/ *
* *
********************************************************************
function: LPC low level routines
- author: Monty <monty@xiph.org>
- modifications by: Monty
- last modification date: Aug 02 1999
********************************************************************/
/* Some of these routines (autocorrelator, LPC coefficient estimator)
- are directly derived from and/or modified from code written by
- Jutta Degener and Carsten Bormann; thus we include their copyright
- below. The entirety of this file is freely redistributable on the
- condition that both of these copyright notices are preserved
- without modification. */
+ are derived from code written by Jutta Degener and Carsten Bormann;
+ thus we include their copyright below. The entirety of this file
+ is freely redistributable on the condition that both of these
+ copyright notices are preserved without modification. */
/* Preserved Copyright: *********************************************/
#include <stdlib.h>
#include <string.h>
#include <math.h>
+#include "os.h"
#include "smallft.h"
#include "lpc.h"
+#include "scales.h"
+#include "misc.h"
-/* This is pared down for Vorbis where we only use LPC to encode
- spectral envelope curves. Thus we only are interested in
- generating the coefficients and recovering the curve from the
- coefficients. Autocorrelation LPC coeff generation algorithm
- invented by N. Levinson in 1947, modified by J. Durbin in 1959. */
+/* Autocorrelation LPC coeff generation algorithm invented by
+ N. Levinson in 1947, modified by J. Durbin in 1959. */
-/* Input : n element envelope curve
- Output: m lpc coefficients, excitation energy */
+/* Input : n elements of time doamin data
+ Output: m lpc coefficients, excitation energy */
-double vorbis_curve_to_lpc(double *curve,int n,double *lpc,int m){
- double aut[m+1],work[n+n],error;
- drft_lookup dl;
+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;
- /* input is a real curve. make it complex-real */
- for(i=0;i<n;i++){
- work[i*2]=curve[i];
- work[i*2+1]=0;
- }
-
- n*=2;
- drft_init(&dl,n);
- drft_backward(&dl,work);
- drft_clear(&dl);
-
- /* 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;){
- double temp=work[i];
- work[i++]=work[j];
- work[j++]=temp;
- }
-
/* autocorrelation, p+1 lag coefficients */
-
j=m+1;
while(j--){
- double d=0;
- for(i=j;i<n;i++)d+=work[i]*work[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(double));
- 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++){
- double 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++){
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.0-r*r;
+ if(i&1)lpc[j]+=lpc[j]*r;
+
+ error*=1.-r*r;
+
}
+ done:
+
+ /* slightly damp the filter */
+ {
+ double g = .99;
+ double damp = g;
+ for(j=0;j<m;j++){
+ lpc[j]*=damp;
+ damp*=g;
+ }
+ }
+
+ 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 */
return error;
}
-/* One can do this the long way by generating the transfer function in
- the time domain and taking the forward FFT of the result. The
- results from direct calculation are cleaner and faster, however */
-
-static double vorbis_lpc_magnitude(double w,double *lpc, int m){
- int k;
- double real=1,imag=0;
- for(k=0;k<m;k++){
- real+=lpc[k]*cos((k+1)*w);
- imag+=lpc[k]*sin((k+1)*w);
- }
- return(1./sqrt(real*real+imag*imag));
-}
-
-/* generate the whole freq response curve on an LPC IIR filter */
+void vorbis_lpc_predict(float *coeff,float *prime,int m,
+ float *data,long n){
-void vorbis_lpc_to_curve(double *curve,int n,double *lpc, double amp,int m){
- int i;
- double w=1./n*M_PI;
- for(i=0;i<n;i++)
- curve[i]=vorbis_lpc_magnitude(i*w,lpc,m)*amp;
-}
+ /* 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 */
-/* find frequency response of LPC filter only at nonsero residue
- points and apply the envelope to the residue */
+ long i,j,o,p;
+ float y;
+ float *work=alloca(sizeof(*work)*(m+n));
-void vorbis_lpc_apply(double *residue,int n,double *lpc, double amp,int m){
- int i;
- double w=1./n*M_PI;
- for(i=0;i<n;i++)
- if(residue[i])
- residue[i]*=vorbis_lpc_magnitude(i*w,lpc,m)*amp;
-}
+ if(!prime)
+ for(i=0;i<m;i++)
+ work[i]=0.f;
+ else
+ for(i=0;i<m;i++)
+ work[i]=prime[i];
+ for(i=0;i<n;i++){
+ y=0;
+ o=i;
+ p=m;
+ for(j=0;j<m;j++)
+ y-=work[o++]*coeff[--p];
+ data[i]=work[o]=y;
+ }
+}