lib/lpc.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
   5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
   6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
   9  * by the Xiph.Org Foundation https://xiph.org/                     *
  10  *                                                                  *
  11  ********************************************************************
  12
  13   function: LPC low level routines
  14
  15  ********************************************************************/
  16
  17 /* Some of these routines (autocorrelator, LPC coefficient estimator)
  18    are derived from code written by Jutta Degener and Carsten Bormann;
  19    thus we include their copyright below.  The entirety of this file
  20    is freely redistributable on the condition that both of these
  21    copyright notices are preserved without modification.  */
  22
  23 /* Preserved Copyright: *********************************************/
  24
  25 /* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
  26 Technische Universita"t Berlin
  27
  28 Any use of this software is permitted provided that this notice is not
  29 removed and that neither the authors nor the Technische Universita"t
  30 Berlin are deemed to have made any representations as to the
  31 suitability of this software for any purpose nor are held responsible
  32 for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
  33 THIS SOFTWARE.
  34
  35 As a matter of courtesy, the authors request to be informed about uses
  36 this software has found, about bugs in this software, and about any
  37 improvements that may be of general interest.
  38
  39 Berlin, 28.11.1994
  40 Jutta Degener
  41 Carsten Bormann
  42
  43 *********************************************************************/
  44
  45 #include <stdlib.h>
  46 #include <string.h>
  47 #include <math.h>
  48 #include "os.h"
  49 #include "smallft.h"
  50 #include "lpc.h"
  51 #include "scales.h"
  52 #include "misc.h"
  53
  54 /* Autocorrelation LPC coeff generation algorithm invented by
  55    N. Levinson in 1947, modified by J. Durbin in 1959. */
  56
  57 /* Input : n elements of time doamin data
  58    Output: m lpc coefficients, excitation energy */
  59
  60 float vorbis_lpc_from_data(float *data,float *lpci,int n,int m){
  61   double *aut=alloca(sizeof(*aut)*(m+1));
  62   double *lpc=alloca(sizeof(*lpc)*(m));
  63   double error;
  64   double epsilon;
  65   int i,j;
  66
  67   /* autocorrelation, p+1 lag coefficients */
  68   j=m+1;
  69   while(j--){
  70     double d=0; /* double needed for accumulator depth */
  71     for(i=j;i<n;i++)d+=(double)data[i]*data[i-j];
  72     aut[j]=d;
  73   }
  74
  75   /* Generate lpc coefficients from autocorr values */
  76
  77   /* set our noise floor to about -100dB */
  78   error=aut[0] * (1. + 1e-10);
  79   epsilon=1e-9*aut[0]+1e-10;
  80
  81   for(i=0;i<m;i++){
  82     double r= -aut[i+1];
  83
  84     if(error<epsilon){
  85       memset(lpc+i,0,(m-i)*sizeof(*lpc));
  86       goto done;
  87     }
  88
  89     /* Sum up this iteration's reflection coefficient; note that in
  90        Vorbis we don't save it.  If anyone wants to recycle this code
  91        and needs reflection coefficients, save the results of 'r' from
  92        each iteration. */
  93
  94     for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
  95     r/=error;
  96
  97     /* Update LPC coefficients and total error */
  98
  99     lpc[i]=r;
 100     for(j=0;j<i/2;j++){
 101       double tmp=lpc[j];
 102
 103       lpc[j]+=r*lpc[i-1-j];
 104       lpc[i-1-j]+=r*tmp;
 105     }
 106     if(i&1)lpc[j]+=lpc[j]*r;
 107
 108     error*=1.-r*r;
 109
 110   }
 111
 112  done:
 113
 114   /* slightly damp the filter */
 115   {
 116     double g = .99;
 117     double damp = g;
 118     for(j=0;j<m;j++){
 119       lpc[j]*=damp;
 120       damp*=g;
 121     }
 122   }
 123
 124   for(j=0;j<m;j++)lpci[j]=(float)lpc[j];
 125
 126   /* we need the error value to know how big an impulse to hit the
 127      filter with later */
 128
 129   return error;
 130 }
 131
 132 void vorbis_lpc_predict(float *coeff,float *prime,int m,
 133                      float *data,long n){
 134
 135   /* in: coeff[0...m-1] LPC coefficients
 136          prime[0...m-1] initial values (allocated size of n+m-1)
 137     out: data[0...n-1] data samples */
 138
 139   long i,j,o,p;
 140   float y;
 141   float *work=alloca(sizeof(*work)*(m+n));
 142
 143   if(!prime)
 144     for(i=0;i<m;i++)
 145       work[i]=0.f;
 146   else
 147     for(i=0;i<m;i++)
 148       work[i]=prime[i];
 149
 150   for(i=0;i<n;i++){
 151     y=0;
 152     o=i;
 153     p=m;
 154     for(j=0;j<m;j++)
 155       y-=work[o++]*coeff[--p];
 156
 157     data[i]=work[o]=y;
 158   }
 159 }