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.18 2000/02/23 09:24:28 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_spectrum(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
 131   n*=2;
 132   drft_backward(&l->fft,work);
 133
 134   /* The autocorrelation will not be circular.  Shift, else we lose
 135      most of the power in the edges. */
 136
 137   for(i=0,j=n/2;i<n/2;){
 138     double temp=work[i];
 139     work[i++]=work[j];
 140     work[j++]=temp;
 141   }
 142
 143   return(vorbis_lpc_from_data(work,lpc,n,m));
 144 }
 145
 146 /* initialize Bark scale and normalization lookups.  We could do this
 147    with static tables, but Vorbis allows a number of possible
 148    combinations, so it's best to do it computationally.
 149
 150    The below is authoritative in terms of defining scale mapping.
 151    Note that the scale depends on the sampling rate as well as the
 152    linear block and mapping sizes */
 153
 154 void lpc_init(lpc_lookup *l,int n, long mapped, long rate, int m){
 155   int i;
 156   double scale;
 157   memset(l,0,sizeof(lpc_lookup));
 158
 159   l->n=n;
 160   l->ln=mapped;
 161   l->m=m;
 162
 163   l->linearmap=malloc(n*sizeof(int));
 164   l->barknorm=malloc(mapped*sizeof(double));
 165
 166   /* we choose a scaling constant so that:
 167      floor(bark(rate/2-1)*C)=mapped-1
 168      floor(bark(rate/2)*C)=mapped */
 169
 170   scale=mapped/toBARK(rate/2.);
 171
 172   /* the mapping from a linear scale to a smaller bark scale is
 173      straightforward.  We do *not* make sure that the linear mapping
 174      does not skip bark-scale bins; the decoder simply skips them and
 175      the encoder may do what it wishes in filling them.  They're
 176      necessary in some mapping combinations to keep the scale spacing
 177      accurate */
 178   {
 179     int last=-1;
 180     for(i=0;i<n;i++){
 181       int val=floor( toBARK((rate/2.)/n*i) *scale); /* bark numbers
 182                                                             represent
 183                                                             band edges */
 184       if(val>=mapped)val=mapped; /* guard against the approximation */
 185       l->linearmap[i]=val;
 186       last=val;
 187     }
 188   }
 189
 190   /* 'Normalization' is just making sure that power isn't lost in the
 191      log scale by virtue of compressing the scale in higher
 192      frequencies.  We figure the weight of bands in proportion to
 193      their linear/bark width ratio below, again, authoritatively.  We
 194      use computed width (not the number of actual bins above) for
 195      smoothness in the scale; they should agree closely */
 196
 197   /* keep it 0. to 1., else the dynamic range starts spreading through
 198      all the squaring... */
 199
 200   for(i=0;i<mapped;i++)
 201     l->barknorm[i]=(fromBARK((i+1)/scale)-fromBARK(i/scale));
 202   for(i=0;i<mapped;i++)
 203     l->barknorm[i]/=l->barknorm[mapped-1];
 204
 205   /* we cheat decoding the LPC spectrum via FFTs */
 206
 207   drft_init(&l->fft,mapped*2);
 208
 209 }
 210
 211 void lpc_clear(lpc_lookup *l){
 212   if(l){
 213     if(l->barknorm)free(l->barknorm);
 214     if(l->linearmap)free(l->linearmap);
 215     drft_clear(&l->fft);
 216   }
 217 }
 218
 219
 220 /* less efficient than the decode side (written for clarity).  We're
 221    not bottlenecked here anyway */
 222 static int frameno=-1;
 223
 224 double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
 225   /* map the input curve to a bark-scale curve for encoding */
 226
 227   int mapped=l->ln;
 228   double *work=alloca(sizeof(double)*mapped);
 229   int i,j,last=0;
 230
 231   frameno++;
 232   _analysis_output("lpc_pre",frameno,curve,l->n);
 233
 234   memset(work,0,sizeof(double)*mapped);
 235
 236   /* Only the decode side is behavior-specced; for now in the encoder,
 237      we select the maximum value of each band as representative (this
 238      helps make sure peaks don't go out of range.  In error terms,
 239      selecting min would make more sense, but the codebook is trained
 240      numerically, so we don't actually lose.  We'd still want to
 241      use the original curve for error and noise estimation */
 242
 243   for(i=0;i<l->n;i++){
 244     int bark=l->linearmap[i];
 245     if(work[bark]<curve[i])work[bark]=curve[i];
 246     if(bark>last+1){
 247       /* If the bark scale is climbing rapidly, some bins may end up
 248          going unused.  This isn't a waste actually; it keeps the
 249          scale resolution even so that the LPC generator has an easy
 250          time.  However, if we leave the bins empty we lose energy.
 251          So, fill 'em in.  The decoder does not do anything with  he
 252          unused bins, so we can fill them anyway we like to end up
 253          with a better spectral curve */
 254
 255       /* we'll always have a bin zero, so we don't need to guard init */
 256       long span=bark-last;
 257       for(j=1;j<span;j++){
 258         double del=(double)j/span;
 259         work[j+last]=work[bark]*del+work[last]*(1.-del);
 260       }
 261     }
 262     last=bark;
 263   }
 264   _analysis_output("lpc_prelog",frameno,work,l->ln);
 265   for(i=0;i<mapped;i++)work[i]*=l->barknorm[i];
 266   _analysis_output("lpc_prelognorm",frameno,work,l->ln);
 267
 268   return vorbis_lpc_from_spectrum(work,lpc,l);
 269 }
 270
 271
 272 /* One can do this the long way by generating the transfer function in
 273    the time domain and taking the forward FFT of the result.  The
 274    results from direct calculation are cleaner and faster.
 275
 276    This version does a linear curve generation and then later
 277    interpolates the log curve from the linear curve.  */
 278
 279 void _vlpc_de_helper(double *curve,double *lpc,double amp,
 280                             lpc_lookup *l){
 281   int i;
 282   memset(curve,0,sizeof(double)*l->ln*2);
 283   if(amp==0)return;
 284
 285   for(i=0;i<l->m;i++){
 286     curve[i*2+1]=lpc[i]/(4*amp);
 287     curve[i*2+2]=-lpc[i]/(4*amp);
 288   }
 289
 290   drft_backward(&l->fft,curve); /* reappropriated ;-) */
 291
 292   {
 293     int l2=l->ln*2;
 294     double unit=1./amp;
 295     curve[0]=(1./(curve[0]*2+unit));
 296     for(i=1;i<l->ln;i++){
 297       double real=(curve[i]+curve[l2-i]);
 298       double imag=(curve[i]-curve[l2-i]);
 299       curve[i]=(1./hypot(real+unit,imag));
 300     }
 301   }
 302 }
 303
 304 /* generate the whole freq response curve of an LPC IIR filter */
 305
 306 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
 307   double *lcurve=alloca(sizeof(double)*(l->ln*2));
 308   int i;
 309
 310   if(amp==0){
 311     memset(curve,0,sizeof(double)*l->n);
 312     return;
 313   }
 314   _vlpc_de_helper(lcurve,lpc,amp,l);
 315   _analysis_output("lpc_lognorm",frameno,lcurve,l->ln);
 316
 317   for(i=0;i<l->ln;i++)lcurve[i]/=l->barknorm[i];
 318   _analysis_output("lpc_log",frameno,lcurve,l->ln);
 319   for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
 320   _analysis_output("lpc",frameno,curve,l->n);
 321
 322 }
 323
 324 /* subtract or add an lpc filter to data.  Vorbis doesn't actually use this. */
 325
 326 void vorbis_lpc_residue(double *coeff,double *prime,int m,
 327                         double *data,long n){
 328
 329   /* in: coeff[0...m-1] LPC coefficients
 330          prime[0...m-1] initial values
 331          data[0...n-1] data samples
 332     out: data[0...n-1] residuals from LPC prediction */
 333
 334   long i,j;
 335   double *work=alloca(sizeof(double)*(m+n));
 336   double y;
 337
 338   if(!prime)
 339     for(i=0;i<m;i++)
 340       work[i]=0;
 341   else
 342     for(i=0;i<m;i++)
 343       work[i]=prime[i];
 344
 345   for(i=0;i<n;i++){
 346     y=0;
 347     for(j=0;j<m;j++)
 348       y-=work[i+j]*coeff[m-j-1];
 349
 350     work[i+m]=data[i];
 351     data[i]-=y;
 352   }
 353 }
 354
 355
 356 void vorbis_lpc_predict(double *coeff,double *prime,int m,
 357                      double *data,long n){
 358
 359   /* in: coeff[0...m-1] LPC coefficients
 360          prime[0...m-1] initial values (allocated size of n+m-1)
 361          data[0...n-1] residuals from LPC prediction
 362     out: data[0...n-1] data samples */
 363
 364   long i,j,o,p;
 365   double y;
 366   double *work=alloca(sizeof(double)*(m+n));
 367
 368   if(!prime)
 369     for(i=0;i<m;i++)
 370       work[i]=0.;
 371   else
 372     for(i=0;i<m;i++)
 373       work[i]=prime[i];
 374
 375   for(i=0;i<n;i++){
 376     y=data[i];
 377     o=i;
 378     p=m;
 379     for(j=0;j<m;j++)
 380       y-=work[o++]*coeff[--p];
 381
 382     data[i]=work[o]=y;
 383   }
 384 }
 385