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.11 1999/12/31 12:35:14 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 "barkmel.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 (note that for a given sample rate
 153    and block size, there's generally a fairly obviously optimal
 154    mapping size */
 155
 156 void lpc_init(lpc_lookup *l,int n, long mapped, long rate, int m){
 157   int i;
 158   double scale;
 159   memset(l,0,sizeof(lpc_lookup));
 160
 161   l->n=n;
 162   l->ln=mapped;
 163   l->m=m;
 164
 165   l->linearmap=malloc(n*sizeof(int));
 166   l->barknorm=malloc(mapped*sizeof(double));
 167
 168   /* we choose a scaling constant so that:
 169      floor(bark(rate-1)*C)=mapped-1
 170      floor(bark(rate)*C)=mapped */
 171
 172   scale=mapped/fBARK(rate);
 173
 174   /* the mapping from a linear scale to a smaller bark scale is
 175      straightforward with a single catch; make sure not to skip any
 176      bark-scale bins.  In order to do this, we assign map_N = min
 177      (map_N-1 + 1, bark(N)) */
 178   {
 179     int last=-1;
 180     for(i=0;i<n;i++){
 181       int val=floor( fBARK(((double)rate)/n*i) *scale); /* bark numbers
 182                                                           represent
 183                                                           band edges */
 184       if(val>=mapped)val=mapped; /* guard against the approximation */
 185       if(val>last+1)val=last+1;
 186       l->linearmap[i]=val;
 187       last=val;
 188     }
 189   }
 190
 191   /* 'Normalization' is just making sure that power isn't lost in the
 192      log scale by virtue of compressing the scale in higher
 193      frequencies.  We figure the weight of bands in proportion to
 194      their linear/bark width ratio below, again, authoritatively.  We
 195      use computed width (not the number of actual bins above) for
 196      smoothness in the scale; they should agree closely unless the
 197      encoder chose parameters poorly (and got a bark scale that would
 198      have had lots of skipped bins) */
 199
 200   for(i=0;i<mapped;i++)
 201     l->barknorm[i]=iBARK((i+1)/scale)-iBARK(i/scale);
 202
 203   /* we cheat decoding the LPC spectrum via FFTs */
 204
 205   drft_init(&l->fft,mapped*2);
 206
 207 }
 208
 209 void lpc_clear(lpc_lookup *l){
 210   if(l){
 211     if(l->barknorm)free(l->barknorm);
 212     if(l->linearmap)free(l->linearmap);
 213     drft_clear(&l->fft);
 214   }
 215 }
 216
 217
 218 /* less efficient than the decode side (written for clarity).  We're
 219    not bottlenecked here anyway */
 220
 221 double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
 222   /* map the input curve to a bark-scale curve for encoding */
 223
 224   int mapped=l->ln;
 225   double *work=alloca(sizeof(double)*mapped);
 226   int i;
 227
 228   memset(work,0,sizeof(double)*mapped);
 229
 230   /* Only the decode side is behavior-specced; for now in the encoder,
 231      we select the maximum value of each band as representative (this
 232      helps make sure peaks don't go out of range.  In error terms,
 233      selecting min would make more sense, but the codebook is trained
 234      numerically, so we don't lose in encoding.  We'd still want to
 235      use the original curve for error and noise estimation */
 236
 237   for(i=0;i<l->n;i++){
 238     int bark=l->linearmap[i];
 239     if(work[bark]<curve[i])work[bark]=curve[i];
 240   }
 241   for(i=0;i<mapped;i++)work[i]*=l->barknorm[i];
 242
 243 #ifdef ANALYSIS
 244   {
 245     int j;
 246     FILE *out;
 247     char buffer[80];
 248     static int frameno=0;
 249
 250     sprintf(buffer,"prelpc%d.m",frameno);
 251     out=fopen(buffer,"w+");
 252     for(j=0;j<l->n;j++)
 253       fprintf(out,"%g\n",curve[j]);
 254     fclose(out);
 255     sprintf(buffer,"preloglpc%d.m",frameno++);
 256     out=fopen(buffer,"w+");
 257     for(j=0;j<l->ln;j++)
 258       fprintf(out,"%g\n",work[j]);
 259     fclose(out);
 260   }
 261 #endif
 262
 263   return vorbis_lpc_from_spectrum(work,lpc,l);
 264 }
 265
 266
 267 /* One can do this the long way by generating the transfer function in
 268    the time domain and taking the forward FFT of the result.  The
 269    results from direct calculation are cleaner and faster.
 270
 271    This version does a linear curve generation and then later
 272    interpolates the log curve from the linear curve.  This could stand
 273    optimization; it could both be more precise as well as not compute
 274    quite a few unused values */
 275
 276 void _vlpc_de_helper(double *curve,double *lpc,double amp,
 277                             lpc_lookup *l){
 278   int i;
 279   memset(curve,0,sizeof(double)*l->ln*2);
 280   if(amp==0)return;
 281
 282   for(i=0;i<l->m;i++){
 283     curve[i*2+1]=lpc[i]/4/amp;
 284     curve[i*2+2]=-lpc[i]/4/amp;
 285   }
 286
 287   drft_backward(&l->fft,curve); /* reappropriated ;-) */
 288
 289   {
 290     int l2=l->ln*2;
 291     double unit=1./amp;
 292     curve[0]=(1./(curve[0]*2+unit));
 293     for(i=1;i<l->ln;i++){
 294       double real=(curve[i]+curve[l2-i]);
 295       double imag=(curve[i]-curve[l2-i]);
 296       curve[i]=(1./hypot(real+unit,imag));
 297     }
 298   }
 299 }
 300
 301
 302 /* generate the whole freq response curve of an LPC IIR filter */
 303
 304 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
 305   double *lcurve=alloca(sizeof(double)*(l->ln*2));
 306   int i;
 307   static int frameno=0;
 308
 309   _vlpc_de_helper(lcurve,lpc,amp,l);
 310
 311 #ifdef ANALYSIS
 312   {
 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   if(amp==0)return;
 326
 327   for(i=0;i<l->ln;i++)lcurve[i]/=l->barknorm[i];
 328   for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
 329
 330 #ifdef ANALYSIS
 331   {
 332     int j;
 333     FILE *out;
 334     char buffer[80];
 335
 336     sprintf(buffer,"lpc%d.m",frameno-1);
 337     out=fopen(buffer,"w+");
 338     for(j=0;j<l->n;j++)
 339       fprintf(out,"%g\n",curve[j]);
 340     fclose(out);
 341   }
 342 #endif
 343 }
 344
 345 void vorbis_lpc_apply(double *residue,double *lpc,double amp,lpc_lookup *l){
 346   double *lcurve=alloca(sizeof(double)*((l->ln+l->n)*2));
 347   int i;
 348   static int frameno=0;
 349
 350   if(amp==0){
 351     memset(residue,0,l->n*sizeof(double));
 352   }else{
 353
 354     _vlpc_de_helper(lcurve,lpc,amp,l);
 355
 356 #ifdef ANALYSIS
 357   {
 358     int j;
 359     FILE *out;
 360     char buffer[80];
 361
 362     sprintf(buffer,"loglpc%d.m",frameno++);
 363     out=fopen(buffer,"w+");
 364     for(j=0;j<l->ln;j++)
 365       fprintf(out,"%g\n",lcurve[j]);
 366     fclose(out);
 367   }
 368 #endif
 369
 370     for(i=0;i<l->ln;i++)lcurve[i]/=l->barknorm[i];
 371     for(i=0;i<l->n;i++)
 372       if(residue[i]!=0)
 373         residue[i]*=lcurve[l->linearmap[i]];
 374   }
 375 }
 376
 377 /* subtract or add an lpc filter to data.  Vorbis doesn't actually use this. */
 378
 379 void vorbis_lpc_residue(double *coeff,double *prime,int m,
 380                         double *data,long n){
 381
 382   /* in: coeff[0...m-1] LPC coefficients
 383          prime[0...m-1] initial values
 384          data[0...n-1] data samples
 385     out: data[0...n-1] residuals from LPC prediction */
 386
 387   long i,j;
 388   double *work=alloca(sizeof(double)*(m+n));
 389   double y;
 390
 391   if(!prime)
 392     for(i=0;i<m;i++)
 393       work[i]=0;
 394   else
 395     for(i=0;i<m;i++)
 396       work[i]=prime[i];
 397
 398   for(i=0;i<n;i++){
 399     y=0;
 400     for(j=0;j<m;j++)
 401       y-=work[i+j]*coeff[m-j-1];
 402
 403     work[i+m]=data[i];
 404     data[i]-=y;
 405   }
 406 }
 407
 408
 409 void vorbis_lpc_predict(double *coeff,double *prime,int m,
 410                      double *data,long n){
 411
 412   /* in: coeff[0...m-1] LPC coefficients
 413          prime[0...m-1] initial values (allocated size of n+m-1)
 414          data[0...n-1] residuals from LPC prediction
 415     out: data[0...n-1] data samples */
 416
 417   long i,j,o,p;
 418   double y;
 419   double *work=alloca(sizeof(double)*(m+n));
 420
 421   if(!prime)
 422     for(i=0;i<m;i++)
 423       work[i]=0.;
 424   else
 425     for(i=0;i<m;i++)
 426       work[i]=prime[i];
 427
 428   for(i=0;i<n;i++){
 429     y=data[i];
 430     o=i;
 431     p=m;
 432     for(j=0;j<m;j++)
 433       y-=work[o++]*coeff[--p];
 434
 435     data[i]=work[o]=y;
 436   }
 437 }
 438
 439