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   int i,j;
  68
  69   /* input is a real curve. make it complex-real */
  70   for(i=0;i<n;i++){
  71     work[i*2]=curve[i];
  72     work[i*2+1]=0;
  73   }
  74
  75   n*=2;
  76   fft_backward(n,work,NULL,NULL);
  77
  78   /* The autocorrelation will not be circular.  Shift, else we lose
  79      most of the power in the edges. */
  80
  81   for(i=0,j=n/2;i<n/2;){
  82     double temp=work[i];
  83     work[i++]=work[j];
  84     work[j++]=temp;
  85   }
  86
  87   /* autocorrelation, p+1 lag coefficients */
  88
  89   j=m+1;
  90   while(j--){
  91     double d=0;
  92     for(i=j;i<n;i++)d+=work[i]*work[i-j];
  93     aut[j]=d;
  94   }
  95
  96   /* Generate lpc coefficients from autocorr values */
  97
  98   error=aut[0];
  99   if(error==0){
 100     memset(lpc,0,m*sizeof(double));
 101     return 0;
 102   }
 103
 104   for(i=0;i<m;i++){
 105     double r=-aut[i+1];
 106
 107     /* Sum up this iteration's reflection coefficient; note that in
 108        Vorbis we don't save it.  If anyone wants to recycle this code
 109        and needs reflection coefficients, save the results of 'r' from
 110        each iteration. */
 111
 112     for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
 113     r/=error;
 114
 115     /* Update LPC coefficients and total error */
 116
 117     lpc[i]=r;
 118     for(j=0;j<i/2;j++){
 119       double tmp=lpc[j];
 120       lpc[j]+=r*lpc[i-1-j];
 121       lpc[i-1-j]+=r*tmp;
 122     }
 123     if(i%2)lpc[j]+=lpc[j]*r;
 124
 125     error*=1.0-r*r;
 126   }
 127
 128   /* we need the error value to know how big an impulse to hit the
 129      filter with later */
 130
 131   return error;
 132 }
 133
 134 /* One can do this the long way by generating the transfer function in
 135    the time domain and taking the forward FFT of the result.  The
 136    results from direct calculation are cleaner and faster, however */
 137
 138 static double vorbis_lpc_magnitude(double w,double *lpc, int m){
 139   int k;
 140   double real=1,imag=0;
 141   for(k=0;k<m;k++){
 142     real+=lpc[k]*cos((k+1)*w);
 143     imag+=lpc[k]*sin((k+1)*w);
 144   }
 145   return(1/sqrt(real*real+imag*imag));
 146 }
 147
 148 /* generate the whole freq response curve on an LPC IIR filter */
 149
 150 void vorbis_lpc_to_curve(double *curve,int n,double *lpc, double amp,int m){
 151   int i;
 152   double w=1./n*M_PI;
 153   for(i=0;i<n;i++)
 154     curve[i]=vorbis_lpc_magnitude(i*w,lpc,m)*amp;
 155 }
 156
 157 /* find frequency response of LPC filter only at nonsero residue
 158    points and apply the envelope to the residue */
 159
 160 void vorbis_lpc_apply(double *residue,int n,double *lpc, double amp,int m){
 161   int i;
 162   double w=1./n*M_PI;
 163   for(i=0;i<n;i++)
 164     if(residue[i])
 165       residue[i]*=vorbis_lpc_magnitude(i*w,lpc,m)*amp;
 166 }
 167
 168