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.23 2000/08/15 09:09:43 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   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     double 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(lpc_lookup));
 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 /* One can do this the long way by generating the transfer function in
 169    the time domain and taking the forward FFT of the result.  The
 170    results from direct calculation are cleaner and faster.
 171
 172    This version does a linear curve generation and then later
 173    interpolates the log curve from the linear curve.  */
 174
 175 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,
 176                          lpc_lookup *l){
 177   int i;
 178   memset(curve,0,sizeof(double)*l->ln*2);
 179   if(amp==0)return;
 180
 181   for(i=0;i<l->m;i++){
 182     curve[i*2+1]=lpc[i]/(4*amp);
 183     curve[i*2+2]=-lpc[i]/(4*amp);
 184   }
 185
 186   drft_backward(&l->fft,curve); /* reappropriated ;-) */
 187
 188   {
 189     int l2=l->ln*2;
 190     double unit=1./amp;
 191     curve[0]=(1./(curve[0]*2+unit));
 192     for(i=1;i<l->ln;i++){
 193       double real=(curve[i]+curve[l2-i]);
 194       double imag=(curve[i]-curve[l2-i]);
 195
 196       double a = real + unit;
 197       curve[i] = 1.0 / FAST_HYPOT(a, imag);
 198     }
 199   }
 200 }
 201
 202 /* subtract or add an lpc filter to data. */
 203
 204 void vorbis_lpc_filter(double *coeff,double *prime,int m,
 205                         double *data,long n,double amp){
 206
 207   /* in: coeff[0...m-1] LPC coefficients
 208          prime[0...m-1] initial values
 209          data[0...n-1] data samples
 210     out: data[0...n-1] residuals from LPC prediction */
 211
 212   long i,j;
 213   double *work=alloca(sizeof(double)*(m+n));
 214   double y;
 215
 216   if(!prime)
 217     for(i=0;i<m;i++)
 218       work[i]=0;
 219   else
 220     for(i=0;i<m;i++)
 221       work[i]=prime[i];
 222
 223   for(i=0;i<n;i++){
 224     y=0;
 225     for(j=0;j<m;j++)
 226       y-=work[i+j]*coeff[m-j-1];
 227
 228     data[i]=work[i+m]=data[i]+y;
 229
 230   }
 231 }
 232
 233 void vorbis_lpc_residue(double *coeff,double *prime,int m,
 234                         double *data,long n){
 235
 236   /* in: coeff[0...m-1] LPC coefficients
 237          prime[0...m-1] initial values
 238          data[0...n-1] data samples
 239     out: data[0...n-1] residuals from LPC prediction */
 240
 241   long i,j;
 242   double *work=alloca(sizeof(double)*(m+n));
 243   double y;
 244
 245   if(!prime)
 246     for(i=0;i<m;i++)
 247       work[i]=0;
 248   else
 249     for(i=0;i<m;i++)
 250       work[i]=prime[i];
 251
 252   for(i=0;i<n;i++){
 253     y=0;
 254     for(j=0;j<m;j++)
 255       y-=work[i+j]*coeff[m-j-1];
 256
 257     work[i+m]=data[i];
 258     data[i]-=y;
 259   }
 260 }
 261
 262 void vorbis_lpc_predict(double *coeff,double *prime,int m,
 263                      double *data,long n){
 264
 265   /* in: coeff[0...m-1] LPC coefficients
 266          prime[0...m-1] initial values (allocated size of n+m-1)
 267          data[0...n-1] residuals from LPC prediction
 268     out: data[0...n-1] data samples */
 269
 270   long i,j,o,p;
 271   double y;
 272   double *work=alloca(sizeof(double)*(m+n));
 273
 274   if(!prime)
 275     for(i=0;i<m;i++)
 276       work[i]=0.;
 277   else
 278     for(i=0;i<m;i++)
 279       work[i]=prime[i];
 280
 281   for(i=0;i<n;i++){
 282     y=data[i];
 283     o=i;
 284     p=m;
 285     for(j=0;j<m;j++)
 286       y-=work[o++]*coeff[--p];
 287
 288     data[i]=work[o]=y;
 289   }
 290 }
 291