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