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.14 2000/01/20 04:43:01 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 <stdio.h>
  49 #include <string.h>
  50 #include <math.h>
  51 #include "os.h"
  52 #include "smallft.h"
  53 #include "lpc.h"
  54 #include "scales.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-1)*C)=mapped-1
 168      floor(bark(rate)*C)=mapped */
 169
 170   scale=mapped/toBARK(rate);
 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(((double)rate)/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   for(i=0;i<mapped;i++)
 198     l->barknorm[i]=fromBARK((i+1)/scale)-fromBARK(i/scale);
 199
 200   /* we cheat decoding the LPC spectrum via FFTs */
 201
 202   drft_init(&l->fft,mapped*2);
 203
 204 }
 205
 206 void lpc_clear(lpc_lookup *l){
 207   if(l){
 208     if(l->barknorm)free(l->barknorm);
 209     if(l->linearmap)free(l->linearmap);
 210     drft_clear(&l->fft);
 211   }
 212 }
 213
 214
 215 /* less efficient than the decode side (written for clarity).  We're
 216    not bottlenecked here anyway */
 217
 218 double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
 219   /* map the input curve to a bark-scale curve for encoding */
 220
 221   int mapped=l->ln;
 222   double *work=alloca(sizeof(double)*mapped);
 223   int i,j,last=0;
 224
 225   memset(work,0,sizeof(double)*mapped);
 226
 227   /* Only the decode side is behavior-specced; for now in the encoder,
 228      we select the maximum value of each band as representative (this
 229      helps make sure peaks don't go out of range.  In error terms,
 230      selecting min would make more sense, but the codebook is trained
 231      numerically, so we don't actually lose.  We'd still want to
 232      use the original curve for error and noise estimation */
 233
 234   for(i=0;i<l->n;i++){
 235     int bark=l->linearmap[i];
 236     if(work[bark]<curve[i])work[bark]=curve[i];
 237     if(bark>last+1){
 238       /* If the bark scale is climbing rapidly, some bins may end up
 239          going unused.  This isn't a waste actually; it keeps the
 240          scale resolution even so that the LPC generator has an easy
 241          time.  However, if we leave the bins empty we lose energy.
 242          So, fill 'em in.  The decoder does not do anything witht he
 243          unused bins, so we can fill them anyway we like to end up
 244          with a better spectral curve */
 245
 246       /* we'll always have a bin zero, so we don't need to guard init */
 247       long span=bark-last;
 248       for(j=1;j<span;j++){
 249         double del=(double)j/span;
 250         work[j+last]=work[bark]*del+work[last]*(1.-del);
 251       }
 252     }
 253     last=bark;
 254   }
 255   for(i=0;i<mapped;i++)work[i]*=l->barknorm[i];
 256
 257 #ifdef ANALYSIS
 258   {
 259     int j;
 260     FILE *out;
 261     char buffer[80];
 262     static int frameno=0;
 263
 264     sprintf(buffer,"prelpc%d.m",frameno);
 265     out=fopen(buffer,"w+");
 266     for(j=0;j<l->n;j++)
 267       fprintf(out,"%g\n",curve[j]);
 268     fclose(out);
 269     sprintf(buffer,"preloglpc%d.m",frameno++);
 270     out=fopen(buffer,"w+");
 271     for(j=0;j<l->ln;j++)
 272       fprintf(out,"%g\n",work[j]);
 273     fclose(out);
 274   }
 275 #endif
 276
 277   return vorbis_lpc_from_spectrum(work,lpc,l);
 278 }
 279
 280
 281 /* One can do this the long way by generating the transfer function in
 282    the time domain and taking the forward FFT of the result.  The
 283    results from direct calculation are cleaner and faster.
 284
 285    This version does a linear curve generation and then later
 286    interpolates the log curve from the linear curve.  This could stand
 287    optimization; it could both be more precise as well as not compute
 288    quite a few unused values */
 289
 290 void _vlpc_de_helper(double *curve,double *lpc,double amp,
 291                             lpc_lookup *l){
 292   int i;
 293   memset(curve,0,sizeof(double)*l->ln*2);
 294   if(amp==0)return;
 295
 296   for(i=0;i<l->m;i++){
 297     curve[i*2+1]=lpc[i]/4/amp;
 298     curve[i*2+2]=-lpc[i]/4/amp;
 299   }
 300
 301   drft_backward(&l->fft,curve); /* reappropriated ;-) */
 302
 303   {
 304     int l2=l->ln*2;
 305     double unit=1./amp;
 306     curve[0]=(1./(curve[0]*2+unit));
 307     for(i=1;i<l->ln;i++){
 308       double real=(curve[i]+curve[l2-i]);
 309       double imag=(curve[i]-curve[l2-i]);
 310       curve[i]=(1./hypot(real+unit,imag));
 311     }
 312   }
 313 }
 314
 315
 316 /* generate the whole freq response curve of an LPC IIR filter */
 317
 318 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
 319   double *lcurve=alloca(sizeof(double)*(l->ln*2));
 320   int i;
 321
 322   _vlpc_de_helper(lcurve,lpc,amp,l);
 323
 324 #ifdef ANALYSIS
 325   {
 326     static int frameno=0;
 327     int j;
 328     FILE *out;
 329     char buffer[80];
 330
 331     sprintf(buffer,"loglpc%d.m",frameno++);
 332     out=fopen(buffer,"w+");
 333     for(j=0;j<l->ln;j++)
 334       fprintf(out,"%g\n",lcurve[j]);
 335     fclose(out);
 336
 337 #endif
 338
 339   if(amp==0)return;
 340
 341   for(i=0;i<l->ln;i++)lcurve[i]/=l->barknorm[i];
 342   for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
 343
 344 #ifdef ANALYSIS
 345
 346     sprintf(buffer,"lpc%d.m",frameno-1);
 347     out=fopen(buffer,"w+");
 348     for(j=0;j<l->n;j++)
 349       fprintf(out,"%g\n",curve[j]);
 350     fclose(out);
 351   }
 352 #endif
 353 }
 354
 355 void vorbis_lpc_apply(double *residue,double *lpc,double amp,lpc_lookup *l){
 356   double *lcurve=alloca(sizeof(double)*((l->ln+l->n)*2));
 357   int i;
 358
 359   if(amp==0){
 360     memset(residue,0,l->n*sizeof(double));
 361   }else{
 362
 363     _vlpc_de_helper(lcurve,lpc,amp,l);
 364
 365     for(i=0;i<l->ln;i++)lcurve[i]/=l->barknorm[i];
 366     for(i=0;i<l->n;i++)
 367       if(residue[i]!=0)
 368         residue[i]*=lcurve[l->linearmap[i]];
 369   }
 370 }
 371
 372 /* subtract or add an lpc filter to data.  Vorbis doesn't actually use this. */
 373
 374 void vorbis_lpc_residue(double *coeff,double *prime,int m,
 375                         double *data,long n){
 376
 377   /* in: coeff[0...m-1] LPC coefficients
 378          prime[0...m-1] initial values
 379          data[0...n-1] data samples
 380     out: data[0...n-1] residuals from LPC prediction */
 381
 382   long i,j;
 383   double *work=alloca(sizeof(double)*(m+n));
 384   double y;
 385
 386   if(!prime)
 387     for(i=0;i<m;i++)
 388       work[i]=0;
 389   else
 390     for(i=0;i<m;i++)
 391       work[i]=prime[i];
 392
 393   for(i=0;i<n;i++){
 394     y=0;
 395     for(j=0;j<m;j++)
 396       y-=work[i+j]*coeff[m-j-1];
 397
 398     work[i+m]=data[i];
 399     data[i]-=y;
 400   }
 401 }
 402
 403
 404 void vorbis_lpc_predict(double *coeff,double *prime,int m,
 405                      double *data,long n){
 406
 407   /* in: coeff[0...m-1] LPC coefficients
 408          prime[0...m-1] initial values (allocated size of n+m-1)
 409          data[0...n-1] residuals from LPC prediction
 410     out: data[0...n-1] data samples */
 411
 412   long i,j,o,p;
 413   double y;
 414   double *work=alloca(sizeof(double)*(m+n));
 415
 416   if(!prime)
 417     for(i=0;i<m;i++)
 418       work[i]=0.;
 419   else
 420     for(i=0;i<m;i++)
 421       work[i]=prime[i];
 422
 423   for(i=0;i<n;i++){
 424     y=data[i];
 425     o=i;
 426     p=m;
 427     for(j=0;j<m;j++)
 428       y-=work[o++]*coeff[--p];
 429
 430     data[i]=work[o]=y;
 431   }
 432 }
 433
 434