diff options
author | Monty <xiphmont@xiph.org> | 2000-03-29 03:49:29 +0000 |
---|---|---|
committer | Monty <xiphmont@xiph.org> | 2000-03-29 03:49:29 +0000 |
commit | 0695c7cbf5d766b7db3c664fa1bb82531c71fa38 (patch) | |
tree | 84e37541e04291d16ed0b09d7294ff15af205534 | |
parent | 2c3a9abf067db21aca6235c7f32c33e7f25e838f (diff) | |
download | libvorbis-git-unlabeled-1.18.2.tar.gz |
Don't want to lose anything while I'm integrating (also don;t want tounlabeled-1.18.2
disturb mainline till I'm done)
Monty
svn path=/branches/unlabeled-1.18.2/vorbis/; revision=286
-rw-r--r-- | lib/lpc.c | 375 |
1 files changed, 375 insertions, 0 deletions
diff --git a/lib/lpc.c b/lib/lpc.c new file mode 100644 index 00000000..ced6b298 --- /dev/null +++ b/lib/lpc.c @@ -0,0 +1,375 @@ +/******************************************************************** + * * + * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. * + * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY * + * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. * + * PLEASE READ THESE TERMS DISTRIBUTING. * + * * + * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000 * + * by Monty <monty@xiph.org> and The XIPHOPHORUS Company * + * http://www.xiph.org/ * + * * + ******************************************************************** + + function: LPC low level routines + last mod: $Id: lpc.c,v 1.18.2.1 2000/03/29 03:49:28 xiphmont Exp $ + + ********************************************************************/ + +/* Some of these routines (autocorrelator, LPC coefficient estimator) + are derived from code written by Jutta Degener and Carsten Bormann; + thus we include their copyright below. The entirety of this file + is freely redistributable on the condition that both of these + copyright notices are preserved without modification. */ + +/* Preserved Copyright: *********************************************/ + +/* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann, +Technische Universita"t Berlin + +Any use of this software is permitted provided that this notice is not +removed and that neither the authors nor the Technische Universita"t +Berlin are deemed to have made any representations as to the +suitability of this software for any purpose nor are held responsible +for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR +THIS SOFTWARE. + +As a matter of courtesy, the authors request to be informed about uses +this software has found, about bugs in this software, and about any +improvements that may be of general interest. + +Berlin, 28.11.1994 +Jutta Degener +Carsten Bormann + +*********************************************************************/ + +#include <stdlib.h> +#include <string.h> +#include <math.h> +#include "os.h" +#include "smallft.h" +#include "lpc.h" +#include "scales.h" +#include "misc.h" + +/* Autocorrelation LPC coeff generation algorithm invented by + N. Levinson in 1947, modified by J. Durbin in 1959. */ + +/* Input : n elements of time doamin data + Output: m lpc coefficients, excitation energy */ + +double vorbis_lpc_from_data(double *data,double *lpc,int n,int m){ + double *aut=alloca(sizeof(double)*(m+1)); + double error; + int i,j; + + /* autocorrelation, p+1 lag coefficients */ + + j=m+1; + while(j--){ + double d=0; + for(i=j;i<n;i++)d+=data[i]*data[i-j]; + aut[j]=d; + } + + /* Generate lpc coefficients from autocorr values */ + + error=aut[0]; + if(error==0){ + memset(lpc,0,m*sizeof(double)); + return 0; + } + + for(i=0;i<m;i++){ + double r=-aut[i+1]; + + /* Sum up this iteration's reflection coefficient; note that in + Vorbis we don't save it. If anyone wants to recycle this code + and needs reflection coefficients, save the results of 'r' from + each iteration. */ + + for(j=0;j<i;j++)r-=lpc[j]*aut[i-j]; + r/=error; + + /* Update LPC coefficients and total error */ + + lpc[i]=r; + for(j=0;j<i/2;j++){ + double tmp=lpc[j]; + lpc[j]+=r*lpc[i-1-j]; + lpc[i-1-j]+=r*tmp; + } + if(i%2)lpc[j]+=lpc[j]*r; + + error*=1.0-r*r; + } + + /* we need the error value to know how big an impulse to hit the + filter with later */ + + return error; +} + +/* Input : n element envelope spectral curve + Output: m lpc coefficients, excitation energy */ + +double vorbis_lpc_from_spectrum(double *curve,double *lpc,lpc_lookup *l){ + int n=l->ln; + int m=l->m; + double *work=alloca(sizeof(double)*(n+n)); + double fscale=.5/n; + int i,j; + + /* input is a real curve. make it complex-real */ + /* This mixes phase, but the LPC generation doesn't care. */ + for(i=0;i<n;i++){ + work[i*2]=curve[i]*fscale; + work[i*2+1]=0; + } + + n*=2; + drft_backward(&l->fft,work); + + /* The autocorrelation will not be circular. Shift, else we lose + most of the power in the edges. */ + + for(i=0,j=n/2;i<n/2;){ + double temp=work[i]; + work[i++]=work[j]; + work[j++]=temp; + } + + return(vorbis_lpc_from_data(work,lpc,n,m)); +} + +/* initialize Bark scale and normalization lookups. We could do this + with static tables, but Vorbis allows a number of possible + combinations, so it's best to do it computationally. + + The below is authoritative in terms of defining scale mapping. + Note that the scale depends on the sampling rate as well as the + linear block and mapping sizes */ + +void lpc_init(lpc_lookup *l,int n, long mapped, long rate, int m){ + int i; + double scale; + memset(l,0,sizeof(lpc_lookup)); + + l->n=n; + l->ln=mapped; + l->m=m; + + l->linearmap=malloc(n*sizeof(int)); + l->barknorm=malloc(mapped*sizeof(double)); + + /* we choose a scaling constant so that: + floor(bark(rate/2-1)*C)=mapped-1 + floor(bark(rate/2)*C)=mapped */ + + scale=mapped/toBARK(rate/2.); + + /* the mapping from a linear scale to a smaller bark scale is + straightforward. We do *not* make sure that the linear mapping + does not skip bark-scale bins; the decoder simply skips them and + the encoder may do what it wishes in filling them. They're + necessary in some mapping combinations to keep the scale spacing + accurate */ + { + int last=-1; + for(i=0;i<n;i++){ + int val=floor( toBARK((rate/2.)/n*i) *scale); /* bark numbers + represent + band edges */ + if(val>=mapped)val=mapped; /* guard against the approximation */ + l->linearmap[i]=val; + last=val; + } + } + + /* 'Normalization' is just making sure that power isn't lost in the + log scale by virtue of compressing the scale in higher + frequencies. We figure the weight of bands in proportion to + their linear/bark width ratio below, again, authoritatively. We + use computed width (not the number of actual bins above) for + smoothness in the scale; they should agree closely */ + + /* keep it 0. to 1., else the dynamic range starts spreading through + all the squaring... */ + + for(i=0;i<mapped;i++) + l->barknorm[i]=(fromBARK((i+1)/scale)-fromBARK(i/scale)); + for(i=0;i<mapped;i++) + l->barknorm[i]/=l->barknorm[mapped-1]; + + /* we cheat decoding the LPC spectrum via FFTs */ + + drft_init(&l->fft,mapped*2); + +} + +void lpc_clear(lpc_lookup *l){ + if(l){ + if(l->barknorm)free(l->barknorm); + if(l->linearmap)free(l->linearmap); + drft_clear(&l->fft); + } +} + + +/* less efficient than the decode side (written for clarity). We're + not bottlenecked here anyway */ +double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){ + /* map the input curve to a bark-scale curve for encoding */ + + int mapped=l->ln; + double *work=alloca(sizeof(double)*mapped); + int i,j,last=0; + + memset(work,0,sizeof(double)*mapped); + + /* Only the decode side is behavior-specced; for now in the encoder, + we select the maximum value of each band as representative (this + helps make sure peaks don't go out of range. In error terms, + selecting min would make more sense, but the codebook is trained + numerically, so we don't actually lose. We'd still want to + use the original curve for error and noise estimation */ + + for(i=0;i<l->n;i++){ + int bark=l->linearmap[i]; + if(work[bark]<curve[i])work[bark]=curve[i]; + if(bark>last+1){ + /* If the bark scale is climbing rapidly, some bins may end up + going unused. This isn't a waste actually; it keeps the + scale resolution even so that the LPC generator has an easy + time. However, if we leave the bins empty we lose energy. + So, fill 'em in. The decoder does not do anything with he + unused bins, so we can fill them anyway we like to end up + with a better spectral curve */ + + /* we'll always have a bin zero, so we don't need to guard init */ + long span=bark-last; + for(j=1;j<span;j++){ + double del=(double)j/span; + work[j+last]=work[bark]*del+work[last]*(1.-del); + } + } + last=bark; + } + /*for(i=0;i<mapped;i++)work[i]*=l->barknorm[i];*/ + + return vorbis_lpc_from_spectrum(work,lpc,l); +} + + +/* One can do this the long way by generating the transfer function in + the time domain and taking the forward FFT of the result. The + results from direct calculation are cleaner and faster. + + This version does a linear curve generation and then later + interpolates the log curve from the linear curve. */ + +void _vlpc_de_helper(double *curve,double *lpc,double amp, + lpc_lookup *l){ + int i; + memset(curve,0,sizeof(double)*l->ln*2); + if(amp==0)return; + + for(i=0;i<l->m;i++){ + curve[i*2+1]=lpc[i]/(4*amp); + curve[i*2+2]=-lpc[i]/(4*amp); + } + + drft_backward(&l->fft,curve); /* reappropriated ;-) */ + + { + int l2=l->ln*2; + double unit=1./amp; + curve[0]=(1./(curve[0]*2+unit)); + for(i=1;i<l->ln;i++){ + double real=(curve[i]+curve[l2-i]); + double imag=(curve[i]-curve[l2-i]); + curve[i]=(1./hypot(real+unit,imag)); + } + } +} + +/* generate the whole freq response curve of an LPC IIR filter */ + +void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){ + double *lcurve=alloca(sizeof(double)*(l->ln*2)); + int i; + + if(amp==0){ + memset(curve,0,sizeof(double)*l->n); + return; + } + _vlpc_de_helper(lcurve,lpc,amp,l); + + /*for(i=0;i<l->ln;i++)lcurve[i]/=l->barknorm[i];*/ + for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]]; + +} + +/* subtract or add an lpc filter to data. Vorbis doesn't actually use this. */ + +void vorbis_lpc_residue(double *coeff,double *prime,int m, + double *data,long n){ + + /* in: coeff[0...m-1] LPC coefficients + prime[0...m-1] initial values + data[0...n-1] data samples + out: data[0...n-1] residuals from LPC prediction */ + + long i,j; + double *work=alloca(sizeof(double)*(m+n)); + double y; + + if(!prime) + for(i=0;i<m;i++) + work[i]=0; + else + for(i=0;i<m;i++) + work[i]=prime[i]; + + for(i=0;i<n;i++){ + y=0; + for(j=0;j<m;j++) + y-=work[i+j]*coeff[m-j-1]; + + work[i+m]=data[i]; + data[i]-=y; + } +} + + +void vorbis_lpc_predict(double *coeff,double *prime,int m, + double *data,long n){ + + /* in: coeff[0...m-1] LPC coefficients + prime[0...m-1] initial values (allocated size of n+m-1) + data[0...n-1] residuals from LPC prediction + out: data[0...n-1] data samples */ + + long i,j,o,p; + double y; + double *work=alloca(sizeof(double)*(m+n)); + + if(!prime) + for(i=0;i<m;i++) + work[i]=0.; + else + for(i=0;i<m;i++) + work[i]=prime[i]; + + for(i=0;i<n;i++){ + y=data[i]; + o=i; + p=m; + for(j=0;j<m;j++) + y-=work[o++]*coeff[--p]; + + data[i]=work[o]=y; + } +} + |