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-1999             *
   9  * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company       *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14   function: LPC low level routines
  15   author: Monty <monty@xiph.org>
  16   modifications by: Monty
  17   last modification date: Oct 11 1999
  18
  19  ********************************************************************/
  20
  21 /* Some of these routines (autocorrelator, LPC coefficient estimator)
  22    are derived from code written by Jutta Degener and Carsten Bormann;
  23    thus we include their copyright below.  The entirety of this file
  24    is freely redistributable on the condition that both of these
  25    copyright notices are preserved without modification.  */
  26
  27 /* Preserved Copyright: *********************************************/
  28
  29 /* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
  30 Technische Universita"t Berlin
  31
  32 Any use of this software is permitted provided that this notice is not
  33 removed and that neither the authors nor the Technische Universita"t
  34 Berlin are deemed to have made any representations as to the
  35 suitability of this software for any purpose nor are held responsible
  36 for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
  37 THIS SOFTWARE.
  38
  39 As a matter of courtesy, the authors request to be informed about uses
  40 this software has found, about bugs in this software, and about any
  41 improvements that may be of general interest.
  42
  43 Berlin, 28.11.1994
  44 Jutta Degener
  45 Carsten Bormann
  46
  47 *********************************************************************/
  48
  49 #include <stdlib.h>
  50 #include <stdio.h>
  51 #include <string.h>
  52 #include <math.h>
  53 #include "smallft.h"
  54 #include "lpc.h"
  55 #include "xlogmap.h"
  56
  57 /* This is pared down for Vorbis where we only use LPC to encode
  58    spectral envelope curves.  Thus we only are interested in
  59    generating the coefficients and recovering the curve from the
  60    coefficients.  Autocorrelation LPC coeff generation algorithm
  61    invented by N. Levinson in 1947, modified by J. Durbin in 1959. */
  62
  63 /* Input : n element envelope curve
  64    Output: m lpc coefficients, excitation energy */
  65
  66 double vorbis_gen_lpc(double *curve,double *lpc,lpc_lookup *l){
  67   int n=l->ln;
  68   int m=l->m;
  69   double aut[m+1],work[n+n],error;
  70   double fscale=.5/n;
  71   int i,j;
  72
  73   /* input is a real curve. make it complex-real */
  74   /* This mixes phase, but the LPC generation doesn't care. */
  75   for(i=0;i<n;i++){
  76     work[i*2]=curve[i]*fscale;
  77     work[i*2+1]=0;
  78   }
  79
  80   n*=2;
  81   drft_backward(&l->fft,work);
  82
  83   /* The autocorrelation will not be circular.  Shift, else we lose
  84      most of the power in the edges. */
  85
  86   for(i=0,j=n/2;i<n/2;){
  87     double temp=work[i];
  88     work[i++]=work[j];
  89     work[j++]=temp;
  90   }
  91
  92   /* autocorrelation, p+1 lag coefficients */
  93
  94   j=m+1;
  95   while(j--){
  96     double d=0;
  97     for(i=j;i<n;i++)d+=work[i]*work[i-j];
  98     aut[j]=d;
  99   }
 100
 101   /* Generate lpc coefficients from autocorr values */
 102
 103   error=aut[0];
 104   if(error==0){
 105     memset(lpc,0,m*sizeof(double));
 106     return 0;
 107   }
 108
 109   for(i=0;i<m;i++){
 110     double r=-aut[i+1];
 111
 112     /* Sum up this iteration's reflection coefficient; note that in
 113        Vorbis we don't save it.  If anyone wants to recycle this code
 114        and needs reflection coefficients, save the results of 'r' from
 115        each iteration. */
 116
 117     for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
 118     r/=error;
 119
 120     /* Update LPC coefficients and total error */
 121
 122     lpc[i]=r;
 123     for(j=0;j<i/2;j++){
 124       double tmp=lpc[j];
 125       lpc[j]+=r*lpc[i-1-j];
 126       lpc[i-1-j]+=r*tmp;
 127     }
 128     if(i%2)lpc[j]+=lpc[j]*r;
 129
 130     error*=1.0-r*r;
 131   }
 132
 133   /* we need the error value to know how big an impulse to hit the
 134      filter with later */
 135
 136   return error;
 137 }
 138
 139 /* On top of this basic LPC infrastructure we introduce two modifications:
 140
 141    1) Filter generation is limited in the resolution of features it
 142    can represent (this is more obvious when the filter is looked at as
 143    a set of LSP coefficients).  Human perception of the audio spectrum
 144    is logarithmic not only in amplitude, but also frequency.  Because
 145    the high frequency features we'll need to encode will be broader
 146    than the low frequency features, filter generation will be
 147    dominated by higher frequencies (when most of the energy is in the
 148    lowest frequencies, and greatest perceived resolution is in the
 149    midrange).  To avoid this effect, Vorbis encodes the frequency
 150    spectrum with a biased log frequency scale. The intent is to
 151    roughly equalize the sizes of the octaves (see xlogmap.h)
 152
 153    2) When we change the frequency scale, we also change the
 154    (apparent) relative energies of the bands; that is, on a log scale
 155    covering 5 octaves, the highest octave goes from being represented
 156    in half the bins, to only 1/32 of the bins.  If the amplitudes
 157    remain the same, we have divided the energy represented by the
 158    highest octave by 16 (as far as Levinson-Durbin is concerned).
 159    This will seriously skew filter generation, which bases calculation
 160    on the mean square error with respect to energy.  Thus, Vorbis
 161    normalizes the amplitudes of the log spectrum frequencies to keep
 162    the relative octave energies correct. */
 163
 164 /* n == size of vector to be used for filter, m == order of filter,
 165    oct == octaves in normalized scale, encode_p == encode (1) or
 166    decode (0) */
 167
 168 void lpc_init(lpc_lookup *l,int n, int mapped, int m, int oct, int encode_p){
 169   double bias=LOG_BIAS(n,oct);
 170   double scale=(float)mapped/(float)oct; /* where n==mapped */
 171   int i;
 172
 173   memset(l,0,sizeof(lpc_lookup));
 174
 175   l->n=n;
 176   l->ln=mapped;
 177   l->m=m;
 178   l->iscale=malloc(n*sizeof(int));
 179   l->norm=malloc(n*sizeof(double));
 180
 181   for(i=0;i<n;i++){
 182     /* how much 'real estate' in the log domain does the bin in the
 183        linear domain represent? */
 184     double logA=LOG_X(i,bias);
 185     double logB=LOG_X(i+1.,bias);
 186     l->norm[i]=logB-logA;  /* this much */
 187   }
 188
 189   /* the scale is encode/decode specific for algebraic simplicity */
 190
 191   if(encode_p){
 192     /* encode */
 193     l->bscale=malloc(n*sizeof(int));
 194     l->escale=malloc(n*sizeof(double));
 195
 196     for(i=0;i<n;i++){
 197       l->escale[i]=LINEAR_X(i/scale,bias);
 198       l->bscale[i]=rint(LOG_X(i,bias)*scale);
 199     }
 200
 201   }
 202   /* decode; encode may use this too */
 203
 204   drft_init(&l->fft,mapped*2);
 205   for(i=0;i<n;i++){
 206     l->iscale[i]=rint(LOG_X(i,bias)/oct*mapped);
 207     if(l->iscale[i]>=l->ln)l->iscale[i]=l->ln-1;
 208   }
 209 }
 210
 211 void lpc_clear(lpc_lookup *l){
 212   if(l){
 213     if(l->bscale)free(l->bscale);
 214     if(l->escale)free(l->escale);
 215     drft_clear(&l->fft);
 216     free(l->iscale);
 217     free(l->norm);
 218   }
 219 }
 220
 221
 222 /* less efficient than the decode side (written for clarity).  We're
 223    not bottlenecked here anyway */
 224
 225 double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
 226   /* map the input curve to a log curve for encoding */
 227
 228   /* for clarity, mapped and n are both represented although setting
 229      'em equal is a decent rule of thumb. The below must be reworked
 230      slightly if mapped != n */
 231
 232   int mapped=l->ln;
 233   double work[mapped];
 234   int i;
 235
 236   /* fairly correct for low frequencies, naieve for high frequencies
 237      (suffers from undersampling) */
 238   for(i=0;i<mapped;i++){
 239     double lin=l->escale[i];
 240     int a=floor(lin);
 241     int b=ceil(lin);
 242     double del=lin-floor(lin);
 243
 244     work[i]=(curve[a]/l->norm[a]*(1.-del)+
 245              curve[b]/l->norm[b]*del);
 246
 247   }
 248
 249   return vorbis_gen_lpc(work,lpc,l);
 250 }
 251
 252
 253 /* One can do this the long way by generating the transfer function in
 254    the time domain and taking the forward FFT of the result.  The
 255    results from direct calculation are cleaner and faster. If one
 256    looks at the below in the context of the calling function, there's
 257    lots of redundant trig, but at least it's clear */
 258
 259 /* This version does a linear curve generation and then later
 260    interpolates the log curve from the linear curve.  This could stand
 261    optimization; it could both be more precise as well as not compute
 262    quite a few unused values */
 263
 264 void _vlpc_de_helper(double *curve,double *lpc,double amp,
 265                             lpc_lookup *l){
 266   int i;
 267   memset(curve,0,sizeof(double)*l->ln*2);
 268   if(amp==0)return;
 269
 270   for(i=0;i<l->m;i++){
 271     curve[i*2+1]=lpc[i]/4/amp;
 272     curve[i*2+2]=-lpc[i]/4/amp;
 273   }
 274
 275   drft_backward(&l->fft,curve); /* reappropriated ;-) */
 276
 277   {
 278     int l2=l->ln*2;
 279     double unit=1./amp;
 280     curve[0]=(1./(curve[0]*2+unit));
 281     for(i=1;i<l->ln;i++){
 282       double real=(curve[i]+curve[l2-i]);
 283       double imag=(curve[i]-curve[l2-i]);
 284       curve[i]=(1./hypot(real+unit,imag));
 285     }
 286   }
 287 }
 288
 289
 290 /* generate the whole freq response curve on an LPC IIR filter */
 291
 292 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
 293   double lcurve[l->ln*2];
 294   int i;
 295
 296   _vlpc_de_helper(lcurve,lpc,amp,l);
 297   if(amp==0)return;
 298
 299   for(i=0;i<l->n;i++)
 300     curve[i]=lcurve[l->iscale[i]]*l->norm[i];
 301 }
 302
 303 void vorbis_lpc_apply(double *residue,double *lpc,double amp,lpc_lookup *l){
 304   double lcurve[l->ln*2];
 305   int i;
 306
 307   if(amp==0){
 308     memset(residue,0,l->n*sizeof(double));
 309   }else{
 310
 311     _vlpc_de_helper(lcurve,lpc,amp,l);
 312
 313     for(i=0;i<l->n;i++)
 314       residue[i]*=lcurve[l->iscale[i]]*l->norm[i];
 315   }
 316 }
 317
 318