lib/lpc.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH    *
   6  * THIS SOURCE. 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.29 2001/02/02 03:51:56 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(float)*(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     float d=0;
  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   if(error==0){
  79     memset(lpc,0,m*sizeof(float));
  80     return 0;
  81   }
  82
  83   for(i=0;i<m;i++){
  84     float r=-aut[i+1];
  85
  86     /* Sum up this iteration's reflection coefficient; note that in
  87        Vorbis we don't save it.  If anyone wants to recycle this code
  88        and needs reflection coefficients, save the results of 'r' from
  89        each iteration. */
  90
  91     for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
  92     r/=error;
  93
  94     /* Update LPC coefficients and total error */
  95
  96     lpc[i]=r;
  97     for(j=0;j<i/2;j++){
  98       float tmp=lpc[j];
  99       lpc[j]+=r*lpc[i-1-j];
 100       lpc[i-1-j]+=r*tmp;
 101     }
 102     if(i%2)lpc[j]+=lpc[j]*r;
 103
 104     error*=1.0f-r*r;
 105   }
 106
 107   /* we need the error value to know how big an impulse to hit the
 108      filter with later */
 109
 110   return error;
 111 }
 112
 113 /* Input : n element envelope spectral curve
 114    Output: m lpc coefficients, excitation energy */
 115
 116 float vorbis_lpc_from_curve(float *curve,float *lpc,lpc_lookup *l){
 117   int n=l->ln;
 118   int m=l->m;
 119   float *work=alloca(sizeof(float)*(n+n));
 120   float fscale=.5f/n;
 121   int i,j;
 122
 123   /* input is a real curve. make it complex-real */
 124   /* This mixes phase, but the LPC generation doesn't care. */
 125   for(i=0;i<n;i++){
 126     work[i*2]=curve[i]*fscale;
 127     work[i*2+1]=0;
 128   }
 129   work[n*2-1]=curve[n-1]*fscale;
 130
 131   n*=2;
 132   drft_backward(&l->fft,work);
 133
 134   /* The autocorrelation will not be circular.  Shift, else we lose
 135      most of the power in the edges. */
 136
 137   for(i=0,j=n/2;i<n/2;){
 138     float temp=work[i];
 139     work[i++]=work[j];
 140     work[j++]=temp;
 141   }
 142
 143   /* we *could* shave speed here by skimping on the edges (thus
 144      speeding up the autocorrelation in vorbis_lpc_from_data) but we
 145      don't right now. */
 146
 147   return(vorbis_lpc_from_data(work,lpc,n,m));
 148 }
 149
 150 void lpc_init(lpc_lookup *l,long mapped, int m){
 151   memset(l,0,sizeof(lpc_lookup));
 152
 153   l->ln=mapped;
 154   l->m=m;
 155
 156   /* we cheat decoding the LPC spectrum via FFTs */
 157   drft_init(&l->fft,mapped*2);
 158
 159 }
 160
 161 void lpc_clear(lpc_lookup *l){
 162   if(l){
 163     drft_clear(&l->fft);
 164   }
 165 }
 166
 167 void vorbis_lpc_predict(float *coeff,float *prime,int m,
 168                      float *data,long n){
 169
 170   /* in: coeff[0...m-1] LPC coefficients
 171          prime[0...m-1] initial values (allocated size of n+m-1)
 172     out: data[0...n-1] data samples */
 173
 174   long i,j,o,p;
 175   float y;
 176   float *work=alloca(sizeof(float)*(m+n));
 177
 178   if(!prime)
 179     for(i=0;i<m;i++)
 180       work[i]=0.f;
 181   else
 182     for(i=0;i<m;i++)
 183       work[i]=prime[i];
 184
 185   for(i=0;i<n;i++){
 186     y=0;
 187     o=i;
 188     p=m;
 189     for(j=0;j<m;j++)
 190       y-=work[o++]*coeff[--p];
 191
 192     data[i]=work[o]=y;
 193   }
 194 }
 195
 196
 197
 198
 199