1 /********************************************************************
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. *
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/ *
12 ********************************************************************
14 function: LPC low level routines
15 last mod: $Id: lpc.c,v 1.21 2000/06/14 01:38:31 xiphmont Exp $
17 ********************************************************************/
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. */
25 /* Preserved Copyright: *********************************************/
27 /* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
28 Technische Universita"t Berlin
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
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.
45 *********************************************************************/
56 /* Autocorrelation LPC coeff generation algorithm invented by
57 N. Levinson in 1947, modified by J. Durbin in 1959. */
59 /* Input : n elements of time doamin data
60 Output: m lpc coefficients, excitation energy */
62 double vorbis_lpc_from_data(double *data,double *lpc,int n,int m){
63 double *aut=alloca(sizeof(double)*(m+1));
67 /* autocorrelation, p+1 lag coefficients */
72 for(i=j;i<n;i++)d+=data[i]*data[i-j];
76 /* Generate lpc coefficients from autocorr values */
80 memset(lpc,0,m*sizeof(double));
87 /* Sum up this iteration's reflection coefficient; note that in
88 Vorbis we don't save it. If anyone wants to recycle this code
89 and needs reflection coefficients, save the results of 'r' from
92 for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
95 /* Update LPC coefficients and total error */
100 lpc[j]+=r*lpc[i-1-j];
103 if(i%2)lpc[j]+=lpc[j]*r;
108 /* we need the error value to know how big an impulse to hit the
114 /* Input : n element envelope spectral curve
115 Output: m lpc coefficients, excitation energy */
117 double vorbis_lpc_from_curve(double *curve,double *lpc,lpc_lookup *l){
120 double *work=alloca(sizeof(double)*(n+n));
124 /* input is a real curve. make it complex-real */
125 /* This mixes phase, but the LPC generation doesn't care. */
127 work[i*2]=curve[i]*fscale;
130 work[n*2-1]=curve[n-1]*fscale;
133 drft_backward(&l->fft,work);
135 /* The autocorrelation will not be circular. Shift, else we lose
136 most of the power in the edges. */
138 for(i=0,j=n/2;i<n/2;){
144 return(vorbis_lpc_from_data(work,lpc,n,m));
147 void lpc_init(lpc_lookup *l,long mapped, int m){
148 memset(l,0,sizeof(lpc_lookup));
153 /* we cheat decoding the LPC spectrum via FFTs */
154 drft_init(&l->fft,mapped*2);
158 void lpc_clear(lpc_lookup *l){
164 /* One can do this the long way by generating the transfer function in
165 the time domain and taking the forward FFT of the result. The
166 results from direct calculation are cleaner and faster.
168 This version does a linear curve generation and then later
169 interpolates the log curve from the linear curve. */
171 void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,
174 memset(curve,0,sizeof(double)*l->ln*2);
178 curve[i*2+1]=lpc[i]/(4*amp);
179 curve[i*2+2]=-lpc[i]/(4*amp);
182 drft_backward(&l->fft,curve); /* reappropriated ;-) */
187 curve[0]=(1./(curve[0]*2+unit));
188 for(i=1;i<l->ln;i++){
189 double real=(curve[i]+curve[l2-i]);
190 double imag=(curve[i]-curve[l2-i]);
192 double a = real + unit;
193 curve[i] = 1.0 / FAST_HYPOT(a, imag);
198 /* subtract or add an lpc filter to data. Vorbis doesn't actually use this. */
200 void vorbis_lpc_residue(double *coeff,double *prime,int m,
201 double *data,long n){
203 /* in: coeff[0...m-1] LPC coefficients
204 prime[0...m-1] initial values
205 data[0...n-1] data samples
206 out: data[0...n-1] residuals from LPC prediction */
209 double *work=alloca(sizeof(double)*(m+n));
222 y-=work[i+j]*coeff[m-j-1];
230 void vorbis_lpc_predict(double *coeff,double *prime,int m,
231 double *data,long n){
233 /* in: coeff[0...m-1] LPC coefficients
234 prime[0...m-1] initial values (allocated size of n+m-1)
235 data[0...n-1] residuals from LPC prediction
236 out: data[0...n-1] data samples */
240 double *work=alloca(sizeof(double)*(m+n));
254 y-=work[o++]*coeff[--p];