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.17 2000/02/09 22:04:13 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
  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
 223 double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
 224   /* map the input curve to a bark-scale curve for encoding */
 225
 226   int mapped=l->ln;
 227   double *work=alloca(sizeof(double)*mapped);
 228   int i,j,last=0;
 229
 230   memset(work,0,sizeof(double)*mapped);
 231
 232   /* Only the decode side is behavior-specced; for now in the encoder,
 233      we select the maximum value of each band as representative (this
 234      helps make sure peaks don't go out of range.  In error terms,
 235      selecting min would make more sense, but the codebook is trained
 236      numerically, so we don't actually lose.  We'd still want to
 237      use the original curve for error and noise estimation */
 238
 239   for(i=0;i<l->n;i++){
 240     int bark=l->linearmap[i];
 241     if(work[bark]<curve[i])work[bark]=curve[i];
 242     if(bark>last+1){
 243       /* If the bark scale is climbing rapidly, some bins may end up
 244          going unused.  This isn't a waste actually; it keeps the
 245          scale resolution even so that the LPC generator has an easy
 246          time.  However, if we leave the bins empty we lose energy.
 247          So, fill 'em in.  The decoder does not do anything with  he
 248          unused bins, so we can fill them anyway we like to end up
 249          with a better spectral curve */
 250
 251       /* we'll always have a bin zero, so we don't need to guard init */
 252       long span=bark-last;
 253       for(j=1;j<span;j++){
 254         double del=(double)j/span;
 255         work[j+last]=work[bark]*del+work[last]*(1.-del);
 256       }
 257     }
 258     last=bark;
 259   }
 260   for(i=0;i<mapped;i++)work[i]*=l->barknorm[i];
 261
 262   return vorbis_lpc_from_spectrum(work,lpc,l);
 263 }
 264
 265
 266 /* One can do this the long way by generating the transfer function in
 267    the time domain and taking the forward FFT of the result.  The
 268    results from direct calculation are cleaner and faster.
 269
 270    This version does a linear curve generation and then later
 271    interpolates the log curve from the linear curve.  */
 272
 273 void _vlpc_de_helper(double *curve,double *lpc,double amp,
 274                             lpc_lookup *l){
 275   int i;
 276   memset(curve,0,sizeof(double)*l->ln*2);
 277   if(amp==0)return;
 278
 279   for(i=0;i<l->m;i++){
 280     curve[i*2+1]=lpc[i]/(4*amp);
 281     curve[i*2+2]=-lpc[i]/(4*amp);
 282   }
 283
 284   drft_backward(&l->fft,curve); /* reappropriated ;-) */
 285
 286   {
 287     int l2=l->ln*2;
 288     double unit=1./amp;
 289     curve[0]=(1./(curve[0]*2+unit));
 290     for(i=1;i<l->ln;i++){
 291       double real=(curve[i]+curve[l2-i]);
 292       double imag=(curve[i]-curve[l2-i]);
 293       curve[i]=(1./hypot(real+unit,imag));
 294     }
 295   }
 296 }
 297
 298 /* generate the whole freq response curve of an LPC IIR filter */
 299
 300 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
 301   double *lcurve=alloca(sizeof(double)*(l->ln*2));
 302   int i;
 303
 304   if(amp==0){
 305     memset(curve,0,sizeof(double)*l->n);
 306     return;
 307   }
 308   _vlpc_de_helper(lcurve,lpc,amp,l);
 309
 310 #ifdef ANALYSIS
 311   {
 312     static int frameno=0;
 313     int j;
 314     FILE *out;
 315     char buffer[80];
 316
 317     sprintf(buffer,"loglpc%d.m",frameno++);
 318     out=fopen(buffer,"w+");
 319     for(j=0;j<l->ln;j++)
 320       fprintf(out,"%g\n",lcurve[j]);
 321     fclose(out);
 322
 323 #endif
 324
 325   for(i=0;i<l->ln;i++)lcurve[i]/=l->barknorm[i];
 326   for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
 327
 328 #ifdef ANALYSIS
 329
 330     sprintf(buffer,"lpc%d.m",frameno-1);
 331     out=fopen(buffer,"w+");
 332     for(j=0;j<l->n;j++)
 333       fprintf(out,"%g\n",curve[j]);
 334     fclose(out);
 335   }
 336 #endif
 337 }
 338
 339 /* subtract or add an lpc filter to data.  Vorbis doesn't actually use this. */
 340
 341 void vorbis_lpc_residue(double *coeff,double *prime,int m,
 342                         double *data,long n){
 343
 344   /* in: coeff[0...m-1] LPC coefficients
 345          prime[0...m-1] initial values
 346          data[0...n-1] data samples
 347     out: data[0...n-1] residuals from LPC prediction */
 348
 349   long i,j;
 350   double *work=alloca(sizeof(double)*(m+n));
 351   double y;
 352
 353   if(!prime)
 354     for(i=0;i<m;i++)
 355       work[i]=0;
 356   else
 357     for(i=0;i<m;i++)
 358       work[i]=prime[i];
 359
 360   for(i=0;i<n;i++){
 361     y=0;
 362     for(j=0;j<m;j++)
 363       y-=work[i+j]*coeff[m-j-1];
 364
 365     work[i+m]=data[i];
 366     data[i]-=y;
 367   }
 368 }
 369
 370
 371 void vorbis_lpc_predict(double *coeff,double *prime,int m,
 372                      double *data,long n){
 373
 374   /* in: coeff[0...m-1] LPC coefficients
 375          prime[0...m-1] initial values (allocated size of n+m-1)
 376          data[0...n-1] residuals from LPC prediction
 377     out: data[0...n-1] data samples */
 378
 379   long i,j,o,p;
 380   double y;
 381   double *work=alloca(sizeof(double)*(m+n));
 382
 383   if(!prime)
 384     for(i=0;i<m;i++)
 385       work[i]=0.;
 386   else
 387     for(i=0;i<m;i++)
 388       work[i]=prime[i];
 389
 390   for(i=0;i<n;i++){
 391     y=data[i];
 392     o=i;
 393     p=m;
 394     for(j=0;j<m;j++)
 395       y-=work[o++]*coeff[--p];
 396
 397     data[i]=work[o]=y;
 398   }
 399 }
 400