1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/libs/vorbis/lpc.c	Sat Feb 01 19:58:19 2014 +0200
     1.3 @@ -0,0 +1,160 @@
     1.4 +/********************************************************************
     1.5 + *                                                                  *
     1.6 + * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
     1.7 + * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
     1.8 + * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
     1.9 + * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
    1.10 + *                                                                  *
    1.11 + * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
    1.12 + * by the Xiph.Org Foundation http://www.xiph.org/                  *
    1.13 + *                                                                  *
    1.14 + ********************************************************************
    1.15 +
    1.16 +  function: LPC low level routines
    1.17 +  last mod: $Id: lpc.c 16227 2009-07-08 06:58:46Z xiphmont $
    1.18 +
    1.19 + ********************************************************************/
    1.20 +
    1.21 +/* Some of these routines (autocorrelator, LPC coefficient estimator)
    1.22 +   are derived from code written by Jutta Degener and Carsten Bormann;
    1.23 +   thus we include their copyright below.  The entirety of this file
    1.24 +   is freely redistributable on the condition that both of these
    1.25 +   copyright notices are preserved without modification.  */
    1.26 +
    1.27 +/* Preserved Copyright: *********************************************/
    1.28 +
    1.29 +/* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
    1.30 +Technische Universita"t Berlin
    1.31 +
    1.32 +Any use of this software is permitted provided that this notice is not
    1.33 +removed and that neither the authors nor the Technische Universita"t
    1.34 +Berlin are deemed to have made any representations as to the
    1.35 +suitability of this software for any purpose nor are held responsible
    1.36 +for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
    1.37 +THIS SOFTWARE.
    1.38 +
    1.39 +As a matter of courtesy, the authors request to be informed about uses
    1.40 +this software has found, about bugs in this software, and about any
    1.41 +improvements that may be of general interest.
    1.42 +
    1.43 +Berlin, 28.11.1994
    1.44 +Jutta Degener
    1.45 +Carsten Bormann
    1.46 +
    1.47 +*********************************************************************/
    1.48 +
    1.49 +#include <stdlib.h>
    1.50 +#include <string.h>
    1.51 +#include <math.h>
    1.52 +#include "os.h"
    1.53 +#include "smallft.h"
    1.54 +#include "lpc.h"
    1.55 +#include "scales.h"
    1.56 +#include "misc.h"
    1.57 +
    1.58 +/* Autocorrelation LPC coeff generation algorithm invented by
    1.59 +   N. Levinson in 1947, modified by J. Durbin in 1959. */
    1.60 +
    1.61 +/* Input : n elements of time doamin data
    1.62 +   Output: m lpc coefficients, excitation energy */
    1.63 +
    1.64 +float vorbis_lpc_from_data(float *data,float *lpci,int n,int m){
    1.65 +  double *aut=alloca(sizeof(*aut)*(m+1));
    1.66 +  double *lpc=alloca(sizeof(*lpc)*(m));
    1.67 +  double error;
    1.68 +  double epsilon;
    1.69 +  int i,j;
    1.70 +
    1.71 +  /* autocorrelation, p+1 lag coefficients */
    1.72 +  j=m+1;
    1.73 +  while(j--){
    1.74 +    double d=0; /* double needed for accumulator depth */
    1.75 +    for(i=j;i<n;i++)d+=(double)data[i]*data[i-j];
    1.76 +    aut[j]=d;
    1.77 +  }
    1.78 +
    1.79 +  /* Generate lpc coefficients from autocorr values */
    1.80 +
    1.81 +  /* set our noise floor to about -100dB */
    1.82 +  error=aut[0] * (1. + 1e-10);
    1.83 +  epsilon=1e-9*aut[0]+1e-10;
    1.84 +
    1.85 +  for(i=0;i<m;i++){
    1.86 +    double r= -aut[i+1];
    1.87 +
    1.88 +    if(error<epsilon){
    1.89 +      memset(lpc+i,0,(m-i)*sizeof(*lpc));
    1.90 +      goto done;
    1.91 +    }
    1.92 +
    1.93 +    /* Sum up this iteration's reflection coefficient; note that in
    1.94 +       Vorbis we don't save it.  If anyone wants to recycle this code
    1.95 +       and needs reflection coefficients, save the results of 'r' from
    1.96 +       each iteration. */
    1.97 +
    1.98 +    for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
    1.99 +    r/=error;
   1.100 +
   1.101 +    /* Update LPC coefficients and total error */
   1.102 +
   1.103 +    lpc[i]=r;
   1.104 +    for(j=0;j<i/2;j++){
   1.105 +      double tmp=lpc[j];
   1.106 +
   1.107 +      lpc[j]+=r*lpc[i-1-j];
   1.108 +      lpc[i-1-j]+=r*tmp;
   1.109 +    }
   1.110 +    if(i&1)lpc[j]+=lpc[j]*r;
   1.111 +
   1.112 +    error*=1.-r*r;
   1.113 +
   1.114 +  }
   1.115 +
   1.116 + done:
   1.117 +
   1.118 +  /* slightly damp the filter */
   1.119 +  {
   1.120 +    double g = .99;
   1.121 +    double damp = g;
   1.122 +    for(j=0;j<m;j++){
   1.123 +      lpc[j]*=damp;
   1.124 +      damp*=g;
   1.125 +    }
   1.126 +  }
   1.127 +
   1.128 +  for(j=0;j<m;j++)lpci[j]=(float)lpc[j];
   1.129 +
   1.130 +  /* we need the error value to know how big an impulse to hit the
   1.131 +     filter with later */
   1.132 +
   1.133 +  return error;
   1.134 +}
   1.135 +
   1.136 +void vorbis_lpc_predict(float *coeff,float *prime,int m,
   1.137 +                     float *data,long n){
   1.138 +
   1.139 +  /* in: coeff[0...m-1] LPC coefficients
   1.140 +         prime[0...m-1] initial values (allocated size of n+m-1)
   1.141 +    out: data[0...n-1] data samples */
   1.142 +
   1.143 +  long i,j,o,p;
   1.144 +  float y;
   1.145 +  float *work=alloca(sizeof(*work)*(m+n));
   1.146 +
   1.147 +  if(!prime)
   1.148 +    for(i=0;i<m;i++)
   1.149 +      work[i]=0.f;
   1.150 +  else
   1.151 +    for(i=0;i<m;i++)
   1.152 +      work[i]=prime[i];
   1.153 +
   1.154 +  for(i=0;i<n;i++){
   1.155 +    y=0;
   1.156 +    o=i;
   1.157 +    p=m;
   1.158 +    for(j=0;j<m;j++)
   1.159 +      y-=work[o++]*coeff[--p];
   1.160 +
   1.161 +    data[i]=work[o]=y;
   1.162 +  }
   1.163 +}