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-2001             *
   9  * by the XIPHOPHORUS Company http://www.xiph.org/                  *
  10
  11  ********************************************************************
  12
  13   function: LPC low level routines
  14   last mod: $Id: lpc.c,v 1.33 2001/12/12 09:45:25 xiphmont Exp $
  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 #include <stdlib.h>
  47 #include <string.h>
  48 #include <math.h>
  49 #include "os.h"
  50 #include "smallft.h"
  51 #include "lpc.h"
  52 #include "scales.h"
  53 #include "misc.h"
  54
  55 /* Autocorrelation LPC coeff generation algorithm invented by
  56    N. Levinson in 1947, modified by J. Durbin in 1959. */
  57
  58 /* Input : n elements of time doamin data
  59    Output: m lpc coefficients, excitation energy */
  60
  61 float vorbis_lpc_from_data(float *data,float *lpc,int n,int m){
  62   float *aut=alloca(sizeof(*aut)*(m+1));
  63   float error;
  64   int i,j;
  65
  66   /* autocorrelation, p+1 lag coefficients */
  67
  68   j=m+1;
  69   while(j--){
  70     double d=0; /* double needed for accumulator depth */
  71     for(i=j;i<n;i++)d+=data[i]*data[i-j];
  72     aut[j]=d;
  73   }
  74
  75   /* Generate lpc coefficients from autocorr values */
  76
  77   error=aut[0];
  78
  79   for(i=0;i<m;i++){
  80     float r= -aut[i+1];
  81
  82     if(error==0){
  83       memset(lpc,0,m*sizeof(*lpc));
  84       return 0;
  85     }
  86
  87     /* Sum up this iteration's reflection coefficient; note that in
  88        Vorbis we don't save it.  If anyone wants to recycle this code
  89        and needs reflection coefficients, save the results of 'r' from
  90        each iteration. */
  91
  92     for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
  93     r/=error;
  94
  95     /* Update LPC coefficients and total error */
  96
  97     lpc[i]=r;
  98     for(j=0;j<i/2;j++){
  99       float tmp=lpc[j];
 100       lpc[j]+=r*lpc[i-1-j];
 101       lpc[i-1-j]+=r*tmp;
 102     }
 103     if(i%2)lpc[j]+=lpc[j]*r;
 104
 105     error*=1.f-r*r;
 106   }
 107
 108   /* we need the error value to know how big an impulse to hit the
 109      filter with later */
 110
 111   return error;
 112 }
 113
 114 /* Input : n element envelope spectral curve
 115    Output: m lpc coefficients, excitation energy */
 116
 117 float vorbis_lpc_from_curve(float *curve,float *lpc,lpc_lookup *l){
 118   int n=l->ln;
 119   int m=l->m;
 120   float *work=alloca(sizeof(*work)*(n+n));
 121   float fscale=.5f/n;
 122   int i,j;
 123
 124   /* input is a real curve. make it complex-real */
 125   /* This mixes phase, but the LPC generation doesn't care. */
 126   for(i=0;i<n;i++){
 127     work[i*2]=curve[i]*fscale;
 128     work[i*2+1]=0;
 129   }
 130   work[n*2-1]=curve[n-1]*fscale;
 131
 132   n*=2;
 133   drft_backward(&l->fft,work);
 134
 135   /* The autocorrelation will not be circular.  Shift, else we lose
 136      most of the power in the edges. */
 137
 138   for(i=0,j=n/2;i<n/2;){
 139     float temp=work[i];
 140     work[i++]=work[j];
 141     work[j++]=temp;
 142   }
 143
 144   /* we *could* shave speed here by skimping on the edges (thus
 145      speeding up the autocorrelation in vorbis_lpc_from_data) but we
 146      don't right now. */
 147
 148   return(vorbis_lpc_from_data(work,lpc,n,m));
 149 }
 150
 151 void lpc_init(lpc_lookup *l,long mapped, int m){
 152   memset(l,0,sizeof(*l));
 153
 154   l->ln=mapped;
 155   l->m=m;
 156
 157   /* we cheat decoding the LPC spectrum via FFTs */
 158   drft_init(&l->fft,mapped*2);
 159
 160 }
 161
 162 void lpc_clear(lpc_lookup *l){
 163   if(l){
 164     drft_clear(&l->fft);
 165   }
 166 }
 167
 168 void vorbis_lpc_predict(float *coeff,float *prime,int m,
 169                      float *data,long n){
 170
 171   /* in: coeff[0...m-1] LPC coefficients
 172          prime[0...m-1] initial values (allocated size of n+m-1)
 173     out: data[0...n-1] data samples */
 174
 175   long i,j,o,p;
 176   float y;
 177   float *work=alloca(sizeof(*work)*(m+n));
 178
 179   if(!prime)
 180     for(i=0;i<m;i++)
 181       work[i]=0.f;
 182   else
 183     for(i=0;i<m;i++)
 184       work[i]=prime[i];
 185
 186   for(i=0;i<n;i++){
 187     y=0;
 188     o=i;
 189     p=m;
 190     for(j=0;j<m;j++)
 191       y-=work[o++]*coeff[--p];
 192
 193     data[i]=work[o]=y;
 194   }
 195 }
 196
 197
 198
 199
 200