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=1./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-.5,bias);
 185     double logB=LOG_X(i+.5,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   {
 206     double w=1./oct*M_PI;
 207     for(i=0;i<n;i++){
 208       l->iscale[i]=rint(LOG_X(i,bias)/oct*mapped);
 209       if(l->iscale[i]>=l->ln)l->iscale[i]=l->ln-1;
 210     }
 211   }
 212 }
 213
 214 void lpc_clear(lpc_lookup *l){
 215   if(l){
 216     if(l->bscale)free(l->bscale);
 217     if(l->escale)free(l->escale);
 218     drft_clear(&l->fft);
 219     free(l->iscale);
 220     free(l->norm);
 221   }
 222 }
 223
 224
 225 /* less efficient than the decode side (written for clarity).  We're
 226    not bottlenecked here anyway */
 227
 228 double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
 229   /* map the input curve to a log curve for encoding */
 230
 231   /* for clarity, mapped and n are both represented although setting
 232      'em equal is a decent rule of thumb. The below must be reworked
 233      slightly if mapped != n */
 234
 235   int mapped=l->ln;
 236   double work[mapped];
 237   int i;
 238
 239   /* fairly correct for low frequencies, naieve for high frequencies
 240      (suffers from undersampling) */
 241   for(i=0;i<mapped;i++){
 242     double lin=l->escale[i];
 243     int a=floor(lin);
 244     int b=ceil(lin);
 245     double del=lin-floor(lin);
 246
 247     work[i]=(curve[a]/l->norm[a]*(1.-del)+
 248              curve[b]/l->norm[b]*del);
 249
 250   }
 251
 252   return vorbis_gen_lpc(work,lpc,l);
 253 }
 254
 255
 256 /* One can do this the long way by generating the transfer function in
 257    the time domain and taking the forward FFT of the result.  The
 258    results from direct calculation are cleaner and faster. If one
 259    looks at the below in the context of the calling function, there's
 260    lots of redundant trig, but at least it's clear */
 261
 262 /* This version does a linear curve generation and then later
 263    interpolates the log curve from the linear curve.  This could stand
 264    optimization; it could both be more precise as well as not compute
 265    quite a few unused values */
 266
 267 static void _vlpc_de_helper(double *curve,double *lpc,double amp,
 268                             lpc_lookup *l){
 269   int i;
 270   memset(curve,0,sizeof(double)*l->ln*2);
 271
 272   for(i=0;i<l->m;i++){
 273     curve[i*2+1]=lpc[i]/4/amp;
 274     curve[i*2+2]=-lpc[i]/4/amp;
 275   }
 276
 277   drft_backward(&l->fft,curve); /* reappropriated ;-) */
 278
 279   {
 280     int l2=l->ln*2;
 281     double unit=1./amp;
 282     curve[0]=(1./(curve[0]+unit));
 283     for(i=1;i<l->ln;i++){
 284       double real=(curve[i]+curve[l2-i]);
 285       double imag=(curve[i]-curve[l2-i]);
 286       curve[i]=(1./hypot(real+unit,imag));
 287     }
 288   }
 289 }
 290
 291
 292 /* generate the whole freq response curve on an LPC IIR filter */
 293
 294 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
 295   double lcurve[l->ln*2];
 296   int i;
 297
 298   _vlpc_de_helper(lcurve,lpc,amp,l);
 299
 300   for(i=0;i<l->n;i++)
 301     curve[i]=lcurve[l->iscale[i]]*l->norm[i];
 302 }
 303
 304 void vorbis_lpc_apply(double *residue,double *lpc,double amp,lpc_lookup *l){
 305   double lcurve[l->ln*2];
 306   int i;
 307
 308   _vlpc_de_helper(lcurve,lpc,amp,l);
 309
 310   for(i=0;i<l->n;i++)
 311     residue[i]*=lcurve[l->iscale[i]]*l->norm[i];
 312 }
 313