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