nuclear@1: /********************************************************************
nuclear@1:  *                                                                  *
nuclear@1:  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
nuclear@1:  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
nuclear@1:  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
nuclear@1:  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
nuclear@1:  *                                                                  *
nuclear@1:  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
nuclear@1:  * by the Xiph.Org Foundation http://www.xiph.org/                  *
nuclear@1:  *                                                                  *
nuclear@1:  ********************************************************************
nuclear@1: 
nuclear@1:   function: LPC low level routines
nuclear@1:   last mod: $Id: lpc.c 16227 2009-07-08 06:58:46Z xiphmont $
nuclear@1: 
nuclear@1:  ********************************************************************/
nuclear@1: 
nuclear@1: /* Some of these routines (autocorrelator, LPC coefficient estimator)
nuclear@1:    are derived from code written by Jutta Degener and Carsten Bormann;
nuclear@1:    thus we include their copyright below.  The entirety of this file
nuclear@1:    is freely redistributable on the condition that both of these
nuclear@1:    copyright notices are preserved without modification.  */
nuclear@1: 
nuclear@1: /* Preserved Copyright: *********************************************/
nuclear@1: 
nuclear@1: /* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
nuclear@1: Technische Universita"t Berlin
nuclear@1: 
nuclear@1: Any use of this software is permitted provided that this notice is not
nuclear@1: removed and that neither the authors nor the Technische Universita"t
nuclear@1: Berlin are deemed to have made any representations as to the
nuclear@1: suitability of this software for any purpose nor are held responsible
nuclear@1: for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
nuclear@1: THIS SOFTWARE.
nuclear@1: 
nuclear@1: As a matter of courtesy, the authors request to be informed about uses
nuclear@1: this software has found, about bugs in this software, and about any
nuclear@1: improvements that may be of general interest.
nuclear@1: 
nuclear@1: Berlin, 28.11.1994
nuclear@1: Jutta Degener
nuclear@1: Carsten Bormann
nuclear@1: 
nuclear@1: *********************************************************************/
nuclear@1: 
nuclear@1: #include <stdlib.h>
nuclear@1: #include <string.h>
nuclear@1: #include <math.h>
nuclear@1: #include "os.h"
nuclear@1: #include "smallft.h"
nuclear@1: #include "lpc.h"
nuclear@1: #include "scales.h"
nuclear@1: #include "misc.h"
nuclear@1: 
nuclear@1: /* Autocorrelation LPC coeff generation algorithm invented by
nuclear@1:    N. Levinson in 1947, modified by J. Durbin in 1959. */
nuclear@1: 
nuclear@1: /* Input : n elements of time doamin data
nuclear@1:    Output: m lpc coefficients, excitation energy */
nuclear@1: 
nuclear@1: float vorbis_lpc_from_data(float *data,float *lpci,int n,int m){
nuclear@1:   double *aut=alloca(sizeof(*aut)*(m+1));
nuclear@1:   double *lpc=alloca(sizeof(*lpc)*(m));
nuclear@1:   double error;
nuclear@1:   double epsilon;
nuclear@1:   int i,j;
nuclear@1: 
nuclear@1:   /* autocorrelation, p+1 lag coefficients */
nuclear@1:   j=m+1;
nuclear@1:   while(j--){
nuclear@1:     double d=0; /* double needed for accumulator depth */
nuclear@1:     for(i=j;i<n;i++)d+=(double)data[i]*data[i-j];
nuclear@1:     aut[j]=d;
nuclear@1:   }
nuclear@1: 
nuclear@1:   /* Generate lpc coefficients from autocorr values */
nuclear@1: 
nuclear@1:   /* set our noise floor to about -100dB */
nuclear@1:   error=aut[0] * (1. + 1e-10);
nuclear@1:   epsilon=1e-9*aut[0]+1e-10;
nuclear@1: 
nuclear@1:   for(i=0;i<m;i++){
nuclear@1:     double r= -aut[i+1];
nuclear@1: 
nuclear@1:     if(error<epsilon){
nuclear@1:       memset(lpc+i,0,(m-i)*sizeof(*lpc));
nuclear@1:       goto done;
nuclear@1:     }
nuclear@1: 
nuclear@1:     /* Sum up this iteration's reflection coefficient; note that in
nuclear@1:        Vorbis we don't save it.  If anyone wants to recycle this code
nuclear@1:        and needs reflection coefficients, save the results of 'r' from
nuclear@1:        each iteration. */
nuclear@1: 
nuclear@1:     for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
nuclear@1:     r/=error;
nuclear@1: 
nuclear@1:     /* Update LPC coefficients and total error */
nuclear@1: 
nuclear@1:     lpc[i]=r;
nuclear@1:     for(j=0;j<i/2;j++){
nuclear@1:       double tmp=lpc[j];
nuclear@1: 
nuclear@1:       lpc[j]+=r*lpc[i-1-j];
nuclear@1:       lpc[i-1-j]+=r*tmp;
nuclear@1:     }
nuclear@1:     if(i&1)lpc[j]+=lpc[j]*r;
nuclear@1: 
nuclear@1:     error*=1.-r*r;
nuclear@1: 
nuclear@1:   }
nuclear@1: 
nuclear@1:  done:
nuclear@1: 
nuclear@1:   /* slightly damp the filter */
nuclear@1:   {
nuclear@1:     double g = .99;
nuclear@1:     double damp = g;
nuclear@1:     for(j=0;j<m;j++){
nuclear@1:       lpc[j]*=damp;
nuclear@1:       damp*=g;
nuclear@1:     }
nuclear@1:   }
nuclear@1: 
nuclear@1:   for(j=0;j<m;j++)lpci[j]=(float)lpc[j];
nuclear@1: 
nuclear@1:   /* we need the error value to know how big an impulse to hit the
nuclear@1:      filter with later */
nuclear@1: 
nuclear@1:   return error;
nuclear@1: }
nuclear@1: 
nuclear@1: void vorbis_lpc_predict(float *coeff,float *prime,int m,
nuclear@1:                      float *data,long n){
nuclear@1: 
nuclear@1:   /* in: coeff[0...m-1] LPC coefficients
nuclear@1:          prime[0...m-1] initial values (allocated size of n+m-1)
nuclear@1:     out: data[0...n-1] data samples */
nuclear@1: 
nuclear@1:   long i,j,o,p;
nuclear@1:   float y;
nuclear@1:   float *work=alloca(sizeof(*work)*(m+n));
nuclear@1: 
nuclear@1:   if(!prime)
nuclear@1:     for(i=0;i<m;i++)
nuclear@1:       work[i]=0.f;
nuclear@1:   else
nuclear@1:     for(i=0;i<m;i++)
nuclear@1:       work[i]=prime[i];
nuclear@1: 
nuclear@1:   for(i=0;i<n;i++){
nuclear@1:     y=0;
nuclear@1:     o=i;
nuclear@1:     p=m;
nuclear@1:     for(j=0;j<m;j++)
nuclear@1:       y-=work[o++]*coeff[--p];
nuclear@1: 
nuclear@1:     data[i]=work[o]=y;
nuclear@1:   }
nuclear@1: }