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: Aug 22 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 memcof(double *data, int n, int m, double *d){
  67   int k,j,i;
  68   double p=0.,wk1[n],wk2[n],wkm[m],xms;
  69
  70   memset(wk1,0,sizeof(wk1));
  71   memset(wk2,0,sizeof(wk2));
  72   memset(wkm,0,sizeof(wkm));
  73
  74   for (j=0;j<n;j++) p += data[j]*data[j];
  75   xms=p/n;
  76
  77   wk1[0]=data[0];
  78   wk2[n-2]=data[n-1];
  79
  80   for (j=2;j<=n-1;j++) {
  81     wk1[j-1]=data[j-1];
  82     wk2[j-2]=data[j-1];
  83   }
  84
  85   for (k=1;k<=m;k++) {
  86     double num=0.,denom=0.;
  87     for (j=1;j<=(n-k);j++) {
  88
  89       num += wk1[j-1]*wk2[j-1];
  90       denom += wk1[j-1]*wk1[j-1] + wk2[j-1]*wk2[j-1];
  91
  92     }
  93
  94     d[k-1]=2.0*num/denom;
  95     xms *= (1.0-d[k-1]*d[k-1]);
  96
  97     for (i=1;i<=(k-1);i++)
  98       d[i-1]=wkm[i-1]-d[k-1]*wkm[k-i-1];
  99
 100     if (k == m) return xms;
 101
 102     for (i=1;i<=k;i++) wkm[i-1]=d[i-1];
 103     for (j=1;j<=(n-k-1);j++) {
 104
 105       wk1[j-1] -= wkm[k-1]*wk2[j-1];
 106       wk2[j-1]=wk2[j]-wkm[k-1]*wk1[j];
 107
 108     }
 109   }
 110 }
 111
 112 static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
 113   double aut[m+1],work[n+n],error;
 114   drft_lookup dl;
 115   int i,j;
 116
 117   /* input is a real curve. make it complex-real */
 118   for(i=0;i<n;i++){
 119     work[i*2]=curve[i];
 120     work[i*2+1]=0;
 121   }
 122
 123   n*=2;
 124   drft_init(&dl,n);
 125   drft_backward(&dl,work);
 126   drft_clear(&dl);
 127
 128   /* The autocorrelation will not be circular.  Shift, else we lose
 129      most of the power in the edges. */
 130
 131   for(i=0,j=n/2;i<n/2;){
 132     double temp=work[i];
 133     work[i++]=work[j];
 134     work[j++]=temp;
 135   }
 136
 137   /* autocorrelation, p+1 lag coefficients */
 138
 139   j=m+1;
 140   while(j--){
 141     double d=0;
 142     for(i=j;i<n;i++)d+=work[i]*work[i-j];
 143     aut[j]=d;
 144   }
 145
 146   /* Generate lpc coefficients from autocorr values */
 147
 148   error=aut[0];
 149   if(error==0){
 150     memset(lpc,0,m*sizeof(double));
 151     return 0;
 152   }
 153
 154   for(i=0;i<m;i++){
 155     double r=-aut[i+1];
 156
 157     /* Sum up this iteration's reflection coefficient; note that in
 158        Vorbis we don't save it.  If anyone wants to recycle this code
 159        and needs reflection coefficients, save the results of 'r' from
 160        each iteration. */
 161
 162     for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
 163     r/=error;
 164
 165     /* Update LPC coefficients and total error */
 166
 167     lpc[i]=r;
 168     for(j=0;j<i/2;j++){
 169       double tmp=lpc[j];
 170       lpc[j]+=r*lpc[i-1-j];
 171       lpc[i-1-j]+=r*tmp;
 172     }
 173     if(i%2)lpc[j]+=lpc[j]*r;
 174
 175     error*=1.0-r*r;
 176   }
 177
 178   /* we need the error value to know how big an impulse to hit the
 179      filter with later */
 180
 181   return error;
 182 }
 183
 184 /* One can do this the long way by generating the transfer function in
 185    the time domain and taking the forward FFT of the result.  The
 186    results from direct calculation are cleaner and faster. If one
 187    looks at the below in the context of the calling function, there's
 188    lots of redundant trig, but at least it's clear */
 189
 190 static double vorbis_lpc_magnitude(double w,double *lpc, int m){
 191   int k;
 192   double real=1,imag=0;
 193   double wn=w;
 194   for(k=0;k<m;k++){
 195     real+=lpc[k]*cos(wn);
 196     imag+=lpc[k]*sin(wn);
 197     wn+=w;
 198   }
 199   return(1./sqrt(real*real+imag*imag));
 200 }
 201
 202 /* On top of this basic LPC infrastructure we introduce two modifications:
 203
 204    1) Filter generation is limited in the resolution of features it
 205    can represent (this is more obvious when the filter is looked at as
 206    a set of LSP coefficients).  Human perception of the audio spectrum
 207    is logarithmic not only in amplitude, but also frequency.  Because
 208    the high frequency features we'll need to encode will be broader
 209    than the low frequency features, filter generation will be
 210    dominated by higher frequencies (when most of the energy is in the
 211    lowest frequencies, and greatest perceived resolution is in the
 212    midrange).  To avoid this effect, Vorbis encodes the frequency
 213    spectrum with a biased log frequency scale. The intent is to
 214    roughly equalize the sizes of the octaves (see xlogmap.h)
 215
 216    2) When we change the frequency scale, we also change the
 217    (apparent) relative energies of the bands; that is, on a log scale
 218    covering 5 octaves, the highest octave goes from being represented
 219    in half the bins, to only 1/32 of the bins.  If the amplitudes
 220    remain the same, we have divided the energy represented by the
 221    highest octave by 16 (as far as Levinson-Durbin is concerned).
 222    This will seriously skew filter generation, which bases calculation
 223    on the mean square error with respect to energy.  Thus, Vorbis
 224    normalizes the amplitudes of the log spectrum frequencies to keep
 225    the relative octave energies correct. */
 226
 227 /* n == size of vector to be used for filter, m == order of filter,
 228    oct == octaves in normalized scale, encode_p == encode (1) or
 229    decode (0) */
 230
 231 double lpc_init(lpc_lookup *l,int n, int m, int oct, int encode_p){
 232   double bias=LOG_BIAS(n,oct);
 233   double scale=(float)n/(float)oct; /* where n==mapped */
 234   int i;
 235
 236   l->n=n;
 237   l->m=m;
 238   l->escale=malloc(n*sizeof(double));
 239   l->dscale=malloc(n*sizeof(double));
 240   l->norm=malloc(n*sizeof(double));
 241
 242   for(i=0;i<n;i++){
 243     /* how much 'real estate' in the log domain does the bin in the
 244        linear domain represent? */
 245     double logA=LOG_X(i-.5,bias);
 246     double logB=LOG_X(i+.5,bias);
 247     l->norm[i]=logB-logA;  /* this much */
 248   }
 249
 250   /* the scale is encode/decode specific for algebraic simplicity */
 251
 252   if(encode_p){
 253     /* encode */
 254
 255     for(i=0;i<n;i++)
 256       l->escale[i]=LINEAR_X(i/scale,bias);
 257
 258   }
 259   /* decode; encode may use this too */
 260
 261   {
 262     double w=1./oct*M_PI;
 263     for(i=0;i<n;i++)
 264       l->dscale[i]=LOG_X(i,bias)*w;
 265   }
 266 }
 267
 268 void lpc_clear(lpc_lookup *l){
 269   if(l){
 270     if(l->escale)free(l->escale);
 271     free(l->dscale);
 272     free(l->norm);
 273   }
 274 }
 275
 276
 277 /* less efficient than the decode side (written for clarity).  We're
 278    not bottlenecked here anyway */
 279
 280 double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
 281   /* map the input curve to a log curve for encoding */
 282
 283   /* for clarity, mapped and n are both represented although setting
 284      'em equal is a decent rule of thumb. */
 285
 286   int n=l->n;
 287   int m=l->m;
 288   int mapped=n;
 289   double work[mapped];
 290   int i;
 291
 292   /* fairly correct for low frequencies, naieve for high frequencies
 293      (suffers from undersampling) */
 294
 295   for(i=0;i<mapped;i++){
 296     double lin=l->escale[i];
 297     int a=floor(lin);
 298     int b=ceil(lin);
 299     double del=lin-floor(lin);
 300
 301     work[i]=(curve[a]/l->norm[a]*(1.-del)+
 302              curve[b]/l->norm[b]*del);
 303
 304   }
 305
 306   memcpy(curve,work,sizeof(work));
 307
 308   return vorbis_gen_lpc(work,mapped,lpc,l->m);
 309 }
 310
 311 /* generate the whole freq response curve on an LPC IIR filter */
 312
 313 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
 314   int i;
 315   for(i=0;i<l->n;i++)
 316     curve[i]=vorbis_lpc_magnitude(l->dscale[i],lpc,l->m)*amp*l->norm[i];
 317 }
 318
 319 /* find frequency response of LPC filter only at nonsero residue
 320    points and apply the envelope to the residue */
 321
 322 void vorbis_lpc_apply(double *residue,double *lpc,double amp,lpc_lookup *l){
 323   int i;
 324   for(i=0;i<l->n;i++)
 325     if(residue[i])
 326       residue[i]*=vorbis_lpc_magnitude(l->dscale[i],lpc,l->m)*amp*l->norm[i];
 327 }
 328
 329