lib/lpc.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
   6  * PLEASE READ THESE TERMS DISTRIBUTING.                            *
   7  *                                                                  *
   8  * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
   9  * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14   function: LPC low level routines
  15   last mod: $Id: lpc.c,v 1.20 2000/05/08 20:49:49 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 /* Some of these routines (autocorrelator, LPC coefficient estimator)
  20    are derived from code written by Jutta Degener and Carsten Bormann;
  21    thus we include their copyright below.  The entirety of this file
  22    is freely redistributable on the condition that both of these
  23    copyright notices are preserved without modification.  */
  24
  25 /* Preserved Copyright: *********************************************/
  26
  27 /* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
  28 Technische Universita"t Berlin
  29
  30 Any use of this software is permitted provided that this notice is not
  31 removed and that neither the authors nor the Technische Universita"t
  32 Berlin are deemed to have made any representations as to the
  33 suitability of this software for any purpose nor are held responsible
  34 for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
  35 THIS SOFTWARE.
  36
  37 As a matter of courtesy, the authors request to be informed about uses
  38 this software has found, about bugs in this software, and about any
  39 improvements that may be of general interest.
  40
  41 Berlin, 28.11.1994
  42 Jutta Degener
  43 Carsten Bormann
  44
  45 *********************************************************************/
  46
  47 #include <stdlib.h>
  48 #include <string.h>
  49 #include <math.h>
  50 #include "os.h"
  51 #include "smallft.h"
  52 #include "lpc.h"
  53 #include "scales.h"
  54 #include "misc.h"
  55
  56 /* Autocorrelation LPC coeff generation algorithm invented by
  57    N. Levinson in 1947, modified by J. Durbin in 1959. */
  58
  59 /* Input : n elements of time doamin data
  60    Output: m lpc coefficients, excitation energy */
  61
  62 double vorbis_lpc_from_data(double *data,double *lpc,int n,int m){
  63   double *aut=alloca(sizeof(double)*(m+1));
  64   double error;
  65   int i,j;
  66
  67   /* autocorrelation, p+1 lag coefficients */
  68
  69   j=m+1;
  70   while(j--){
  71     double d=0;
  72     for(i=j;i<n;i++)d+=data[i]*data[i-j];
  73     aut[j]=d;
  74   }
  75
  76   /* Generate lpc coefficients from autocorr values */
  77
  78   error=aut[0];
  79   if(error==0){
  80     memset(lpc,0,m*sizeof(double));
  81     return 0;
  82   }
  83
  84   for(i=0;i<m;i++){
  85     double r=-aut[i+1];
  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       double 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.0-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 double vorbis_lpc_from_curve(double *curve,double *lpc,lpc_lookup *l){
 118   int n=l->ln;
 119   int m=l->m;
 120   double *work=alloca(sizeof(double)*(n+n));
 121   double fscale=.5/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
 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     double temp=work[i];
 139     work[i++]=work[j];
 140     work[j++]=temp;
 141   }
 142
 143   return(vorbis_lpc_from_data(work,lpc,n,m));
 144 }
 145
 146 void lpc_init(lpc_lookup *l,long mapped, int m){
 147   memset(l,0,sizeof(lpc_lookup));
 148
 149   l->ln=mapped;
 150   l->m=m;
 151
 152   /* we cheat decoding the LPC spectrum via FFTs */
 153   drft_init(&l->fft,mapped*2);
 154
 155 }
 156
 157 void lpc_clear(lpc_lookup *l){
 158   if(l){
 159     drft_clear(&l->fft);
 160   }
 161 }
 162
 163 /* One can do this the long way by generating the transfer function in
 164    the time domain and taking the forward FFT of the result.  The
 165    results from direct calculation are cleaner and faster.
 166
 167    This version does a linear curve generation and then later
 168    interpolates the log curve from the linear curve.  */
 169
 170 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,
 171                          lpc_lookup *l){
 172   int i;
 173   memset(curve,0,sizeof(double)*l->ln*2);
 174   if(amp==0)return;
 175
 176   for(i=0;i<l->m;i++){
 177     curve[i*2+1]=lpc[i]/(4*amp);
 178     curve[i*2+2]=-lpc[i]/(4*amp);
 179   }
 180
 181   drft_backward(&l->fft,curve); /* reappropriated ;-) */
 182
 183   {
 184     int l2=l->ln*2;
 185     double unit=1./amp;
 186     curve[0]=(1./(curve[0]*2+unit));
 187     for(i=1;i<l->ln;i++){
 188       double real=(curve[i]+curve[l2-i]);
 189       double imag=(curve[i]-curve[l2-i]);
 190
 191       double a = real + unit;
 192       curve[i] = 1.0 / FAST_HYPOT(a, imag);
 193     }
 194   }
 195 }
 196
 197 /* subtract or add an lpc filter to data.  Vorbis doesn't actually use this. */
 198
 199 void vorbis_lpc_residue(double *coeff,double *prime,int m,
 200                         double *data,long n){
 201
 202   /* in: coeff[0...m-1] LPC coefficients
 203          prime[0...m-1] initial values
 204          data[0...n-1] data samples
 205     out: data[0...n-1] residuals from LPC prediction */
 206
 207   long i,j;
 208   double *work=alloca(sizeof(double)*(m+n));
 209   double y;
 210
 211   if(!prime)
 212     for(i=0;i<m;i++)
 213       work[i]=0;
 214   else
 215     for(i=0;i<m;i++)
 216       work[i]=prime[i];
 217
 218   for(i=0;i<n;i++){
 219     y=0;
 220     for(j=0;j<m;j++)
 221       y-=work[i+j]*coeff[m-j-1];
 222
 223     work[i+m]=data[i];
 224     data[i]-=y;
 225   }
 226 }
 227
 228
 229 void vorbis_lpc_predict(double *coeff,double *prime,int m,
 230                      double *data,long n){
 231
 232   /* in: coeff[0...m-1] LPC coefficients
 233          prime[0...m-1] initial values (allocated size of n+m-1)
 234          data[0...n-1] residuals from LPC prediction
 235     out: data[0...n-1] data samples */
 236
 237   long i,j,o,p;
 238   double y;
 239   double *work=alloca(sizeof(double)*(m+n));
 240
 241   if(!prime)
 242     for(i=0;i<m;i++)
 243       work[i]=0.;
 244   else
 245     for(i=0;i<m;i++)
 246       work[i]=prime[i];
 247
 248   for(i=0;i<n;i++){
 249     y=data[i];
 250     o=i;
 251     p=m;
 252     for(j=0;j<m;j++)
 253       y-=work[o++]*coeff[--p];
 254
 255     data[i]=work[o]=y;
 256   }
 257 }
 258