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