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: Aug 02 1999
  18
  19  ********************************************************************/
  20
  21 /* Some of these routines (autocorrelator, LPC coefficient estimator)
  22    are directly derived from and/or modified from code written by
  23    Jutta Degener and Carsten Bormann; thus we include their copyright
  24    below.  The entirety of this file is freely redistributable on the
  25    condition that both of these copyright notices are preserved
  26    without modification.  */
  27
  28 /* Preserved Copyright: *********************************************/
  29
  30 /* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
  31 Technische Universita"t Berlin
  32
  33 Any use of this software is permitted provided that this notice is not
  34 removed and that neither the authors nor the Technische Universita"t
  35 Berlin are deemed to have made any representations as to the
  36 suitability of this software for any purpose nor are held responsible
  37 for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
  38 THIS SOFTWARE.
  39
  40 As a matter of courtesy, the authors request to be informed about uses
  41 this software has found, about bugs in this software, and about any
  42 improvements that may be of general interest.
  43
  44 Berlin, 28.11.1994
  45 Jutta Degener
  46 Carsten Bormann
  47
  48 *********************************************************************/
  49
  50 #include <stdlib.h>
  51 #include <string.h>
  52 #include <math.h>
  53 #include "smallft.h"
  54 #include "lpc.h"
  55
  56 /* This is pared down for Vorbis where we only use LPC to encode
  57    spectral envelope curves.  Thus we only are interested in
  58    generating the coefficients and recovering the curve from the
  59    coefficients.  Autocorrelation LPC coeff generation algorithm
  60    invented by N. Levinson in 1947, modified by J. Durbin in 1959. */
  61
  62 /*  Input : n element envelope curve
  63     Output: m lpc coefficients, excitation energy */
  64
  65 double vorbis_curve_to_lpc(double *curve,int n,double *lpc,int m){
  66   double aut[m+1],work[n+n],error;
  67   drft_lookup dl;
  68   int i,j;
  69
  70   /* input is a real curve. make it complex-real */
  71   for(i=0;i<n;i++){
  72     work[i*2]=curve[i];
  73     work[i*2+1]=0;
  74   }
  75
  76   n*=2;
  77   drft_init(&dl,n);
  78   drft_backward(&dl,work);
  79   drft_clear(&dl);
  80
  81   /* The autocorrelation will not be circular.  Shift, else we lose
  82      most of the power in the edges. */
  83
  84   for(i=0,j=n/2;i<n/2;){
  85     double temp=work[i];
  86     work[i++]=work[j];
  87     work[j++]=temp;
  88   }
  89
  90   /* autocorrelation, p+1 lag coefficients */
  91
  92   j=m+1;
  93   while(j--){
  94     double d=0;
  95     for(i=j;i<n;i++)d+=work[i]*work[i-j];
  96     aut[j]=d;
  97   }
  98
  99   /* Generate lpc coefficients from autocorr values */
 100
 101   error=aut[0];
 102   if(error==0){
 103     memset(lpc,0,m*sizeof(double));
 104     return 0;
 105   }
 106
 107   for(i=0;i<m;i++){
 108     double r=-aut[i+1];
 109
 110     /* Sum up this iteration's reflection coefficient; note that in
 111        Vorbis we don't save it.  If anyone wants to recycle this code
 112        and needs reflection coefficients, save the results of 'r' from
 113        each iteration. */
 114
 115     for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
 116     r/=error;
 117
 118     /* Update LPC coefficients and total error */
 119
 120     lpc[i]=r;
 121     for(j=0;j<i/2;j++){
 122       double tmp=lpc[j];
 123       lpc[j]+=r*lpc[i-1-j];
 124       lpc[i-1-j]+=r*tmp;
 125     }
 126     if(i%2)lpc[j]+=lpc[j]*r;
 127
 128     error*=1.0-r*r;
 129   }
 130
 131   /* we need the error value to know how big an impulse to hit the
 132      filter with later */
 133
 134   return error;
 135 }
 136
 137 /* One can do this the long way by generating the transfer function in
 138    the time domain and taking the forward FFT of the result.  The
 139    results from direct calculation are cleaner and faster, however */
 140
 141 static double vorbis_lpc_magnitude(double w,double *lpc, int m){
 142   int k;
 143   double real=1,imag=0;
 144   for(k=0;k<m;k++){
 145     real+=lpc[k]*cos((k+1)*w);
 146     imag+=lpc[k]*sin((k+1)*w);
 147   }
 148   return(1./sqrt(real*real+imag*imag));
 149 }
 150
 151 /* generate the whole freq response curve on an LPC IIR filter */
 152
 153 void vorbis_lpc_to_curve(double *curve,int n,double *lpc, double amp,int m){
 154   int i;
 155   double w=1./n*M_PI;
 156   for(i=0;i<n;i++)
 157     curve[i]=vorbis_lpc_magnitude(i*w,lpc,m)*amp;
 158 }
 159
 160 /* find frequency response of LPC filter only at nonsero residue
 161    points and apply the envelope to the residue */
 162
 163 void vorbis_lpc_apply(double *residue,int n,double *lpc, double amp,int m){
 164   int i;
 165   double w=1./n*M_PI;
 166   for(i=0;i<n;i++)
 167     if(residue[i])
 168       residue[i]*=vorbis_lpc_magnitude(i*w,lpc,m)*amp;
 169 }
 170
 171