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.16 2000/02/06 13:39:42 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   if(amp==0){
 323     memset(curve,0,sizeof(double)*l->n);
 324     return;
 325   }
 326   _vlpc_de_helper(lcurve,lpc,amp,l);
 327
 328 #ifdef ANALYSIS
 329   {
 330     static int frameno=0;
 331     int j;
 332     FILE *out;
 333     char buffer[80];
 334
 335     sprintf(buffer,"loglpc%d.m",frameno++);
 336     out=fopen(buffer,"w+");
 337     for(j=0;j<l->ln;j++)
 338       fprintf(out,"%g\n",lcurve[j]);
 339     fclose(out);
 340
 341 #endif
 342
 343   for(i=0;i<l->ln;i++)lcurve[i]/=l->barknorm[i];
 344   for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
 345
 346 #ifdef ANALYSIS
 347
 348     sprintf(buffer,"lpc%d.m",frameno-1);
 349     out=fopen(buffer,"w+");
 350     for(j=0;j<l->n;j++)
 351       fprintf(out,"%g\n",curve[j]);
 352     fclose(out);
 353   }
 354 #endif
 355 }
 356
 357 /* subtract or add an lpc filter to data.  Vorbis doesn't actually use this. */
 358
 359 void vorbis_lpc_residue(double *coeff,double *prime,int m,
 360                         double *data,long n){
 361
 362   /* in: coeff[0...m-1] LPC coefficients
 363          prime[0...m-1] initial values
 364          data[0...n-1] data samples
 365     out: data[0...n-1] residuals from LPC prediction */
 366
 367   long i,j;
 368   double *work=alloca(sizeof(double)*(m+n));
 369   double y;
 370
 371   if(!prime)
 372     for(i=0;i<m;i++)
 373       work[i]=0;
 374   else
 375     for(i=0;i<m;i++)
 376       work[i]=prime[i];
 377
 378   for(i=0;i<n;i++){
 379     y=0;
 380     for(j=0;j<m;j++)
 381       y-=work[i+j]*coeff[m-j-1];
 382
 383     work[i+m]=data[i];
 384     data[i]-=y;
 385   }
 386 }
 387
 388
 389 void vorbis_lpc_predict(double *coeff,double *prime,int m,
 390                      double *data,long n){
 391
 392   /* in: coeff[0...m-1] LPC coefficients
 393          prime[0...m-1] initial values (allocated size of n+m-1)
 394          data[0...n-1] residuals from LPC prediction
 395     out: data[0...n-1] data samples */
 396
 397   long i,j,o,p;
 398   double y;
 399   double *work=alloca(sizeof(double)*(m+n));
 400
 401   if(!prime)
 402     for(i=0;i<m;i++)
 403       work[i]=0.;
 404   else
 405     for(i=0;i<m;i++)
 406       work[i]=prime[i];
 407
 408   for(i=0;i<n;i++){
 409     y=data[i];
 410     o=i;
 411     p=m;
 412     for(j=0;j<m;j++)
 413       y-=work[o++]*coeff[--p];
 414
 415     data[i]=work[o]=y;
 416   }
 417 }
 418
 419