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-1999             *
   9  * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company       *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14   function: LPC low level routines
  15   author: Monty <monty@xiph.org>
  16   modifications by: Monty
  17   last modification date: Nov 16 1999
  18
  19  ********************************************************************/
  20
  21 /* Some of these routines (autocorrelator, LPC coefficient estimator)
  22    are derived from code written by Jutta Degener and Carsten Bormann;
  23    thus we include their copyright below.  The entirety of this file
  24    is freely redistributable on the condition that both of these
  25    copyright notices are preserved without modification.  */
  26
  27 /* Preserved Copyright: *********************************************/
  28
  29 /* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
  30 Technische Universita"t Berlin
  31
  32 Any use of this software is permitted provided that this notice is not
  33 removed and that neither the authors nor the Technische Universita"t
  34 Berlin are deemed to have made any representations as to the
  35 suitability of this software for any purpose nor are held responsible
  36 for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
  37 THIS SOFTWARE.
  38
  39 As a matter of courtesy, the authors request to be informed about uses
  40 this software has found, about bugs in this software, and about any
  41 improvements that may be of general interest.
  42
  43 Berlin, 28.11.1994
  44 Jutta Degener
  45 Carsten Bormann
  46
  47 *********************************************************************/
  48
  49 #include <stdlib.h>
  50 #include <stdio.h>
  51 #include <string.h>
  52 #include <math.h>
  53 #include "os.h"
  54 #include "smallft.h"
  55 #include "lpc.h"
  56 #include "xlogmap.h"
  57
  58 /* This is pared down for Vorbis. Autocorrelation LPC coeff generation
  59    algorithm invented by N. Levinson in 1947, modified by J. Durbin in
  60    1959. */
  61
  62 /* Input : n elements of time doamin data
  63    Output: m lpc coefficients, excitation energy */
  64
  65
  66 double vorbis_lpc_from_data(double *data,double *lpc,int n,int m){
  67   double *aut=alloca(sizeof(double)*(m+1));
  68   double error;
  69   int i,j;
  70
  71   /* autocorrelation, p+1 lag coefficients */
  72
  73   j=m+1;
  74   while(j--){
  75     double d=0;
  76     for(i=j;i<n;i++)d+=data[i]*data[i-j];
  77     aut[j]=d;
  78   }
  79
  80   /* Generate lpc coefficients from autocorr values */
  81
  82   error=aut[0];
  83   if(error==0){
  84     memset(lpc,0,m*sizeof(double));
  85     return 0;
  86   }
  87
  88   for(i=0;i<m;i++){
  89     double r=-aut[i+1];
  90
  91     /* Sum up this iteration's reflection coefficient; note that in
  92        Vorbis we don't save it.  If anyone wants to recycle this code
  93        and needs reflection coefficients, save the results of 'r' from
  94        each iteration. */
  95
  96     for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
  97     r/=error;
  98
  99     /* Update LPC coefficients and total error */
 100
 101     lpc[i]=r;
 102     for(j=0;j<i/2;j++){
 103       double tmp=lpc[j];
 104       lpc[j]+=r*lpc[i-1-j];
 105       lpc[i-1-j]+=r*tmp;
 106     }
 107     if(i%2)lpc[j]+=lpc[j]*r;
 108
 109     error*=1.0-r*r;
 110   }
 111
 112   /* we need the error value to know how big an impulse to hit the
 113      filter with later */
 114
 115   return error;
 116 }
 117
 118 /* Input : n element envelope spectral curve
 119    Output: m lpc coefficients, excitation energy */
 120
 121 double vorbis_lpc_from_spectrum(double *curve,double *lpc,lpc_lookup *l){
 122   int n=l->ln;
 123   int m=l->m;
 124   double *work=alloca(sizeof(double)*(n+n));
 125   double fscale=.5/n;
 126   int i,j;
 127
 128   /* input is a real curve. make it complex-real */
 129   /* This mixes phase, but the LPC generation doesn't care. */
 130   for(i=0;i<n;i++){
 131     work[i*2]=curve[i]*fscale;
 132     work[i*2+1]=0;
 133   }
 134
 135   n*=2;
 136   drft_backward(&l->fft,work);
 137
 138   /* The autocorrelation will not be circular.  Shift, else we lose
 139      most of the power in the edges. */
 140
 141   for(i=0,j=n/2;i<n/2;){
 142     double temp=work[i];
 143     work[i++]=work[j];
 144     work[j++]=temp;
 145   }
 146
 147   return(vorbis_lpc_from_data(work,lpc,n,m));
 148 }
 149
 150 /* On top of this basic LPC infrastructure we introduce two modifications:
 151
 152    1) Filter generation is limited in the resolution of features it
 153    can represent (this is more obvious when the filter is looked at as
 154    a set of LSP coefficients).  Human perception of the audio spectrum
 155    is logarithmic not only in amplitude, but also frequency.  Because
 156    the high frequency features we'll need to encode will be broader
 157    than the low frequency features, filter generation will be
 158    dominated by higher frequencies (when most of the energy is in the
 159    lowest frequencies, and greatest perceived resolution is in the
 160    midrange).  To avoid this effect, Vorbis encodes the frequency
 161    spectrum with a biased log frequency scale. The intent is to
 162    roughly equalize the sizes of the octaves (see xlogmap.h)
 163
 164    2) When we change the frequency scale, we also change the
 165    (apparent) relative energies of the bands; that is, on a log scale
 166    covering 5 octaves, the highest octave goes from being represented
 167    in half the bins, to only 1/32 of the bins.  If the amplitudes
 168    remain the same, we have divided the energy represented by the
 169    highest octave by 16 (as far as Levinson-Durbin is concerned).
 170    This will seriously skew filter generation, which bases calculation
 171    on the mean square error with respect to energy.  Thus, Vorbis
 172    normalizes the amplitudes of the log spectrum frequencies to keep
 173    the relative octave energies correct. */
 174
 175 /* n == size of vector to be used for filter, m == order of filter,
 176    oct == octaves in normalized scale, encode_p == encode (1) or
 177    decode (0) */
 178
 179 void lpc_init(lpc_lookup *l,int n, int mapped, int m, int oct, int encode_p){
 180   double bias=LOG_BIAS(n,oct);
 181   double scale=(float)mapped/(float)oct; /* where n==mapped */
 182   int i;
 183
 184   memset(l,0,sizeof(lpc_lookup));
 185
 186   l->n=n;
 187   l->ln=mapped;
 188   l->m=m;
 189   l->iscale=malloc(n*sizeof(int));
 190   l->ifrac=malloc(n*sizeof(double));
 191   l->norm=malloc(n*sizeof(double));
 192
 193   for(i=0;i<n;i++){
 194     /* how much 'real estate' in the log domain does the bin in the
 195        linear domain represent? */
 196     double logA=LOG_X(i,bias);
 197     double logB=LOG_X(i+1.,bias);
 198     l->norm[i]=logB-logA;  /* this much */
 199   }
 200
 201   /* the scale is encode/decode specific for algebraic simplicity */
 202
 203   if(encode_p){
 204     /* encode */
 205     l->escale=malloc(mapped*sizeof(double));
 206     l->uscale=malloc(n*sizeof(int));
 207
 208     /* undersample guard */
 209     for(i=0;i<n;i++){
 210       l->uscale[i]=rint(LOG_X(i,bias)/oct*mapped);
 211     }
 212
 213     for(i=0;i<mapped;i++){
 214       l->escale[i]=LINEAR_X(i/scale,bias);
 215       l->uscale[(int)(floor(l->escale[i]))]=-1;
 216       l->uscale[(int)(ceil(l->escale[i]))]=-1;
 217     }
 218
 219
 220   }
 221   /* decode; encode may use this too */
 222
 223   drft_init(&l->fft,mapped*2);
 224   for(i=0;i<n;i++){
 225     double is=LOG_X(i,bias)/oct*mapped;
 226     if(is<0.)is=0.;
 227
 228     l->iscale[i]=floor(is);
 229     if(l->iscale[i]>=l->ln-1)l->iscale[i]=l->ln-2;
 230
 231     l->ifrac[i]=is-floor(is);
 232     if(l->ifrac[i]>1.)l->ifrac[i]=1.;
 233
 234   }
 235 }
 236
 237 void lpc_clear(lpc_lookup *l){
 238   if(l){
 239     if(l->escale)free(l->escale);
 240     drft_clear(&l->fft);
 241     free(l->iscale);
 242     free(l->ifrac);
 243     free(l->norm);
 244   }
 245 }
 246
 247
 248 /* less efficient than the decode side (written for clarity).  We're
 249    not bottlenecked here anyway */
 250
 251 double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
 252   /* map the input curve to a log curve for encoding */
 253
 254   /* for clarity, mapped and n are both represented although setting
 255      'em equal is a decent rule of thumb. The below must be reworked
 256      slightly if mapped != n */
 257
 258   int mapped=l->ln;
 259   double *work=alloca(sizeof(double)*mapped);
 260   int i;
 261
 262   /* fairly correct for low frequencies, naieve for high frequencies
 263      (suffers from undersampling) */
 264   for(i=0;i<mapped;i++){
 265     double lin=l->escale[i];
 266     int a=floor(lin);
 267     int b=ceil(lin);
 268     double del=lin-floor(lin);
 269
 270     work[i]=(curve[a]/l->norm[a]*(1.-del)+
 271              curve[b]/l->norm[b]*del);
 272
 273   }
 274
 275   /*  for(i=0;i<l->n;i++)
 276     if(l->uscale[i]>0)
 277     if(work[l->uscale[i]]<curve[i])work[l->uscale[i]]=curve[i];*/
 278
 279 #ifdef ANALYSIS
 280   {
 281     int j;
 282     FILE *out;
 283     char buffer[80];
 284     static int frameno=0;
 285
 286     sprintf(buffer,"preloglpc%d.m",frameno++);
 287     out=fopen(buffer,"w+");
 288     for(j=0;j<l->ln;j++)
 289       fprintf(out,"%g\n",work[j]);
 290     fclose(out);
 291   }
 292 #endif
 293
 294   return vorbis_lpc_from_spectrum(work,lpc,l);
 295 }
 296
 297
 298 /* One can do this the long way by generating the transfer function in
 299    the time domain and taking the forward FFT of the result.  The
 300    results from direct calculation are cleaner and faster.
 301
 302    This version does a linear curve generation and then later
 303    interpolates the log curve from the linear curve.  This could stand
 304    optimization; it could both be more precise as well as not compute
 305    quite a few unused values */
 306
 307 void _vlpc_de_helper(double *curve,double *lpc,double amp,
 308                             lpc_lookup *l){
 309   int i;
 310   memset(curve,0,sizeof(double)*l->ln*2);
 311   if(amp==0)return;
 312
 313   for(i=0;i<l->m;i++){
 314     curve[i*2+1]=lpc[i]/4/amp;
 315     curve[i*2+2]=-lpc[i]/4/amp;
 316   }
 317
 318   drft_backward(&l->fft,curve); /* reappropriated ;-) */
 319
 320   {
 321     int l2=l->ln*2;
 322     double unit=1./amp;
 323     curve[0]=(1./(curve[0]*2+unit));
 324     for(i=1;i<l->ln;i++){
 325       double real=(curve[i]+curve[l2-i]);
 326       double imag=(curve[i]-curve[l2-i]);
 327       curve[i]=(1./hypot(real+unit,imag));
 328     }
 329   }
 330 }
 331
 332
 333 /* generate the whole freq response curve on an LPC IIR filter */
 334
 335 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
 336   double *lcurve=alloca(sizeof(double)*(l->ln*2));
 337   int i;
 338
 339   _vlpc_de_helper(lcurve,lpc,amp,l);
 340
 341 #ifdef ANALYSIS
 342   {
 343     int j;
 344     FILE *out;
 345     char buffer[80];
 346     static int frameno=0;
 347
 348     sprintf(buffer,"loglpc%d.m",frameno++);
 349     out=fopen(buffer,"w+");
 350     for(j=0;j<l->ln;j++)
 351       fprintf(out,"%g\n",lcurve[j]);
 352     fclose(out);
 353   }
 354 #endif
 355
 356   if(amp==0)return;
 357
 358   for(i=0;i<l->n;i++){
 359     int ii=l->iscale[i];
 360     curve[i]=((1.-l->ifrac[i])*lcurve[ii]+
 361               l->ifrac[i]*lcurve[ii+1])*l->norm[i];
 362   }
 363
 364 }
 365
 366 void vorbis_lpc_apply(double *residue,double *lpc,double amp,lpc_lookup *l){
 367   double *lcurve=alloca(sizeof(double)*((l->ln+l->n)*2));
 368   int i;
 369
 370   if(amp==0){
 371     memset(residue,0,l->n*sizeof(double));
 372   }else{
 373
 374     _vlpc_de_helper(lcurve,lpc,amp,l);
 375
 376     for(i=0;i<l->n;i++){
 377       if(residue[i]!=0){
 378         int ii=l->iscale[i];
 379         residue[i]*=((1.-l->ifrac[i])*lcurve[ii]+
 380                      l->ifrac[i]*lcurve[ii+1])*l->norm[i];
 381       }
 382     }
 383   }
 384 }
 385
 386 /* subtract or add an lpc filter to data.  Vorbis doesn't actually use this. */
 387
 388 void vorbis_lpc_residue(double *coeff,double *prime,int m,
 389                         double *data,long n){
 390
 391   /* in: coeff[0...m-1] LPC coefficients
 392          prime[0...m-1] initial values
 393          data[0...n-1] data samples
 394     out: data[0...n-1] residuals from LPC prediction */
 395
 396   long i,j;
 397   double *work=alloca(sizeof(double)*(m+n));
 398   double y;
 399
 400   if(!prime)
 401     for(i=0;i<m;i++)
 402       work[i]=0;
 403   else
 404     for(i=0;i<m;i++)
 405       work[i]=prime[i];
 406
 407   for(i=0;i<n;i++){
 408     y=0;
 409     for(j=0;j<m;j++)
 410       y-=work[i+j]*coeff[m-j-1];
 411
 412     work[i+m]=data[i];
 413     data[i]-=y;
 414   }
 415 }
 416
 417
 418 void vorbis_lpc_predict(double *coeff,double *prime,int m,
 419                      double *data,long n){
 420
 421   /* in: coeff[0...m-1] LPC coefficients
 422          prime[0...m-1] initial values (allocated size of n+m-1)
 423          data[0...n-1] residuals from LPC prediction
 424     out: data[0...n-1] data samples */
 425
 426   long i,j,o,p;
 427   double y;
 428   double *work=alloca(sizeof(double)*(m+n));
 429
 430   if(!prime)
 431     for(i=0;i<m;i++)
 432       work[i]=0.;
 433   else
 434     for(i=0;i<m;i++)
 435       work[i]=prime[i];
 436
 437   for(i=0;i<n;i++){
 438     y=data[i];
 439     o=i;
 440     p=m;
 441     for(j=0;j<m;j++)
 442       y-=work[o++]*coeff[--p];
 443
 444     data[i]=work[o]=y;
 445   }
 446 }
 447
 448